1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
/* (c) PolKA project at Inria Lorraine
written by Paul Zimmermann, February 1999
returns 0 iff b = c
a positive value iff b > c
a negative value iff b < c
More precisely, in case b and c are of same sign, the absolute value
of the result is one plus the absolute difference between the exponents
of b and c, i.e. one plus the number of bits shifts to align b and c
(this value is useful in mpfr_sub).
*/
#include <stdio.h>
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
#include "mpfr.h"
int mpfr_cmp ( mpfr_srcptr b, mpfr_srcptr c )
{
long int s, diff_exp;
unsigned long bn, cn;
mp_limb_t *bp, *cp;
s = SIGN(b) * SIGN(c);
if (s<0) return(SIGN(b));
/* now signs are equal */
diff_exp = EXP(b)-EXP(c);
s = (SIGN(b)>0) ? 1 : -1;
if (diff_exp>0) return(s*(1+diff_exp));
else if (diff_exp<0) return(s*(-1+diff_exp));
/* both signs and exponents are equal */
bn = (PREC(b)-1)/mp_bits_per_limb+1;
cn = (PREC(c)-1)/mp_bits_per_limb+1;
bp = MANT(b); cp = MANT(c);
while (bn && cn) {
if (bp[--bn] != cp[--cn])
return((bp[bn]>cp[cn]) ? s : -s);
}
if (bn) { while (bn) if (bp[--bn]) return(s); }
else if (cn) while (cn) if (cp[--cn]) return(-s);
return 0;
}
/* returns the number of cancelled bits when one subtracts abs(c) from abs(b).
Assumes b>=c, which implies EXP(b)>=EXP(c).
if b=c, returns prec(b).
*/
int mpfr_cmp2 ( mpfr_srcptr b, mpfr_srcptr c )
{
long int d, bn, cn, k;
mp_limb_t *bp, *cp, t, u=0, cc=0;
if (NOTZERO(c)==0) return 0;
d = EXP(b)-EXP(c);
k = 0; /* result can be d or d+1 if d>1, or >= d otherwise */
#ifdef DEBUG
if (d<0) { printf("assumption EXP(b)<EXP(c) violated\n"); exit(1); }
#endif
bn = (PREC(b)-1)/mp_bits_per_limb;
cn = (PREC(c)-1)/mp_bits_per_limb;
bp = MANT(b); cp = MANT(c);
/* subtracts c from b from most significant to less significant limbs,
and first determines first non zero limb difference */
if (d) {
cc = bp[bn--];
if (d<mp_bits_per_limb)
cc -= cp[cn]>>d; /* cannot be zero since b is normalized */
}
else { /* d=0 */
while (bn>=0 && cn>=0 && (cc=(bp[bn--]-cp[cn--]))==0) {
k+=mp_bits_per_limb;
}
if (cc==0) { /* either bn<0 or cn<0 */
while (bn>=0 && (cc=bp[bn--])==0) k+=mp_bits_per_limb;
}
/* now bn<0 or cc<>0 */
if (cc==0 && bn<0) return(PREC(b));
}
/* the first non-zero limb difference is cc, and the number
of cancelled bits in the upper limbs is k */
count_leading_zeros(u, cc);
k += u;
if (cc != (1<<(mp_bits_per_limb-u-1))) return k;
/* now cc is an exact power of two */
cc = mp_bits_per_limb-u;
while (bn>=0) {
/* computes limb bn of difference.
cc is previous borrow: either 0 or 1.
*/
t = bp[bn--];
if (d<mp_bits_per_limb) {
if (d) {
u = cp[cn--]<<(mp_bits_per_limb-d);
if (cn>=0) u+=(cp[cn]>>d); else return(k);
}
else u = cp[cn--];
if (t>u) return(k); /* no borrow possible */
t -= u;
if (t>1) { /* t=0 when diff=0, t=1 when b=0 and c=~0 */
count_leading_zeros(u, t);
if (u!=~(mp_limb_t)0) return(k+cc+u); /* borrow = 1 */
else cc+=mp_bits_per_limb;
}
}
else {
d -= mp_bits_per_limb;
}
}
return k;
}
|
