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-2023 Free Software Foundation, Inc.
6 Contributed by the AriC and Caramba 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 https://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 unspecified, rem is NaN.
44 If x is Inf, whatever y: *quo is unspecified, 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 unspecified, 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 int tiny = 0;
63
64 MPFR_ASSERTD (rnd_q == MPFR_RNDN || rnd_q == MPFR_RNDZ);
65
66 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y)))
67 {
68 if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x)
69 || MPFR_IS_ZERO (y))
70 {
71 /* for remquo, *quo is unspecified */
72 MPFR_SET_NAN (rem);
73 MPFR_RET_NAN;
74 }
75 else /* either y is Inf and x is 0 or non-special,
76 or x is 0 and y is non-special,
77 in both cases the quotient is zero. */
78 {
79 if (quo)
80 *quo = 0;
81 return mpfr_set (rem, x, rnd);
82 }
83 }
84
85 /* now neither x nor y is NaN, Inf or zero */
86
87 mpz_init (mx);
88 mpz_init (my);
89 mpz_init (r);
90
91 ex = mpfr_get_z_2exp (mx, x); /* x = mx*2^ex */
92 ey = mpfr_get_z_2exp (my, y); /* y = my*2^ey */
93
94 /* to get rid of sign problems, we compute it separately:
95 quo(-x,-y) = quo(x,y), rem(-x,-y) = -rem(x,y)
96 quo(-x,y) = -quo(x,y), rem(-x,y) = -rem(x,y)
97 thus quo = sign(x/y)*quo(|x|,|y|), rem = sign(x)*rem(|x|,|y|) */
98 sign = (signx == MPFR_SIGN (y)) ? 1 : -1;
99 mpz_abs (mx, mx);
100 mpz_abs (my, my);
101 q_is_odd = 0;
102
103 /* Divide my by 2^k if possible to make operations mod my easier.
104 Since my comes from a regular MPFR number, due to the constraints on the
105 exponent and the precision, there can be no integer overflow below. */
106 {
107 mpfr_exp_t k = mpz_scan1 (my, 0);
108 ey += k;
109 mpz_fdiv_q_2exp (my, my, k);
110 }
111
112 if (ex <= ey)
113 {
114 /* q = x/y = mx/(my*2^(ey-ex)) */
115
116 /* First detect cases where q=0, to avoid creating a huge number
117 my*2^(ey-ex): if sx = mpz_sizeinbase (mx, 2) and sy =
118 mpz_sizeinbase (my, 2), we have x < 2^(ex + sx) and
119 y >= 2^(ey + sy - 1), thus if ex + sx <= ey + sy - 1
120 the quotient is 0 */
121 if (ex + (mpfr_exp_t) mpz_sizeinbase (mx, 2) <
122 ey + (mpfr_exp_t) mpz_sizeinbase (my, 2))
123 {
124 tiny = 1;
125 mpz_set (r, mx);
126 mpz_set_ui (mx, 0);
127 }
128 else
129 {
130 mpz_mul_2exp (my, my, ey - ex); /* divide mx by my*2^(ey-ex) */
131
132 /* since mx > 0 and my > 0, we can use mpz_tdiv_qr in all cases */
133 mpz_tdiv_qr (mx, r, mx, my);
134 /* 0 <= |r| <= |my|, r has the same sign as mx */
135 }
136
137 if (rnd_q == MPFR_RNDN)
138 q_is_odd = mpz_tstbit (mx, 0);
139 if (quo) /* mx is the quotient */
140 {
141 mpz_tdiv_r_2exp (mx, mx, WANTED_BITS);
142 *quo = mpz_get_si (mx);
143 }
144 }
145 else /* ex > ey */
146 {
147 if (quo) /* remquo case */
148 /* for remquo, to get the low WANTED_BITS more bits of the quotient,
149 we first compute R = X mod Y*2^WANTED_BITS, where X and Y are
150 defined below. Then the low WANTED_BITS of the quotient are
151 floor(R/Y). */
152 mpz_mul_2exp (my, my, WANTED_BITS); /* 2^WANTED_BITS*Y */
153
154 else if (rnd_q == MPFR_RNDN) /* remainder case */
155 /* Let X = mx*2^(ex-ey) and Y = my. Then both X and Y are integers.
156 Assume X = R mod Y, then x = X*2^ey = R*2^ey mod (Y*2^ey=y).
157 To be able to perform the rounding, we need the least significant
158 bit of the quotient, i.e., one more bit in the remainder,
159 which is obtained by dividing by 2Y. */
160 mpz_mul_2exp (my, my, 1); /* 2Y */
161
162 /* Warning: up to GMP 6.2.0, mpz_powm_ui is not optimized when BASE^EXP
163 has about the same size as MOD, in which case it should first compute
164 BASE^EXP exactly, then reduce it modulo MOD:
165 https://gmplib.org/list-archives/gmp-bugs/2020-February/004736.html
166 Thus when 2^(ex-ey) is less than my^3, we use this algorithm. */
167 if (ex - ey > 3 * mpz_sizeinbase (my, 2))
168 {
169 mpz_set_ui (r, 2);
170 mpz_powm_ui (r, r, ex - ey, my); /* 2^(ex-ey) mod my */
171 }
172 else
173 mpz_ui_pow_ui (r, 2, ex - ey);
174 mpz_mul (r, r, mx);
175 mpz_mod (r, r, my);
176
177 if (quo) /* now 0 <= r < 2^WANTED_BITS*Y */
178 {
179 mpz_fdiv_q_2exp (my, my, WANTED_BITS); /* back to Y */
180 mpz_tdiv_qr (mx, r, r, my);
181 /* oldr = mx*my + newr */
182 *quo = mpz_get_si (mx);
183 q_is_odd = *quo & 1;
184 }
185 else if (rnd_q == MPFR_RNDN) /* now 0 <= r < 2Y in the remainder case */
186 {
187 mpz_fdiv_q_2exp (my, my, 1); /* back to Y */
188 /* least significant bit of q */
189 q_is_odd = mpz_cmpabs (r, my) >= 0;
190 if (q_is_odd)
191 mpz_sub (r, r, my);
192 }
193 /* now 0 <= |r| < |my|, and if needed,
194 q_is_odd is the least significant bit of q */
195 }
196
197 if (mpz_cmp_ui (r, 0) == 0)
198 {
199 inex = mpfr_set_ui (rem, 0, MPFR_RNDN);
200 /* take into account sign of x */
201 if (signx < 0)
202 mpfr_neg (rem, rem, MPFR_RNDN);
203 }
204 else
205 {
206 if (rnd_q == MPFR_RNDN)
207 {
208 /* FIXME: the comparison 2*r < my could be done more efficiently
209 at the mpn level */
210 mpz_mul_2exp (r, r, 1);
211 /* if tiny=1, we should compare r with my*2^(ey-ex) */
212 if (tiny)
213 {
214 if (ex + (mpfr_exp_t) mpz_sizeinbase (r, 2) <
215 ey + (mpfr_exp_t) mpz_sizeinbase (my, 2))
216 compare = 0; /* r*2^ex < my*2^ey */
217 else
218 {
219 mpz_mul_2exp (my, my, ey - ex);
220 compare = mpz_cmpabs (r, my);
221 }
222 }
223 else
224 compare = mpz_cmpabs (r, my);
225 mpz_fdiv_q_2exp (r, r, 1);
226 compare = ((compare > 0) ||
227 ((rnd_q == MPFR_RNDN) && (compare == 0) && q_is_odd));
228 /* if compare != 0, we need to subtract my to r, and add 1 to quo */
229 if (compare)
230 {
231 mpz_sub (r, r, my);
232 if (quo && (rnd_q == MPFR_RNDN))
233 *quo += 1;
234 }
235 }
236 /* take into account sign of x */
237 if (signx < 0)
238 mpz_neg (r, r);
239 inex = mpfr_set_z_2exp (rem, r, ex > ey ? ey : ex, rnd);
240 }
241
242 if (quo)
243 *quo *= sign;
244
245 mpz_clear (mx);
246 mpz_clear (my);
247 mpz_clear (r);
248
249 return inex;
250 }
251
252 int
253 mpfr_remainder (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
254 {
255 return mpfr_rem1 (rem, (long *) 0, MPFR_RNDN, x, y, rnd);
256 }
257
258 int
259 mpfr_remquo (mpfr_ptr rem, long *quo,
260 mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
261 {
262 return mpfr_rem1 (rem, quo, MPFR_RNDN, x, y, rnd);
263 }
264
265 int
266 mpfr_fmod (mpfr_ptr rem, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd)
267 {
268 return mpfr_rem1 (rem, (long *) 0, MPFR_RNDZ, x, y, rnd);
269 }
270
271 int
272 mpfr_fmodquo (mpfr_ptr rem, long *quo, mpfr_srcptr x, mpfr_srcptr y,
273 mpfr_rnd_t rnd)
274 {
275 return mpfr_rem1 (rem, quo, MPFR_RNDZ, x, y, rnd);
276 }