1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers.
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library.
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 2.1 of the License, or (at your
11 option) any later version.
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.
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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
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.
30 Assumes neither of b or c is NaN, +/- infinity, or +/- 0.
32 In other terms, if |b| != |c|, mpfr_cmp2 (b, c) returns
33 EXP(max(|b|,|c|)) - EXP(|b| - |c|).
37 mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mp_prec_t *cancel)
39 mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0;
41 mp_exp_unsigned_t diff_exp;
45 /* b=c should not happen, since cmp2 is called only from agm
46 (with different variables), and from sub1 (if same b=c, then
47 sub1sp would be called instead */
48 MPFR_ASSERTD (b != c);
50 /* the cases b=0 or c=0 are also treated apart in agm and sub
52 MPFR_ASSERTD (MPFR_IS_PURE_FP(b));
53 MPFR_ASSERTD (MPFR_IS_PURE_FP(c));
55 if (MPFR_GET_EXP (b) >= MPFR_GET_EXP (c))
58 diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
63 bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
64 cn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB; /* # of limbs of c minus 1 */
66 if (MPFR_UNLIKELY( diff_exp == 0 ))
68 while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn])
72 res += BITS_PER_MP_LIMB;
75 if (MPFR_UNLIKELY (bn < 0))
77 if (MPFR_LIKELY (cn < 0)) /* b = c */
86 if (MPFR_UNLIKELY (cn < 0))
87 /* c discards exactly the upper part of b */
91 MPFR_ASSERTD (bn >= 0);
95 if (--bn < 0) /* b = c */
97 res += BITS_PER_MP_LIMB;
100 count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
105 MPFR_ASSERTD (bn >= 0);
106 MPFR_ASSERTD (cn >= 0);
107 MPFR_ASSERTD (bp[bn] != cp[cn]);
113 tp = bp; bp = cp; cp = tp;
114 tn = bn; bn = cn; cn = tn;
118 } /* MPFR_EXP(b) >= MPFR_EXP(c) */
119 else /* MPFR_EXP(b) < MPFR_EXP(c) */
122 diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (c) - MPFR_GET_EXP (b);
127 bn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
128 cn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
131 /* now we have removed the identical upper limbs of b and c
132 (can happen only when diff_exp = 0), and after the possible
133 swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0,
134 diff_exp = EXP(b) - EXP(c).
137 if (MPFR_LIKELY (diff_exp < BITS_PER_MP_LIMB))
139 cc = cp[cn] >> diff_exp;
140 /* warning: a shift by BITS_PER_MP_LIMB may give wrong results */
142 lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
146 diff_exp -= BITS_PER_MP_LIMB; /* cc = 0 */
148 dif = bp[bn--] - cc; /* necessarily dif >= 1 */
149 MPFR_ASSERTD(dif >= 1);
151 /* now high_dif = 0, dif >= 1, lastc is the neglected part of cp[cn+1] */
153 while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0)
154 && (high_dif == 0) && (dif == 1)))
155 { /* dif=1 implies diff_exp = 0 or 1 */
156 bb = (bn >= 0) ? bp[bn] : 0;
164 else /* diff_exp = 1 */
167 lastc = cp[cn] << (BITS_PER_MP_LIMB - 1);
172 high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1);
175 res += BITS_PER_MP_LIMB;
178 /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
180 if (MPFR_UNLIKELY (high_dif != 0)) /* high_dif == 1 */
189 else /* high_dif == 0 */
193 count_leading_zeros(z, dif); /* dif > 1 here */
195 if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - z - 1))))
196 { /* dif is not a power of two */
202 /* now result is res + (low(b) < low(c)) */
203 while (MPFR_UNLIKELY (bn >= 0 && (cn >= 0 || lastc != 0)))
205 if (diff_exp >= BITS_PER_MP_LIMB)
206 diff_exp -= BITS_PER_MP_LIMB;
212 cc += cp[cn] >> diff_exp;
214 lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
222 *cancel = res + (bp[bn] < cc);
234 while (cn >= 0 && cp[cn] == 0)