Annotation of OpenXM_contrib/gmp/tests/refmpf.c, Revision 1.1.1.1
1.1 ohara 1: /* Reference floating point routines.
2:
3: Copyright 1996, 2001 Free Software Foundation, Inc.
4:
5: This file is part of the GNU MP Library.
6:
7: The GNU MP 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 GNU MP 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 GNU MP 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 "tests.h"
25:
26:
27: void
28: refmpf_add (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
29: {
30: mp_size_t hi, lo, size;
31: mp_ptr ut, vt, wt;
32: int neg;
33: mp_exp_t exp;
34: mp_limb_t cy;
35: TMP_DECL (mark);
36:
37: TMP_MARK (mark);
38:
39: if (SIZ (u) == 0)
40: {
41: size = ABSIZ (v);
42: wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
43: MPN_COPY (wt, PTR (v), size);
44: exp = EXP (v);
45: neg = SIZ (v) < 0;
46: goto done;
47: }
48: if (SIZ (v) == 0)
49: {
50: size = ABSIZ (u);
51: wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
52: MPN_COPY (wt, PTR (u), size);
53: exp = EXP (u);
54: neg = SIZ (u) < 0;
55: goto done;
56: }
57: if ((SIZ (u) ^ SIZ (v)) < 0)
58: {
59: mpf_t tmp;
60: SIZ (tmp) = -SIZ (v);
61: EXP (tmp) = EXP (v);
62: PTR (tmp) = PTR (v);
63: refmpf_sub (w, u, tmp);
64: return;
65: }
66: neg = SIZ (u) < 0;
67:
68: /* Compute the significance of the hi and lo end of the result. */
69: hi = MAX (EXP (u), EXP (v));
70: lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
71: size = hi - lo;
72: ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
73: vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
74: wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
75: MPN_ZERO (ut, size);
76: MPN_ZERO (vt, size);
77: {int off;
78: off = size + (EXP (u) - hi) - ABSIZ (u);
79: MPN_COPY (ut + off, PTR (u), ABSIZ (u));
80: off = size + (EXP (v) - hi) - ABSIZ (v);
81: MPN_COPY (vt + off, PTR (v), ABSIZ (v));
82: }
83:
84: cy = mpn_add_n (wt, ut, vt, size);
85: wt[size] = cy;
86: size += cy;
87: exp = hi + cy;
88:
89: done:
90: if (size > PREC (w))
91: {
92: wt += size - PREC (w);
93: size = PREC (w);
94: }
95: MPN_COPY (PTR (w), wt, size);
96: SIZ (w) = neg == 0 ? size : -size;
97: EXP (w) = exp;
98: TMP_FREE (mark);
99: }
100:
101: void
102: refmpf_sub (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
103: {
104: mp_size_t hi, lo, size;
105: mp_ptr ut, vt, wt;
106: int neg;
107: mp_exp_t exp;
108: TMP_DECL (mark);
109:
110: TMP_MARK (mark);
111:
112: if (SIZ (u) == 0)
113: {
114: size = ABSIZ (v);
115: wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
116: MPN_COPY (wt, PTR (v), size);
117: exp = EXP (v);
118: neg = SIZ (v) > 0;
119: goto done;
120: }
121: if (SIZ (v) == 0)
122: {
123: size = ABSIZ (u);
124: wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
125: MPN_COPY (wt, PTR (u), size);
126: exp = EXP (u);
127: neg = SIZ (u) < 0;
128: goto done;
129: }
130: if ((SIZ (u) ^ SIZ (v)) < 0)
131: {
132: mpf_t tmp;
133: SIZ (tmp) = -SIZ (v);
134: EXP (tmp) = EXP (v);
135: PTR (tmp) = PTR (v);
136: refmpf_add (w, u, tmp);
137: if (SIZ (u) < 0)
138: mpf_neg (w, w);
139: return;
140: }
141: neg = SIZ (u) < 0;
142:
143: /* Compute the significance of the hi and lo end of the result. */
144: hi = MAX (EXP (u), EXP (v));
145: lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
146: size = hi - lo;
147: ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
148: vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
149: wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
150: MPN_ZERO (ut, size);
151: MPN_ZERO (vt, size);
152: {int off;
153: off = size + (EXP (u) - hi) - ABSIZ (u);
154: MPN_COPY (ut + off, PTR (u), ABSIZ (u));
155: off = size + (EXP (v) - hi) - ABSIZ (v);
156: MPN_COPY (vt + off, PTR (v), ABSIZ (v));
157: }
158:
159: if (mpn_cmp (ut, vt, size) >= 0)
160: mpn_sub_n (wt, ut, vt, size);
161: else
162: {
163: mpn_sub_n (wt, vt, ut, size);
164: neg ^= 1;
165: }
166: exp = hi;
167: while (size != 0 && wt[size - 1] == 0)
168: {
169: size--;
170: exp--;
171: }
172:
173: done:
174: if (size > PREC (w))
175: {
176: wt += size - PREC (w);
177: size = PREC (w);
178: }
179: MPN_COPY (PTR (w), wt, size);
180: SIZ (w) = neg == 0 ? size : -size;
181: EXP (w) = exp;
182: TMP_FREE (mark);
183: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to refmpf.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / tests