1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers. 2 3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4 Contributed by the Arenaire and Caramel 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_uexp_t diff_exp; 42 mpfr_prec_t res = 0; 43 int sign; 44 45 /* b=c should not happen, since cmp2 is called only from agm (with 46 different variables) and from sub1 (if b=c, then sub1sp would be 47 called instead). So, no need for a particular optimization here. */ 48 49 /* the cases b=0 or c=0 are also treated apart in agm and sub 50 (which calls sub1) */ 51 MPFR_ASSERTD (MPFR_IS_PURE_FP(b)); 52 MPFR_ASSERTD (MPFR_IS_PURE_FP(c)); 53 54 if (MPFR_GET_EXP (b) >= MPFR_GET_EXP (c)) 55 { 56 sign = 1; 57 diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c); 58 59 bp = MPFR_MANT(b); 60 cp = MPFR_MANT(c); 61 62 bn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS; 63 cn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; /* # of limbs of c minus 1 */ 64 65 if (MPFR_UNLIKELY( diff_exp == 0 )) 66 { 67 while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn]) 68 { 69 bn--; 70 cn--; 71 res += GMP_NUMB_BITS; 72 } 73 74 if (MPFR_UNLIKELY (bn < 0)) 75 { 76 if (MPFR_LIKELY (cn < 0)) /* b = c */ 77 return 0; 78 79 bp = cp; 80 bn = cn; 81 cn = -1; 82 sign = -1; 83 } 84 85 if (MPFR_UNLIKELY (cn < 0)) 86 /* c discards exactly the upper part of b */ 87 { 88 unsigned int z; 89 90 MPFR_ASSERTD (bn >= 0); 91 92 while (bp[bn] == 0) 93 { 94 if (--bn < 0) /* b = c */ 95 return 0; 96 res += GMP_NUMB_BITS; 97 } 98 99 count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */ 100 *cancel = res + z; 101 return sign; 102 } 103 104 MPFR_ASSERTD (bn >= 0); 105 MPFR_ASSERTD (cn >= 0); 106 MPFR_ASSERTD (bp[bn] != cp[cn]); 107 if (bp[bn] < cp[cn]) 108 { 109 mp_limb_t *tp; 110 mp_size_t tn; 111 112 tp = bp; bp = cp; cp = tp; 113 tn = bn; bn = cn; cn = tn; 114 sign = -1; 115 } 116 } 117 } /* MPFR_EXP(b) >= MPFR_EXP(c) */ 118 else /* MPFR_EXP(b) < MPFR_EXP(c) */ 119 { 120 sign = -1; 121 diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (c) - MPFR_GET_EXP (b); 122 123 bp = MPFR_MANT(c); 124 cp = MPFR_MANT(b); 125 126 bn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; 127 cn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS; 128 } 129 130 /* now we have removed the identical upper limbs of b and c 131 (can happen only when diff_exp = 0), and after the possible 132 swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0, 133 diff_exp = EXP(b) - EXP(c). 134 */ 135 136 if (MPFR_LIKELY (diff_exp < GMP_NUMB_BITS)) 137 { 138 cc = cp[cn] >> diff_exp; 139 /* warning: a shift by GMP_NUMB_BITS may give wrong results */ 140 if (diff_exp) 141 lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp); 142 cn--; 143 } 144 else 145 diff_exp -= GMP_NUMB_BITS; /* cc = 0 */ 146 147 dif = bp[bn--] - cc; /* necessarily dif >= 1 */ 148 MPFR_ASSERTD(dif >= 1); 149 150 /* now high_dif = 0, dif >= 1, lastc is the neglected part of cp[cn+1] */ 151 152 while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0) 153 && (high_dif == 0) && (dif == 1))) 154 { /* dif=1 implies diff_exp = 0 or 1 */ 155 bb = (bn >= 0) ? bp[bn] : 0; 156 cc = lastc; 157 if (cn >= 0) 158 { 159 if (diff_exp == 0) 160 { 161 cc += cp[cn]; 162 } 163 else /* diff_exp = 1 */ 164 { 165 cc += cp[cn] >> 1; 166 lastc = cp[cn] << (GMP_NUMB_BITS - 1); 167 } 168 } 169 else 170 lastc = 0; 171 high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1); 172 bn--; 173 cn--; 174 res += GMP_NUMB_BITS; 175 } 176 177 /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */ 178 179 if (MPFR_UNLIKELY (high_dif != 0)) /* high_dif == 1 */ 180 { 181 res--; 182 MPFR_ASSERTD (res >= 0); 183 if (dif != 0) 184 { 185 *cancel = res; 186 return sign; 187 } 188 } 189 else /* high_dif == 0 */ 190 { 191 unsigned int z; 192 193 count_leading_zeros(z, dif); /* dif > 1 here */ 194 res += z; 195 if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (GMP_NUMB_BITS - z - 1)))) 196 { /* dif is not a power of two */ 197 *cancel = res; 198 return sign; 199 } 200 } 201 202 /* now result is res + (low(b) < low(c)) */ 203 while (MPFR_UNLIKELY (bn >= 0 && (cn >= 0 || lastc != 0))) 204 { 205 if (diff_exp >= GMP_NUMB_BITS) 206 diff_exp -= GMP_NUMB_BITS; 207 else 208 { 209 cc = lastc; 210 if (cn >= 0) 211 { 212 cc += cp[cn] >> diff_exp; 213 if (diff_exp != 0) 214 lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp); 215 } 216 else 217 lastc = 0; 218 cn--; 219 } 220 if (bp[bn] != cc) 221 { 222 *cancel = res + (bp[bn] < cc); 223 return sign; 224 } 225 bn--; 226 } 227 228 if (bn < 0) 229 { 230 if (lastc != 0) 231 res++; 232 else 233 { 234 while (cn >= 0 && cp[cn] == 0) 235 cn--; 236 if (cn >= 0) 237 res++; 238 } 239 } 240 241 *cancel = res; 242 return sign; 243 } 244