![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to cmp2.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr|
Return to cmp2.c CVS log | Up to [local] / OpenXM_contrib / gmp / mpfr |
1.1 ohara 1: /* mpfr_cmp2 -- exponent shift when subtracting two numbers.
2:
3: Copyright 1999, 2000, 2001 Free Software Foundation.
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 Lesser General Public License as published by
9: the Free Software Foundation; either version 2.1 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 Lesser General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Lesser 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 "gmp.h"
23: #include "gmp-impl.h"
24: #include "longlong.h"
25: #include "mpfr.h"
26: #include "mpfr-impl.h"
27:
28: /* Returns the number of canceled bits when one subtracts |c| from |b|
29: if |b| != |c|, and the sign.
30:
31: Assumes neither of b or c is NaN or +/- infinity.
32:
33: In other terms, if |b| != |c|, mpfr_cmp2 (b, c) returns
34: EXP(max(|b|,|c|)) - EXP(|b| - |c|).
35: */
36:
37: int
38: mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mp_prec_t *cancel)
39: {
40: mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0;
41: mp_size_t bn, cn;
42: mp_exp_unsigned_t diff_exp;
43: mp_prec_t res = 0;
44: int sign;
45:
46: MPFR_ASSERTN(MPFR_IS_FP(b));
47: MPFR_ASSERTN(MPFR_IS_FP(c));
48:
49: /* Optimized case x - x */
50: if (b == c)
51: return 0;
52:
53: if (MPFR_IS_ZERO(b))
54: {
55: if (MPFR_IS_ZERO(c))
56: return 0;
57:
58: *cancel = 0;
59: return -1;
60: }
61:
62: if (MPFR_IS_ZERO(c))
63: {
64: *cancel = 0;
65: return 1;
66: }
67:
68: if (MPFR_EXP(b) >= MPFR_EXP(c))
69: {
70: sign = 1;
71: diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c);
72:
73: bp = MPFR_MANT(b);
74: cp = MPFR_MANT(c);
75:
76: bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
77: cn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
78:
79: if (diff_exp == 0)
80: {
81: while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn])
82: {
83: bn--;
84: cn--;
85: res += BITS_PER_MP_LIMB;
86: }
87:
88: if (bn < 0)
89: {
90: if (cn < 0) /* b = c */
91: return 0;
92:
93: bp = cp;
94: bn = cn;
95: cn = -1;
96: sign = -1;
97: }
98:
99: if (cn < 0) /* c discards exactly the upper part of b */
100: {
101: unsigned int z;
102:
103: MPFR_ASSERTN(bn >= 0);
104:
105: while (bp[bn] == 0)
106: {
107: if (--bn < 0) /* b = c */
108: return 0;
109: res += BITS_PER_MP_LIMB;
110: }
111:
112: count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
113: *cancel = res + z;
114: return sign;
115: }
116:
117: MPFR_ASSERTN(bn >= 0);
118: MPFR_ASSERTN(cn >= 0);
119: MPFR_ASSERTN(bp[bn] != cp[cn]);
120: if (bp[bn] < cp[cn])
121: {
122: mp_limb_t *tp;
123: mp_size_t tn;
124:
125: tp = bp; bp = cp; cp = tp;
126: tn = bn; bn = cn; cn = tn;
127: sign = -1;
128: }
129: }
130: } /* MPFR_EXP(b) >= MPFR_EXP(c) */
131: else /* MPFR_EXP(b) < MPFR_EXP(c) */
132: {
133: sign = -1;
134: diff_exp = (mp_exp_unsigned_t) MPFR_EXP(c) - MPFR_EXP(b);
135:
136: bp = MPFR_MANT(c);
137: cp = MPFR_MANT(b);
138:
139: bn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
140: cn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
141: }
142:
143: /* now we have removed the identical upper limbs of b and c
144: (can happen only when diff_exp = 0), and after the possible
145: swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0 */
146:
147: if (diff_exp < BITS_PER_MP_LIMB)
148: {
149: cc = cp[cn] >> diff_exp;
150: /* warning: a shift by BITS_PER_MP_LIMB may give wrong results */
151: if (diff_exp)
152: lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
153: cn--;
154: }
155: else
156: diff_exp -= BITS_PER_MP_LIMB;
157:
158: dif = bp[bn--] - cc; /* necessarily dif >= 1 */
159:
160: while ((cn >= 0 || lastc != 0) && (high_dif == 0) && (dif == 1))
161: { /* dif=1 implies diff_exp = 0 or 1 */
162: bb = (bn >= 0) ? bp[bn] : 0;
163: cc = lastc;
164: if (cn >= 0)
165: {
166: if (diff_exp == 0)
167: {
168: cc += cp[cn];
169: }
170: else /* diff_exp = 1 */
171: {
172: cc += cp[cn] >> 1;
173: lastc = cp[cn] << (BITS_PER_MP_LIMB - 1);
174: }
175: }
176: else
177: lastc = 0;
178: high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1);
179: bn--;
180: cn--;
181: res += BITS_PER_MP_LIMB;
182: }
183:
184: /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
185:
186: if (high_dif != 0) /* necessarily high_dif = 1 */
187: {
188: res--;
189: if (dif != 0)
190: {
191: *cancel = res;
192: return sign;
193: }
194: }
195: else /* high_dif = 0 */
196: {
197: unsigned int z;
198:
199: count_leading_zeros(z, dif); /* dif > 1 here */
200: res += z;
201: if (dif != (MP_LIMB_T_ONE << (BITS_PER_MP_LIMB - z - 1)))
202: { /* dif is not a power of two */
203: *cancel = res;
204: return sign;
205: }
206: }
207:
208: /* now result is res + (low(b) < low(c)) */
209: while (bn >= 0 && (cn >= 0 || lastc != 0))
210: {
211: if (diff_exp >= BITS_PER_MP_LIMB)
212: diff_exp -= BITS_PER_MP_LIMB;
213: else
214: {
215: cc = lastc;
216: if (cn >= 0)
217: {
218: cc += cp[cn] >> diff_exp;
219: if (diff_exp != 0)
220: lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
221: }
222: else
223: lastc = 0;
224: cn--;
225: }
226: if (bp[bn] != cc)
227: {
228: *cancel = res + (bp[bn] < cc);
229: return sign;
230: }
231: bn--;
232: }
233:
234: if (bn < 0)
235: {
236: if (lastc != 0)
237: res++;
238: else
239: {
240: while (cn >= 0 && cp[cn] == 0)
241: cn--;
242: if (cn >= 0)
243: res++;
244: }
245: }
246:
247: *cancel = res;
248: return sign;
249: }