1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers. 2 3 Copyright 1999-2004, 2006-2018 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 http://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) returns 33 EXP(max(|b|,|c|)) - EXP(|b| - |c|). 34 */ 35 36 int 37 mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel) 38 { 39 mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0; 40 mp_size_t bn, cn; 41 mpfr_exp_t sdiff_exp; 42 mpfr_uexp_t diff_exp; 43 mpfr_prec_t res = 0; 44 int sign; 45 46 /* b=c should not happen, since cmp2 is called only from agm (with 47 different variables) and from sub1 (if b=c, then sub1sp would be 48 called instead). So, no need for a particular optimization here. */ 49 50 /* the cases b=0 or c=0 are also treated apart in agm and sub 51 (which calls sub1) */ 52 MPFR_ASSERTD (MPFR_IS_PURE_UBF (b)); 53 MPFR_ASSERTD (MPFR_IS_PURE_UBF (c)); 54 55 sdiff_exp = MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)) ? 56 mpfr_ubf_diff_exp (b, c) : MPFR_GET_EXP (b) - MPFR_GET_EXP (c); 57 58 /* The returned result is restricted to [MPFR_EXP_MIN,MPFR_EXP_MAX], 59 which is the range of the mpfr_exp_t type. But under the condition 60 below, the value of cancel will not be affected. */ 61 MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX); 62 63 if (sdiff_exp >= 0) 64 { 65 sign = 1; 66 diff_exp = sdiff_exp; 67 68 bp = MPFR_MANT(b); 69 cp = MPFR_MANT(c); 70 71 bn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS; 72 cn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; /* # of limbs of c minus 1 */ 73 74 if (diff_exp == 0) 75 { 76 while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn]) 77 { 78 bn--; 79 cn--; 80 res += GMP_NUMB_BITS; 81 } 82 83 if (MPFR_UNLIKELY (bn < 0)) 84 { 85 if (MPFR_LIKELY (cn < 0)) /* |b| = |c| */ 86 return 0; 87 88 /* b has been read entirely, but not c. Replace b by c for the 89 symmetric case below (only the sign differs if not 0). */ 90 bp = cp; 91 bn = cn; 92 cn = -1; /* to enter the following "if" */ 93 sign = -1; 94 } 95 96 if (MPFR_UNLIKELY (cn < 0)) 97 /* c discards exactly the upper part of b */ 98 { 99 int z; 100 101 MPFR_ASSERTD (bn >= 0); 102 103 while (bp[bn] == 0) 104 { 105 if (--bn < 0) /* |b| = |c| */ 106 return 0; 107 res += GMP_NUMB_BITS; 108 } 109 110 count_leading_zeros (z, bp[bn]); /* bp[bn] <> 0 */ 111 *cancel = res + z; 112 return sign; 113 } 114 115 MPFR_ASSERTD (bn >= 0); 116 MPFR_ASSERTD (cn >= 0); 117 MPFR_ASSERTD (bp[bn] != cp[cn]); 118 if (bp[bn] < cp[cn]) 119 { 120 mp_limb_t *tp; 121 mp_size_t tn; 122 123 tp = bp; bp = cp; cp = tp; 124 tn = bn; bn = cn; cn = tn; 125 sign = -1; 126 } 127 } 128 } /* MPFR_EXP(b) >= MPFR_EXP(c) */ 129 else /* MPFR_EXP(b) < MPFR_EXP(c) */ 130 { 131 sign = -1; 132 diff_exp = - (mpfr_uexp_t) sdiff_exp; 133 134 bp = MPFR_MANT(c); 135 cp = MPFR_MANT(b); 136 137 bn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; 138 cn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS; 139 } 140 141 /* now we have removed the identical upper limbs of b and c 142 (can happen only when diff_exp = 0), and after the possible 143 swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0, 144 diff_exp = EXP(b) - EXP(c). 145 */ 146 147 if (MPFR_LIKELY (diff_exp < GMP_NUMB_BITS)) 148 { 149 cc = cp[cn] >> diff_exp; 150 /* warning: a shift by GMP_NUMB_BITS is not allowed by ISO C */ 151 if (diff_exp) 152 lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp); 153 cn--; 154 } 155 else 156 diff_exp -= GMP_NUMB_BITS; /* cc = 0 */ 157 158 dif = bp[bn--] - cc; /* necessarily dif >= 1 */ 159 MPFR_ASSERTD(dif >= 1); 160 161 /* now high_dif = 0, dif >= 1, lastc is the neglected part of cp[cn+1] */ 162 163 while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0) 164 && (high_dif == 0) && (dif == 1))) 165 { /* dif=1 implies diff_exp = 0 or 1 */ 166 MPFR_ASSERTD (diff_exp <= 1); 167 bb = (bn >= 0) ? bp[bn] : 0; 168 cc = lastc; 169 if (cn >= 0) 170 { 171 if (diff_exp == 0) 172 { 173 cc += cp[cn]; 174 } 175 else 176 { 177 MPFR_ASSERTD (diff_exp == 1); 178 cc += cp[cn] >> 1; 179 lastc = cp[cn] << (GMP_NUMB_BITS - 1); 180 } 181 } 182 else 183 lastc = 0; 184 high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1); 185 bn--; 186 cn--; 187 res += GMP_NUMB_BITS; 188 } 189 190 /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */ 191 192 if (MPFR_UNLIKELY (high_dif != 0)) /* high_dif == 1 */ 193 { 194 res--; 195 MPFR_ASSERTD (res >= 0); 196 if (dif != 0) 197 { 198 *cancel = res; 199 return sign; 200 } 201 } 202 else /* high_dif == 0 */ 203 { 204 int z; 205 206 count_leading_zeros (z, dif); /* dif > 1 here */ 207 res += z; 208 if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (GMP_NUMB_BITS - z - 1)))) 209 { /* dif is not a power of two */ 210 *cancel = res; 211 return sign; 212 } 213 } 214 215 /* now result is res + (low(b) < low(c)) */ 216 while (bn >= 0 && (cn >= 0 || lastc != 0)) 217 { 218 if (diff_exp >= GMP_NUMB_BITS) 219 diff_exp -= GMP_NUMB_BITS; 220 else 221 { 222 cc = lastc; 223 if (cn >= 0) 224 { 225 cc += cp[cn] >> diff_exp; 226 if (diff_exp != 0) 227 lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp); 228 } 229 else 230 lastc = 0; 231 cn--; 232 } 233 if (bp[bn] != cc) 234 { 235 *cancel = res + (bp[bn] < cc); 236 return sign; 237 } 238 bn--; 239 } 240 241 if (bn < 0) 242 { 243 if (lastc != 0) 244 res++; 245 else 246 { 247 while (cn >= 0 && cp[cn] == 0) 248 cn--; 249 if (cn >= 0) 250 res++; 251 } 252 } 253 254 *cancel = res; 255 return sign; 256 } 257