xref: /dragonfly/contrib/mpfr/src/rem1.c (revision 25a2db75)
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
4 
5 Copyright 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
6 Contributed by the AriC and Caramel projects, INRIA.
7 
8 This file is part of the GNU MPFR Library.
9 
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 3 of the License, or (at your
13 option) any later version.
14 
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.
19 
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.LESSER.  If not, see
22 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
23 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 
25 # include "mpfr-impl.h"
26 
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)
29 
30 /*
31   rem1 works as follows:
32   The first rounding mode rnd_q indicate if we are actually computing
33   a fmod (MPFR_RNDZ) or a remainder/remquo (MPFR_RNDN).
34 
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.
42 
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.
48 
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.
53  */
54 
55 static int
56 mpfr_rem1 (mpfr_ptr rem, long *quo, mpfr_rnd_t rnd_q,
57            mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
58 {
59   mpfr_exp_t ex, ey;
60   int compare, inex, q_is_odd, sign, signx = MPFR_SIGN (x);
61   mpz_t mx, my, r;
62 
63   MPFR_ASSERTD (rnd_q == MPFR_RNDN || rnd_q == MPFR_RNDZ);
64 
65   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
66     {
67       if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
68           || MPFR_IS_ZERO (y))
69         {
70           /* for remquo, quo is undefined */
71           MPFR_SET_NAN (rem);
72           MPFR_RET_NAN;
73         }
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. */
77         {
78           if (quo)
79             *quo = 0;
80           return mpfr_set (rem, x, rnd);
81         }
82     }
83 
84   /* now neither x nor y is NaN, Inf or zero */
85 
86   mpz_init (mx);
87   mpz_init (my);
88   mpz_init (r);
89 
90   ex = mpfr_get_z_2exp (mx, x);  /* x = mx*2^ex */
91   ey = mpfr_get_z_2exp (my, y);  /* y = my*2^ey */
92 
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;
101 
102   /* divide my by 2^k if possible to make operations mod my easier */
103   {
104     unsigned long k = mpz_scan1 (my, 0);
105     ey += k;
106     mpz_fdiv_q_2exp (my, my, k);
107   }
108 
109   if (ex <= ey)
110     {
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 == MPFR_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);
119 
120       if (rnd_q == MPFR_RNDN)
121         q_is_odd = mpz_tstbit (mx, 0);
122       if (quo)                  /* mx is the quotient */
123         {
124           mpz_tdiv_r_2exp (mx, mx, WANTED_BITS);
125           *quo = mpz_get_si (mx);
126         }
127     }
128   else                          /* ex > ey */
129     {
130       if (quo) /* remquo case */
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 
137       else if (rnd_q == MPFR_RNDN) /* remainder case */
138         /* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers.
139            Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y).
140            To be able to perform the rounding, we need the least significant
141            bit of the quotient, i.e., one more bit in the remainder,
142            which is obtained by dividing by 2Y. */
143         mpz_mul_2exp (my, my, 1);       /* 2Y */
144 
145       mpz_set_ui (r, 2);
146       mpz_powm_ui (r, r, ex - ey, my);  /* 2^(ex-ey) mod my */
147       mpz_mul (r, r, mx);
148       mpz_mod (r, r, my);
149 
150       if (quo)                  /* now 0 <= r < 2^WANTED_BITS*Y */
151         {
152           mpz_fdiv_q_2exp (my, my, WANTED_BITS);   /* back to Y */
153           mpz_tdiv_qr (mx, r, r, my);
154           /* oldr = mx*my + newr */
155           *quo = mpz_get_si (mx);
156           q_is_odd = *quo & 1;
157         }
158       else if (rnd_q == MPFR_RNDN) /* now 0 <= r < 2Y in the remainder case */
159         {
160           mpz_fdiv_q_2exp (my, my, 1);     /* back to Y */
161           /* least significant bit of q */
162           q_is_odd = mpz_cmpabs (r, my) >= 0;
163           if (q_is_odd)
164             mpz_sub (r, r, my);
165         }
166       /* now 0 <= |r| < |my|, and if needed,
167          q_is_odd is the least significant bit of q */
168     }
169 
170   if (mpz_cmp_ui (r, 0) == 0)
171     {
172       inex = mpfr_set_ui (rem, 0, MPFR_RNDN);
173       /* take into account sign of x */
174       if (signx < 0)
175         mpfr_neg (rem, rem, MPFR_RNDN);
176     }
177   else
178     {
179       if (rnd_q == MPFR_RNDN)
180         {
181           /* FIXME: the comparison 2*r < my could be done more efficiently
182              at the mpn level */
183           mpz_mul_2exp (r, r, 1);
184           compare = mpz_cmpabs (r, my);
185           mpz_fdiv_q_2exp (r, r, 1);
186           compare = ((compare > 0) ||
187                      ((rnd_q == MPFR_RNDN) && (compare == 0) && q_is_odd));
188           /* if compare != 0, we need to subtract my to r, and add 1 to quo */
189           if (compare)
190             {
191               mpz_sub (r, r, my);
192               if (quo && (rnd_q == MPFR_RNDN))
193                 *quo += 1;
194             }
195         }
196       /* take into account sign of x */
197       if (signx < 0)
198         mpz_neg (r, r);
199       inex = mpfr_set_z_2exp (rem, r, ex > ey ? ey : ex, rnd);
200     }
201 
202   if (quo)
203     *quo *= sign;
204 
205   mpz_clear (mx);
206   mpz_clear (my);
207   mpz_clear (r);
208 
209   return inex;
210 }
211 
212 int
213 mpfr_remainder (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
214 {
215   return mpfr_rem1 (rem, (long *) 0, MPFR_RNDN, x, y, rnd);
216 }
217 
218 int
219 mpfr_remquo (mpfr_ptr rem, long *quo,
220              mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
221 {
222   return mpfr_rem1 (rem, quo, MPFR_RNDN, x, y, rnd);
223 }
224 
225 int
226 mpfr_fmod (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
227 {
228   return mpfr_rem1 (rem, (long *) 0, MPFR_RNDZ, x, y, rnd);
229 }
230