Annotation of OpenXM_contrib/gmp/mpfr/cmp.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpfr_cmp -- compare two floating-point numbers
2:
3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
4:
5: This file is part of the MPFR Library.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
8: it under the terms of the GNU Library General Public License as published by
9: the Free Software Foundation; either version 2 of the License, or (at your
10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Library General Public License
18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include <stdio.h>
23: #include "gmp.h"
24: #include "gmp-impl.h"
25: #include "longlong.h"
26: #include "mpfr.h"
27:
28: /* returns 0 iff b = c
29: a positive value iff b > c
30: a negative value iff b < c
31:
32: More precisely, in case b and c are of same sign, the absolute value
33: of the result is one plus the absolute difference between the exponents
34: of b and c, i.e. one plus the number of bits shifts to align b and c
35: (this value is useful in mpfr_sub).
36:
37: */
38:
39: /* #define DEBUG */
40:
41: /* compares b and sign(s)*c */
42: int
43: #if __STDC__
44: mpfr_cmp3 ( mpfr_srcptr b, mpfr_srcptr c, long int s)
45: #else
46: mpfr_cmp3(b, c, s)
47: mpfr_srcptr b;
48: mpfr_srcptr c;
49: long int s;
50: #endif
51: {
52: long int diff_exp;
53: unsigned long bn, cn;
54: mp_limb_t *bp, *cp;
55:
56: if (!NOTZERO(b) && !NOTZERO(c)) { return 0; }
57: s = s*SIGN(b)*SIGN(c);
58: if (s<0) return(SIGN(b));
59:
60: /* now signs are equal */
61:
62: diff_exp = EXP(b)-EXP(c);
63: s = (SIGN(b)>0) ? 1 : -1;
64:
65: if (diff_exp>0) return(s*(1+diff_exp));
66: else if (diff_exp<0) return(s*(-1+diff_exp));
67: /* both signs and exponents are equal */
68:
69: bn = (PREC(b)-1)/mp_bits_per_limb+1;
70: cn = (PREC(c)-1)/mp_bits_per_limb+1;
71: bp = MANT(b); cp = MANT(c);
72:
73: while (bn && cn) {
74: if (bp[--bn] != cp[--cn])
75: return((bp[bn]>cp[cn]) ? s : -s);
76: }
77:
78: if (bn) { while (bn) if (bp[--bn]) return(s); }
79: else if (cn) while (cn) if (cp[--cn]) return(-s);
80:
81: return 0;
82: }
83:
84: /* returns the number of cancelled bits when one subtracts abs(c) from abs(b).
85: Assumes b>=c, which implies EXP(b)>=EXP(c).
86: if b=c, returns prec(b).
87: */
88: int
89: #if __STDC__
90: mpfr_cmp2 ( mpfr_srcptr b, mpfr_srcptr c )
91: #else
92: mpfr_cmp2(b, c)
93: mpfr_srcptr b;
94: mpfr_srcptr c;
95: #endif
96: {
97: long int d, bn, cn, k, z;
98: mp_limb_t *bp, *cp, t, u=0, cc=0;
99:
100: #ifdef DEBUG
101: printf("b="); mpfr_print_raw(b); putchar('\n');
102: printf("c="); mpfr_print_raw(c); putchar('\n');
103: #endif
104: if (NOTZERO(c)==0) return (NOTZERO(b)) ? 0 : PREC(b);
105: d = EXP(b)-EXP(c);
106: k = 0; /* result can be d or d+1 if d>1, or >= d otherwise */
107: /* k is the number of identical bits in the high part,
108: then z is the number of possibly cancelled bits */
109: #ifdef DEBUG
110: if (d<0) { printf("assumption EXP(b)<EXP(c) violated\n"); exit(1); }
111: #endif
112: bn = (PREC(b)-1)/mp_bits_per_limb;
113: cn = (PREC(c)-1)/mp_bits_per_limb;
114: bp = MANT(b); cp = MANT(c);
115: /* subtract c from b from most significant to less significant limbs,
116: and first determines first non zero limb difference */
117: if (d)
118: {
119: cc = bp[bn--];
120: if (d<mp_bits_per_limb)
121: cc -= cp[cn]>>d;
122: }
123: else { /* d=0 */
124: while (bn>=0 && cn>=0 && (cc=(bp[bn--]-cp[cn--]))==0) {
125: k+=mp_bits_per_limb;
126: }
127:
128: if (cc==0) { /* either bn<0 or cn<0 */
129: while (bn>=0 && (cc=bp[bn--])==0) k+=mp_bits_per_limb;
130: }
131: /* now bn<0 or cc<>0 */
132: if (cc==0 && bn<0) return(PREC(b));
133: }
134:
135: /* the first non-zero limb difference is cc, and the number
136: of cancelled bits in the upper limbs is k */
137:
138: count_leading_zeros(u, cc);
139: k += u;
140:
141: if (cc != (1<<(mp_bits_per_limb-u-1))) return k;
142:
143: /* now cc is an exact power of two */
144: if (cc != 1)
145: /* We just need to compare the following limbs */
146: /* until two of them differ. The result is either k or k+1. */
147: {
148: /* First flush all the unmatched limbs of b ; they all have to
149: be 0 in order for the process to go on */
150: while (bn >= 0)
151: {
152: if (cn < 0) { return k; }
153: t = bp[bn--];
154: if (d < mp_bits_per_limb)
155: {
156: if (d)
157: {
158: u = cp[cn--] << (mp_bits_per_limb - d);
159: if (cn >= 0) u+=(cp[cn]>>d);
160: }
161: else u = cp[cn--];
162:
163: if (t > u || (t == u && cn < 0)) return k;
164: if (t < u) return k+1;
165: }
166: else
167: if (t) return k; else d -= mp_bits_per_limb;
168: }
169:
170: /* bn < 0; if some limb of c is nonzero, return k+1, otherwise return k*/
171:
172: if (cn>=0 && (cp[cn--] << (mp_bits_per_limb - d))) { return k+1; }
173:
174: while (cn >= 0)
175: if (cp[cn--]) return k+1;
176: return k;
177: }
178:
179: /* cc = 1. Too bad. */
180: z = 0; /* number of possibly cancelled bits - 1 */
181: /* thus result is either k if low(b) >= low(c)
182: or k+z+1 if low(b) < low(c) */
183: if (d > mp_bits_per_limb) { return k; }
184:
185: while (bn >= 0)
186: {
187: if (cn < 0) { return k; }
188:
189: if (d)
190: {
191: u = cp[cn--] << (mp_bits_per_limb - d);
192: if (cn >= 0) u+=(cp[cn]>>d);
193: }
194: else u = cp[cn--];
195:
196: /* bp[bn--] > cp[cn--] : no borrow possible, k unchanged
197: bp[bn--] = cp[cn--] : need to consider next limbs
198: bp[bn--] < cp[cn--] : borrow
199: */
200: if ((cc = bp[bn--]) != u) {
201: if (cc>u) return k;
202: else { count_leading_zeros(u, cc-u); return k + 1 + z + u; }
203: }
204: else z += mp_bits_per_limb;
205: }
206:
207: if (cn >= 0)
208: count_leading_zeros(cc, ~(cp[cn--] << (mp_bits_per_limb - d)));
209: else { cc = 0; }
210:
211: k += cc;
212: if (cc < d) return k;
213:
214: while (cn >= 0 && !~cp[cn--]) { z += mp_bits_per_limb; }
215: if (cn >= 0) { count_leading_zeros(cc, ~cp[cn--]); return (k + z + cc); }
216:
217: return k; /* We **need** that the nonsignificant limbs of c are set
218: to zero there */
219: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to cmp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr