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