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
;
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 */
135 if (MPFR_LIKELY (diff_exp
< BITS_PER_MP_LIMB
))
137 cc
= cp
[cn
] >> diff_exp
;
138 /* warning: a shift by BITS_PER_MP_LIMB may give wrong results */
140 lastc
= cp
[cn
] << (BITS_PER_MP_LIMB
- diff_exp
);
144 diff_exp
-= BITS_PER_MP_LIMB
;
146 dif
= bp
[bn
--] - cc
; /* necessarily dif >= 1 */
148 while (MPFR_UNLIKELY ((cn
>= 0 || lastc
!= 0)
149 && (high_dif
== 0) && (dif
== 1)))
150 { /* dif=1 implies diff_exp = 0 or 1 */
151 bb
= (bn
>= 0) ? bp
[bn
] : 0;
159 else /* diff_exp = 1 */
162 lastc
= cp
[cn
] << (BITS_PER_MP_LIMB
- 1);
167 high_dif
= 1 - mpn_sub_n (&dif
, &bb
, &cc
, 1);
170 res
+= BITS_PER_MP_LIMB
;
173 /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
175 if (MPFR_UNLIKELY (high_dif
!= 0)) /* high_dif == 1 */
184 else /* high_dif == 0 */
188 count_leading_zeros(z
, dif
); /* dif > 1 here */
190 if (MPFR_LIKELY(dif
!= (MPFR_LIMB_ONE
<< (BITS_PER_MP_LIMB
- z
- 1))))
191 { /* dif is not a power of two */
197 /* now result is res + (low(b) < low(c)) */
198 while (MPFR_UNLIKELY (bn
>= 0 && (cn
>= 0 || lastc
!= 0)))
200 if (diff_exp
>= BITS_PER_MP_LIMB
)
201 diff_exp
-= BITS_PER_MP_LIMB
;
207 cc
+= cp
[cn
] >> diff_exp
;
209 lastc
= cp
[cn
] << (BITS_PER_MP_LIMB
- diff_exp
);
217 *cancel
= res
+ (bp
[bn
] < cc
);
229 while (cn
>= 0 && cp
[cn
] == 0)