dma: rework config parsing
[dragonfly.git] / contrib / mpfr / rem1.c
blobc433849bd4433e1284dfd1331a5add16b2f82f6d
1 /* mpfr_rem1 -- internal function
2 mpfr_fmod -- compute the floating-point remainder of x/y
3 mpfr_remquo and mpfr_remainder -- argument reduction functions
5 Copyright 2007, 2008, 2009 Free Software Foundation, Inc.
6 Contributed by the Arenaire and Cacao projects, INRIA.
8 This file is part of the GNU MPFR Library.
10 The GNU MPFR Library is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 2.1 of the License, or (at your
13 option) any later version.
15 The GNU MPFR Library is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 License for more details.
20 You should have received a copy of the GNU Lesser General Public License
21 along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
22 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
23 MA 02110-1301, USA. */
25 # include "mpfr-impl.h"
27 /* we return as many bits as we can, keeping just one bit for the sign */
28 # define WANTED_BITS (sizeof(long) * CHAR_BIT - 1)
31 rem1 works as follows:
32 The first rounding mode rnd_q indicate if we are actually computing
33 a fmod (GMP_RNDZ) or a remainder/remquo (GMP_RNDN).
35 Let q = x/y rounded to an integer in the direction rnd_q.
36 Put x - q*y in rem, rounded according to rnd.
37 If quo is not null, the value stored in *quo has the sign of q,
38 and agrees with q with the 2^n low order bits.
39 In other words, *quo = q (mod 2^n) and *quo q >= 0.
40 If rem is zero, then it has the sign of x.
41 The returned 'int' is the inexact flag giving the place of rem wrt x - q*y.
43 If x or y is NaN: *quo is undefined, rem is NaN.
44 If x is Inf, whatever y: *quo is undefined, rem is NaN.
45 If y is Inf, x not NaN nor Inf: *quo is 0, rem is x.
46 If y is 0, whatever x: *quo is undefined, rem is NaN.
47 If x is 0, whatever y (not NaN nor 0): *quo is 0, rem is x.
49 Otherwise if x and y are neither NaN, Inf nor 0, q is always defined,
50 thus *quo is.
51 Since |x - q*y| <= y/2, no overflow is possible.
52 Only an underflow is possible when y is very small.
55 static int
56 mpfr_rem1 (mpfr_ptr rem, long *quo, mp_rnd_t rnd_q,
57 mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
59 mp_exp_t ex, ey;
60 int compare, inex, q_is_odd, sign, signx = MPFR_SIGN (x);
61 mpz_t mx, my, r;
63 MPFR_ASSERTD (rnd_q == GMP_RNDN || rnd_q == GMP_RNDZ);
65 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
67 if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
68 || MPFR_IS_ZERO (y))
70 /* for remquo, quo is undefined */
71 MPFR_SET_NAN (rem);
72 MPFR_RET_NAN;
74 else /* either y is Inf and x is 0 or non-special,
75 or x is 0 and y is non-special,
76 in both cases the quotient is zero. */
78 if (quo)
79 *quo = 0;
80 return mpfr_set (rem, x, rnd);
84 /* now neither x nor y is NaN, Inf or zero */
86 mpz_init (mx);
87 mpz_init (my);
88 mpz_init (r);
90 ex = mpfr_get_z_exp (mx, x); /* x = mx*2^ex */
91 ey = mpfr_get_z_exp (my, y); /* y = my*2^ey */
93 /* to get rid of sign problems, we compute it separately:
94 quo(-x,-y) = quo(x,y), rem(-x,-y) = -rem(x,y)
95 quo(-x,y) = -quo(x,y), rem(-x,y) = -rem(x,y)
96 thus quo = sign(x/y)*quo(|x|,|y|), rem = sign(x)*rem(|x|,|y|) */
97 sign = (signx == MPFR_SIGN (y)) ? 1 : -1;
98 mpz_abs (mx, mx);
99 mpz_abs (my, my);
100 q_is_odd = 0;
102 /* divide my by 2^k if possible to make operations mod my easier */
104 unsigned long k = mpz_scan1 (my, 0);
105 ey += k;
106 mpz_div_2exp (my, my, k);
109 if (ex <= ey)
111 /* q = x/y = mx/(my*2^(ey-ex)) */
112 mpz_mul_2exp (my, my, ey - ex); /* divide mx by my*2^(ey-ex) */
113 if (rnd_q == GMP_RNDZ)
114 /* 0 <= |r| <= |my|, r has the same sign as mx */
115 mpz_tdiv_qr (mx, r, mx, my);
116 else
117 /* 0 <= |r| <= |my|, r has the same sign as my */
118 mpz_fdiv_qr (mx, r, mx, my);
120 if (rnd_q == GMP_RNDN)
121 q_is_odd = mpz_tstbit (mx, 0);
122 if (quo) /* mx is the quotient */
124 mpz_tdiv_r_2exp (mx, mx, WANTED_BITS);
125 *quo = mpz_get_si (mx);
128 else /* ex > ey */
130 if (quo)
131 /* for remquo, to get the low WANTED_BITS more bits of the quotient,
132 we first compute R = X mod Y*2^WANTED_BITS, where X and Y are
133 defined below. Then the low WANTED_BITS of the quotient are
134 floor(R/Y). */
135 mpz_mul_2exp (my, my, WANTED_BITS); /* 2^WANTED_BITS*Y */
136 else
137 /* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers.
138 Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y).
139 To be able to perform the rounding, we need the least significant
140 bit of the quotient, i.e., one more bit in the remainder,
141 which is obtained by dividing by 2Y. */
142 mpz_mul_2exp (my, my, 1); /* 2Y */
144 mpz_set_ui (r, 2);
145 mpz_powm_ui (r, r, ex - ey, my); /* 2^(ex-ey) mod my */
146 mpz_mul (r, r, mx);
147 mpz_mod (r, r, my);
149 if (quo) /* now 0 <= r < 2^WANTED_BITS*Y */
151 mpz_div_2exp (my, my, WANTED_BITS); /* back to Y */
152 mpz_tdiv_qr (mx, r, r, my);
153 /* oldr = mx*my + newr */
154 *quo = mpz_get_si (mx);
155 q_is_odd = *quo & 1;
157 else /* now 0 <= r < 2Y */
159 mpz_div_2exp (my, my, 1); /* back to Y */
160 if (rnd_q == GMP_RNDN)
162 /* least significant bit of q */
163 q_is_odd = mpz_cmpabs (r, my) >= 0;
164 if (q_is_odd)
165 mpz_sub (r, r, my);
168 /* now 0 <= |r| < |my|, and if needed,
169 q_is_odd is the least significant bit of q */
172 if (mpz_cmp_ui (r, 0) == 0)
173 inex = mpfr_set_ui (rem, 0, GMP_RNDN);
174 else
176 if (rnd_q == GMP_RNDN)
178 /* FIXME: the comparison 2*r < my could be done more efficiently
179 at the mpn level */
180 mpz_mul_2exp (r, r, 1);
181 compare = mpz_cmpabs (r, my);
182 mpz_div_2exp (r, r, 1);
183 compare = ((compare > 0) ||
184 ((rnd_q == GMP_RNDN) && (compare == 0) && q_is_odd));
185 /* if compare != 0, we need to subtract my to r, and add 1 to quo */
186 if (compare)
188 mpz_sub (r, r, my);
189 if (quo && (rnd_q == GMP_RNDN))
190 *quo += 1;
193 inex = mpfr_set_z (rem, r, rnd);
194 /* if ex > ey, rem should be multiplied by 2^ey, else by 2^ex */
195 MPFR_EXP (rem) += (ex > ey) ? ey : ex;
198 if (quo)
199 *quo *= sign;
201 /* take into account sign of x */
202 if (signx < 0)
204 mpfr_neg (rem, rem, GMP_RNDN);
205 inex = -inex;
208 mpz_clear (mx);
209 mpz_clear (my);
210 mpz_clear (r);
212 return inex;
216 mpfr_remainder (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
218 return mpfr_rem1 (rem, (long *) 0, GMP_RNDN, x, y, rnd);
222 mpfr_remquo (mpfr_ptr rem, long *quo,
223 mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
225 return mpfr_rem1 (rem, quo, GMP_RNDN, x, y, rnd);
229 mpfr_fmod (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd)
231 return mpfr_rem1 (rem, (long *) 0, GMP_RNDZ, x, y, rnd);