1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers. 2 3 Copyright 1999-2004, 2006-2023 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 24 #define MPFR_NEED_LONGLONG_H 25 #include "mpfr-impl.h" 26 27 /* If |b| != |c|, puts the number of canceled bits when one subtracts |c| 28 from |b| in *cancel. Returns the sign of the difference (-1, 0, 1). 29 30 Assumes neither of b or c is NaN, +/- infinity, or +/- 0. 31 32 In other terms, if |b| != |c|, mpfr_cmp2 (b, c) stores 33 EXP(max(|b|,|c|)) - EXP(|b| - |c|) in *cancel. 34 35 One necessarily has 0 <= cancel <= max(PREC(b),PREC(c)), so that this 36 value is representable in a mpfr_prec_t. Note that in the code, the 37 maximum intermediate value is cancel + 1, but since MPFR_PREC_MAX is 38 not the maximum value of mpfr_prec_t, there is no integer overflow. 39 */ 40 41 int 42 mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) 43 { 44 mp_limb_t *bp, *cp, bb, cc, lastc, dif; 45 int high_dif; /* manipulated like a boolean */ 46 mp_size_t bn, cn; 47 mpfr_exp_t sdiff_exp; 48 mpfr_uexp_t diff_exp; 49 mpfr_prec_t res = 0; /* will be the number of canceled bits (*cancel) */ 50 int sign; 51 52 /* b=c should not happen, since cmp2 is called only from agm (with 53 different variables) and from sub1 (if b=c, then sub1sp would be 54 called instead). So, no need for a particular optimization here. */ 55 56 /* the cases b=0 or c=0 are also treated apart in agm and sub 57 (which calls sub1) */ 58 MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); 59 MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); 60 61 sdiff_exp = MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)) ? 62 mpfr_ubf_diff_exp (b, c) : MPFR_GET_EXP (b) - MPFR_GET_EXP (c); 63 64 /* The returned result is saturated to [MPFR_EXP_MIN,MPFR_EXP_MAX], 65 which is the range of the mpfr_exp_t type. But under the condition 66 below, since |MPFR_EXP_MIN| >= MPFR_EXP_MAX, the value of cancel 67 will not be affected: by symmetry (as done in the code), assume 68 |b| >= |c|; if EXP(b) - EXP(c) >= MPFR_EXP_MAX, then |c| < ulp(b), 69 so that the value of cancel is 0, or 1 if |b| is a power of 2, 70 whatever the exact value of EXP(b) - EXP(c). */ 71 MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX); 72 73 if (sdiff_exp >= 0) 74 { 75 sign = 1; /* assumes |b| > |c|; will be changed if not. */ 76 diff_exp = sdiff_exp; 77 78 bp = MPFR_MANT(b); 79 cp = MPFR_MANT(c); 80 81 /* index of the most significant limb of b and c */ 82 bn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS; 83 cn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; 84 85 /* If diff_exp != 0, i.e. diff_exp > 0, then |b| > |c|. Otherwise... */ 86 if (diff_exp == 0) 87 { 88 /* Skip the identical most significant limbs, adding GMP_NUMB_BITS 89 to the number of canceled bits at each iteration. */ 90 while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn]) 91 { 92 bn--; 93 cn--; 94 res += GMP_NUMB_BITS; 95 } 96 97 if (MPFR_UNLIKELY (bn < 0)) 98 { 99 if (MPFR_LIKELY (cn < 0)) /* |b| = |c| */ 100 return 0; 101 102 /* b has been read entirely, but not c. Thus |b| <= |c|. 103 Swap (bp,bn) and (cp,cn), and take the opposite sign 104 for the symmetric case below (simulating a swap). 105 Note: cp will not be used, thus is not assigned; and 106 "cn = -1;" is necessary to enter the following "if" 107 (probably less confusing than a "goto"). */ 108 bp = cp; 109 bn = cn; 110 cn = -1; 111 sign = -1; 112 } 113 114 if (MPFR_UNLIKELY (cn < 0)) 115 /* c discards exactly the upper part of b */ 116 { 117 int z; 118 119 MPFR_ASSERTD (bn >= 0); 120 121 /* Skip null limbs of b (= non-represented null limbs of c), 122 adding GMP_NUMB_BITS to the number of canceled bits at 123 each iteration. */ 124 while (bp[bn] == 0) 125 { 126 if (--bn < 0) /* |b| = |c| */ 127 return 0; 128 res += GMP_NUMB_BITS; 129 } 130 131 count_leading_zeros (z, bp[bn]); /* bp[bn] != 0 */ 132 *cancel = res + z; 133 return sign; 134 } 135 136 MPFR_ASSERTD (bn >= 0); 137 MPFR_ASSERTD (cn >= 0); 138 MPFR_ASSERTD (bp[bn] != cp[cn]); 139 140 /* |b| != |c|. If |b| < |c|: swap (bp,bn) and (cp,cn), 141 and take the opposite sign. */ 142 if (bp[bn] < cp[cn]) 143 { 144 mp_limb_t *tp; 145 mp_size_t tn; 146 147 tp = bp; bp = cp; cp = tp; 148 tn = bn; bn = cn; cn = tn; 149 sign = -1; 150 } 151 } 152 } /* MPFR_EXP(b) >= MPFR_EXP(c) */ 153 else /* MPFR_EXP(b) < MPFR_EXP(c) */ 154 { 155 /* We necessarily have |b| < |c|. Simulate a swap by reading the 156 parameters so that |(bp,bn)| > |(cp,cn)|. */ 157 158 sign = -1; 159 diff_exp = - (mpfr_uexp_t) sdiff_exp; 160 161 bp = MPFR_MANT(c); 162 cp = MPFR_MANT(b); 163 164 bn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; 165 cn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS; 166 } 167 168 /* Now we have removed the identical upper limbs of b and c 169 (when diff_exp = 0), and after the possible swap, we have |b| > |c|, 170 where b is represented by (bp,bn) and c is represented by (cp,cn). 171 The value diff_exp = EXP(b) - EXP(c) can be regarded as the number 172 of leading zeros of c, when aligned with b. */ 173 174 /* When a limb of c is read from memory, the part that is not taken 175 into account for the operation with a limb bp[bn] of b will be put 176 in lastc, shifted to the leftmost part (for alignment with b): 177 [-------- bp[bn] --------][------- bp[bn-1] -------] 178 [-- old_lastc --][-------- cp[cn] --------] 179 [-- new_lastc --] 180 Note: if diff_exp == 0, then lastc will always remain 0. */ 181 lastc = 0; 182 183 /* Compute the next limb difference, which cannot be 0 (dif >= 1). */ 184 185 if (MPFR_LIKELY (diff_exp < GMP_NUMB_BITS)) 186 { 187 cc = cp[cn] >> diff_exp; 188 /* warning: a shift by GMP_NUMB_BITS is not allowed by ISO C */ 189 if (diff_exp != 0) 190 lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp); 191 cn--; 192 } 193 else 194 { 195 cc = 0; 196 diff_exp -= GMP_NUMB_BITS; /* remove GMP_NUMB_BITS leading zeros */ 197 } 198 199 MPFR_ASSERTD (bp[bn] >= cc); /* no borrow out in subtraction below */ 200 dif = bp[bn--] - cc; 201 MPFR_ASSERTD (dif >= 1); 202 high_dif = 0; 203 204 /* The current difference, here and later, is expressed under the form 205 [high_dif][dif], where high_dif is 0 or 1, and dif is a limb. 206 Here, since we have computed a difference of limbs (with b >= c), 207 high_dif = 0. */ 208 209 /* One needs to accumulate canceled bits for the remaining case where 210 b and c are close to each other due to a long borrow propagation: 211 b = [common part]1000...000[low(b)] 212 c = [common part]0111...111[low(c)] 213 After eliminating the common part above, we have computed a difference 214 of the most significant parts, which has been stored in [high_dif][dif] 215 with high_dif = 0. We will loop as long as the currently computed 216 difference [high_dif][dif] = 1 (it is >= 1 by construction). The 217 computation of the difference will be: 218 1bbb...bbb 219 - ccc...ccc 220 where the leading 1 before bbb...bbb corresponds to [high_dif][dif] 221 at the beginning of the loop. We will exit the loop also when c has 222 entirely been taken into account as cancellation is no longer possible 223 in this case (it is no longer possible to cancel the leading 1). 224 Note: We can enter the loop only with diff_exp = 0 (with a non-empty 225 common part, partly or entirely removed) or with diff_exp = 1 (with 226 an empty common part). Indeed, if diff_exp > 1, then no limbs have 227 been skipped, so that bp[bn] had its MSB equal to 1 and the most two 228 significant bits of cc are 0, which implies that dif > 1. */ 229 230 while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0) 231 && high_dif == 0 && dif == 1)) 232 { 233 /* Since we consider the next limb, we assume a cancellation of 234 GMP_NUMB_BITS (the new exponent of the difference now being the 235 one of the MSB of the next limb). But if the leading 1 remains 236 1 in the difference (i.e. high_dif = 1 at the end of the loop), 237 then we will need to decrease res. */ 238 res += GMP_NUMB_BITS; 239 MPFR_ASSERTD (diff_exp <= 1); /* see comment before the loop */ 240 bb = bn >= 0 ? bp[bn--] : 0; /* next limb of b or non-represented 0 */ 241 if (MPFR_UNLIKELY (cn < 0)) 242 { 243 cc = lastc; 244 lastc = 0; 245 } 246 else if (diff_exp == 0) 247 { 248 cc = cp[cn--]; 249 } 250 else 251 { 252 MPFR_ASSERTD (diff_exp == 1); 253 MPFR_ASSERTD (lastc == 0 || lastc == MPFR_LIMB_HIGHBIT); 254 cc = lastc + (cp[cn] >> 1); 255 lastc = cp[cn--] << (GMP_NUMB_BITS - 1); 256 } 257 dif = bb - cc; 258 high_dif = bb >= cc; 259 } 260 261 /* Now, c has entirely been taken into account or [high_dif][dif] > 1. 262 In any case, [high_dif][dif] >= 1 by construction. 263 First, we determine the currently number of canceled bits, 264 corresponding to the exponent of the current difference. 265 The trailing bits of c, if any, can still decrease the exponent of 266 the difference when [high_dif][dif] is a power of two, but since 267 [high_dif][dif] > 1 in this case, by not more than 1. */ 268 269 if (high_dif != 0) /* high_dif == 1 */ 270 { 271 res--; /* see comment at the beginning of the above loop */ 272 MPFR_ASSERTD (res >= 0); 273 /* Terminate if [high_dif][dif] is not a power of two. */ 274 if (MPFR_LIKELY (dif != 0)) 275 goto end; 276 } 277 else /* high_dif == 0 */ 278 { 279 int z; 280 281 MPFR_ASSERTD (dif >= 1); /* [high_dif][dif] >= 1 */ 282 count_leading_zeros (z, dif); 283 res += z; 284 /* Terminate if [high_dif][dif] is not a power of two. */ 285 if (MPFR_LIKELY (NOT_POW2 (dif))) 286 goto end; 287 } 288 289 /* Now, the result will be res + (low(b) < low(c)). */ 290 291 /* If c has entirely been taken into account, it can no longer modify 292 the current result. */ 293 if (cn < 0 && lastc == 0) 294 goto end; 295 296 for (; bn >= 0 ; bn--) 297 { 298 if (diff_exp >= GMP_NUMB_BITS) 299 { 300 diff_exp -= GMP_NUMB_BITS; 301 MPFR_ASSERTD (cc == 0); 302 } 303 else if (MPFR_UNLIKELY (cn < 0)) 304 { 305 cc = lastc; 306 lastc = 0; 307 } 308 else if (diff_exp == 0) 309 { 310 cc = cp[cn--]; 311 } 312 else 313 { 314 MPFR_ASSERTD (diff_exp >= 1 && diff_exp < GMP_NUMB_BITS); 315 cc = lastc + (cp[cn] >> diff_exp); 316 lastc = cp[cn--] << (GMP_NUMB_BITS - diff_exp); 317 } 318 319 if (bp[bn] != cc) 320 { 321 res += bp[bn] < cc; 322 goto end; 323 } 324 } 325 326 /* b has entirely been read. Determine whether the trailing part of c 327 is non-zero. */ 328 329 if (lastc != 0) 330 res++; 331 else 332 { 333 while (cn >= 0 && cp[cn] == 0) 334 cn--; 335 if (cn >= 0) 336 res++; 337 } 338 339 end: 340 *cancel = res; 341 return sign; 342 } 343