1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers.
3 Copyright 1999-2004, 2006-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel 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 3 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.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. */
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 (-1, 0, 1).
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
, mpfr_prec_t
*cancel
)
39 mp_limb_t
*bp
, *cp
, bb
, cc
= 0, lastc
= 0, dif
, high_dif
= 0;
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. */
49 /* the cases b=0 or c=0 are also treated apart in agm and sub
51 MPFR_ASSERTD (MPFR_IS_PURE_FP(b
));
52 MPFR_ASSERTD (MPFR_IS_PURE_FP(c
));
54 if (MPFR_GET_EXP (b
) >= MPFR_GET_EXP (c
))
57 diff_exp
= (mpfr_uexp_t
) MPFR_GET_EXP (b
) - MPFR_GET_EXP (c
);
62 bn
= (MPFR_PREC(b
) - 1) / GMP_NUMB_BITS
;
63 cn
= (MPFR_PREC(c
) - 1) / GMP_NUMB_BITS
; /* # of limbs of c minus 1 */
65 if (MPFR_UNLIKELY( diff_exp
== 0 ))
67 while (bn
>= 0 && cn
>= 0 && bp
[bn
] == cp
[cn
])
74 if (MPFR_UNLIKELY (bn
< 0))
76 if (MPFR_LIKELY (cn
< 0)) /* b = c */
85 if (MPFR_UNLIKELY (cn
< 0))
86 /* c discards exactly the upper part of b */
90 MPFR_ASSERTD (bn
>= 0);
94 if (--bn
< 0) /* b = c */
99 count_leading_zeros(z
, bp
[bn
]); /* bp[bn] <> 0 */
104 MPFR_ASSERTD (bn
>= 0);
105 MPFR_ASSERTD (cn
>= 0);
106 MPFR_ASSERTD (bp
[bn
] != cp
[cn
]);
112 tp
= bp
; bp
= cp
; cp
= tp
;
113 tn
= bn
; bn
= cn
; cn
= tn
;
117 } /* MPFR_EXP(b) >= MPFR_EXP(c) */
118 else /* MPFR_EXP(b) < MPFR_EXP(c) */
121 diff_exp
= (mpfr_uexp_t
) MPFR_GET_EXP (c
) - MPFR_GET_EXP (b
);
126 bn
= (MPFR_PREC(c
) - 1) / GMP_NUMB_BITS
;
127 cn
= (MPFR_PREC(b
) - 1) / GMP_NUMB_BITS
;
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).
136 if (MPFR_LIKELY (diff_exp
< GMP_NUMB_BITS
))
138 cc
= cp
[cn
] >> diff_exp
;
139 /* warning: a shift by GMP_NUMB_BITS may give wrong results */
141 lastc
= cp
[cn
] << (GMP_NUMB_BITS
- diff_exp
);
145 diff_exp
-= GMP_NUMB_BITS
; /* cc = 0 */
147 dif
= bp
[bn
--] - cc
; /* necessarily dif >= 1 */
148 MPFR_ASSERTD(dif
>= 1);
150 /* now high_dif = 0, dif >= 1, lastc is the neglected part of cp[cn+1] */
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;
163 else /* diff_exp = 1 */
166 lastc
= cp
[cn
] << (GMP_NUMB_BITS
- 1);
171 high_dif
= 1 - mpn_sub_n (&dif
, &bb
, &cc
, 1);
174 res
+= GMP_NUMB_BITS
;
177 /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
179 if (MPFR_UNLIKELY (high_dif
!= 0)) /* high_dif == 1 */
182 MPFR_ASSERTD (res
>= 0);
189 else /* high_dif == 0 */
193 count_leading_zeros(z
, dif
); /* dif > 1 here */
195 if (MPFR_LIKELY(dif
!= (MPFR_LIMB_ONE
<< (GMP_NUMB_BITS
- 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
>= GMP_NUMB_BITS
)
206 diff_exp
-= GMP_NUMB_BITS
;
212 cc
+= cp
[cn
] >> diff_exp
;
214 lastc
= cp
[cn
] << (GMP_NUMB_BITS
- diff_exp
);
222 *cancel
= res
+ (bp
[bn
] < cc
);
234 while (cn
>= 0 && cp
[cn
] == 0)