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 } 277