1/* mpfr_cmp2 -- exponent shift when subtracting two numbers. 2 3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 36int 37mpfr_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