efee5258bc
The MPFR library is a C library for multiple-precision floating-point computations with exact rounding (also called correct rounding). It is based on the GMP multiple-precision library and should replace the MPF class in further releases of GMP. GCC >= 4.2 requires MPFR.
230 lines
7.6 KiB
C
230 lines
7.6 KiB
C
/* mpfr_rem1 -- internal function
|
|
mpfr_fmod -- compute the floating-point remainder of x/y
|
|
mpfr_remquo and mpfr_remainder -- argument reduction functions
|
|
|
|
Copyright 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
|
|
Contributed by the Arenaire and Cacao projects, INRIA.
|
|
|
|
This file is part of the GNU MPFR Library.
|
|
|
|
The GNU MPFR Library is free software; you can redistribute it and/or modify
|
|
it under the terms of the GNU Lesser General Public License as published by
|
|
the Free Software Foundation; either version 3 of the License, or (at your
|
|
option) any later version.
|
|
|
|
The GNU MPFR Library is distributed in the hope that it will be useful, but
|
|
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
|
|
License for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public License
|
|
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
|
|
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
|
|
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
|
|
|
|
# include "mpfr-impl.h"
|
|
|
|
/* we return as many bits as we can, keeping just one bit for the sign */
|
|
# define WANTED_BITS (sizeof(long) * CHAR_BIT - 1)
|
|
|
|
/*
|
|
rem1 works as follows:
|
|
The first rounding mode rnd_q indicate if we are actually computing
|
|
a fmod (MPFR_RNDZ) or a remainder/remquo (MPFR_RNDN).
|
|
|
|
Let q = x/y rounded to an integer in the direction rnd_q.
|
|
Put x - q*y in rem, rounded according to rnd.
|
|
If quo is not null, the value stored in *quo has the sign of q,
|
|
and agrees with q with the 2^n low order bits.
|
|
In other words, *quo = q (mod 2^n) and *quo q >= 0.
|
|
If rem is zero, then it has the sign of x.
|
|
The returned 'int' is the inexact flag giving the place of rem wrt x - q*y.
|
|
|
|
If x or y is NaN: *quo is undefined, rem is NaN.
|
|
If x is Inf, whatever y: *quo is undefined, rem is NaN.
|
|
If y is Inf, x not NaN nor Inf: *quo is 0, rem is x.
|
|
If y is 0, whatever x: *quo is undefined, rem is NaN.
|
|
If x is 0, whatever y (not NaN nor 0): *quo is 0, rem is x.
|
|
|
|
Otherwise if x and y are neither NaN, Inf nor 0, q is always defined,
|
|
thus *quo is.
|
|
Since |x - q*y| <= y/2, no overflow is possible.
|
|
Only an underflow is possible when y is very small.
|
|
*/
|
|
|
|
static int
|
|
mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q,
|
|
mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
|
|
{
|
|
mpfr_exp_t ex, ey;
|
|
int compare, inex, q_is_odd, sign, signx = MPFR_SIGN (x);
|
|
mpz_t mx, my, r;
|
|
|
|
MPFR_ASSERTD (rnd_q == MPFR_RNDN || rnd_q == MPFR_RNDZ);
|
|
|
|
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
|
|
{
|
|
if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
|
|
|| MPFR_IS_ZERO (y))
|
|
{
|
|
/* for remquo, quo is undefined */
|
|
MPFR_SET_NAN (rem);
|
|
MPFR_RET_NAN;
|
|
}
|
|
else /* either y is Inf and x is 0 or non-special,
|
|
or x is 0 and y is non-special,
|
|
in both cases the quotient is zero. */
|
|
{
|
|
if (quo)
|
|
*quo = 0;
|
|
return mpfr_set (rem, x, rnd);
|
|
}
|
|
}
|
|
|
|
/* now neither x nor y is NaN, Inf or zero */
|
|
|
|
mpz_init (mx);
|
|
mpz_init (my);
|
|
mpz_init (r);
|
|
|
|
ex = mpfr_get_z_2exp (mx, x); /* x = mx*2^ex */
|
|
ey = mpfr_get_z_2exp (my, y); /* y = my*2^ey */
|
|
|
|
/* to get rid of sign problems, we compute it separately:
|
|
quo(-x,-y) = quo(x,y), rem(-x,-y) = -rem(x,y)
|
|
quo(-x,y) = -quo(x,y), rem(-x,y) = -rem(x,y)
|
|
thus quo = sign(x/y)*quo(|x|,|y|), rem = sign(x)*rem(|x|,|y|) */
|
|
sign = (signx == MPFR_SIGN (y)) ? 1 : -1;
|
|
mpz_abs (mx, mx);
|
|
mpz_abs (my, my);
|
|
q_is_odd = 0;
|
|
|
|
/* divide my by 2^k if possible to make operations mod my easier */
|
|
{
|
|
unsigned long k = mpz_scan1 (my, 0);
|
|
ey += k;
|
|
mpz_fdiv_q_2exp (my, my, k);
|
|
}
|
|
|
|
if (ex <= ey)
|
|
{
|
|
/* q = x/y = mx/(my*2^(ey-ex)) */
|
|
mpz_mul_2exp (my, my, ey - ex); /* divide mx by my*2^(ey-ex) */
|
|
if (rnd_q == MPFR_RNDZ)
|
|
/* 0 <= |r| <= |my|, r has the same sign as mx */
|
|
mpz_tdiv_qr (mx, r, mx, my);
|
|
else
|
|
/* 0 <= |r| <= |my|, r has the same sign as my */
|
|
mpz_fdiv_qr (mx, r, mx, my);
|
|
|
|
if (rnd_q == MPFR_RNDN)
|
|
q_is_odd = mpz_tstbit (mx, 0);
|
|
if (quo) /* mx is the quotient */
|
|
{
|
|
mpz_tdiv_r_2exp (mx, mx, WANTED_BITS);
|
|
*quo = mpz_get_si (mx);
|
|
}
|
|
}
|
|
else /* ex > ey */
|
|
{
|
|
if (quo) /* remquo case */
|
|
/* for remquo, to get the low WANTED_BITS more bits of the quotient,
|
|
we first compute R = X mod Y*2^WANTED_BITS, where X and Y are
|
|
defined below. Then the low WANTED_BITS of the quotient are
|
|
floor(R/Y). */
|
|
mpz_mul_2exp (my, my, WANTED_BITS); /* 2^WANTED_BITS*Y */
|
|
|
|
else if (rnd_q == MPFR_RNDN) /* remainder case */
|
|
/* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers.
|
|
Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y).
|
|
To be able to perform the rounding, we need the least significant
|
|
bit of the quotient, i.e., one more bit in the remainder,
|
|
which is obtained by dividing by 2Y. */
|
|
mpz_mul_2exp (my, my, 1); /* 2Y */
|
|
|
|
mpz_set_ui (r, 2);
|
|
mpz_powm_ui (r, r, ex - ey, my); /* 2^(ex-ey) mod my */
|
|
mpz_mul (r, r, mx);
|
|
mpz_mod (r, r, my);
|
|
|
|
if (quo) /* now 0 <= r < 2^WANTED_BITS*Y */
|
|
{
|
|
mpz_fdiv_q_2exp (my, my, WANTED_BITS); /* back to Y */
|
|
mpz_tdiv_qr (mx, r, r, my);
|
|
/* oldr = mx*my + newr */
|
|
*quo = mpz_get_si (mx);
|
|
q_is_odd = *quo & 1;
|
|
}
|
|
else if (rnd_q == MPFR_RNDN) /* now 0 <= r < 2Y in the remainder case */
|
|
{
|
|
mpz_fdiv_q_2exp (my, my, 1); /* back to Y */
|
|
/* least significant bit of q */
|
|
q_is_odd = mpz_cmpabs (r, my) >= 0;
|
|
if (q_is_odd)
|
|
mpz_sub (r, r, my);
|
|
}
|
|
/* now 0 <= |r| < |my|, and if needed,
|
|
q_is_odd is the least significant bit of q */
|
|
}
|
|
|
|
if (mpz_cmp_ui (r, 0) == 0)
|
|
{
|
|
inex = mpfr_set_ui (rem, 0, MPFR_RNDN);
|
|
/* take into account sign of x */
|
|
if (signx < 0)
|
|
mpfr_neg (rem, rem, MPFR_RNDN);
|
|
}
|
|
else
|
|
{
|
|
if (rnd_q == MPFR_RNDN)
|
|
{
|
|
/* FIXME: the comparison 2*r < my could be done more efficiently
|
|
at the mpn level */
|
|
mpz_mul_2exp (r, r, 1);
|
|
compare = mpz_cmpabs (r, my);
|
|
mpz_fdiv_q_2exp (r, r, 1);
|
|
compare = ((compare > 0) ||
|
|
((rnd_q == MPFR_RNDN) && (compare == 0) && q_is_odd));
|
|
/* if compare != 0, we need to subtract my to r, and add 1 to quo */
|
|
if (compare)
|
|
{
|
|
mpz_sub (r, r, my);
|
|
if (quo && (rnd_q == MPFR_RNDN))
|
|
*quo += 1;
|
|
}
|
|
}
|
|
/* take into account sign of x */
|
|
if (signx < 0)
|
|
mpz_neg (r, r);
|
|
inex = mpfr_set_z_2exp (rem, r, ex > ey ? ey : ex, rnd);
|
|
}
|
|
|
|
if (quo)
|
|
*quo *= sign;
|
|
|
|
mpz_clear (mx);
|
|
mpz_clear (my);
|
|
mpz_clear (r);
|
|
|
|
return inex;
|
|
}
|
|
|
|
int
|
|
mpfr_remainder (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
|
|
{
|
|
return mpfr_rem1 (rem, (long *) 0, MPFR_RNDN, x, y, rnd);
|
|
}
|
|
|
|
int
|
|
mpfr_remquo (mpfr_ptr rem, long *quo,
|
|
mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
|
|
{
|
|
return mpfr_rem1 (rem, quo, MPFR_RNDN, x, y, rnd);
|
|
}
|
|
|
|
int
|
|
mpfr_fmod (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
|
|
{
|
|
return mpfr_rem1 (rem, (long *) 0, MPFR_RNDZ, x, y, rnd);
|
|
}
|