1 /* Reference floating point routines.
2
3 Copyright 1996, 2001, 2004, 2005 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., 51 Franklin Street, Fifth Floor, Boston,
20 MA 02110-1301, USA. */
21
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include "mpir.h"
26 #include "gmp-impl.h"
27 #include "tests.h"
28
29
30 void
refmpf_add(mpf_ptr w,mpf_srcptr u,mpf_srcptr v)31 refmpf_add (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
32 {
33 mp_size_t hi, lo, size;
34 mp_ptr ut, vt, wt;
35 int neg;
36 mp_exp_t exp;
37 mp_limb_t cy;
38 TMP_DECL;
39
40 TMP_MARK;
41
42 if (SIZ (u) == 0)
43 {
44 size = ABSIZ (v);
45 wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
46 MPN_COPY (wt, PTR (v), size);
47 exp = EXP (v);
48 neg = SIZ (v) < 0;
49 goto done;
50 }
51 if (SIZ (v) == 0)
52 {
53 size = ABSIZ (u);
54 wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
55 MPN_COPY (wt, PTR (u), size);
56 exp = EXP (u);
57 neg = SIZ (u) < 0;
58 goto done;
59 }
60 if ((SIZ (u) ^ SIZ (v)) < 0)
61 {
62 mpf_t tmp;
63 SIZ (tmp) = -SIZ (v);
64 EXP (tmp) = EXP (v);
65 PTR (tmp) = PTR (v);
66 refmpf_sub (w, u, tmp);
67 return;
68 }
69 neg = SIZ (u) < 0;
70
71 /* Compute the significance of the hi and lo end of the result. */
72 hi = MAX (EXP (u), EXP (v));
73 lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
74 size = hi - lo;
75 ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
76 vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
77 wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
78 MPN_ZERO (ut, size);
79 MPN_ZERO (vt, size);
80 {int off;
81 off = size + (EXP (u) - hi) - ABSIZ (u);
82 MPN_COPY (ut + off, PTR (u), ABSIZ (u));
83 off = size + (EXP (v) - hi) - ABSIZ (v);
84 MPN_COPY (vt + off, PTR (v), ABSIZ (v));
85 }
86
87 cy = mpn_add_n (wt, ut, vt, size);
88 wt[size] = cy;
89 size += cy;
90 exp = hi + cy;
91
92 done:
93 if (size > PREC (w))
94 {
95 wt += size - PREC (w);
96 size = PREC (w);
97 }
98 MPN_COPY (PTR (w), wt, size);
99 SIZ (w) = neg == 0 ? size : -size;
100 EXP (w) = exp;
101 TMP_FREE;
102 }
103
104
105 /* Add 1 "unit in last place" (ie. in the least significant limb) to f.
106 f cannot be zero, since that has no well-defined "last place".
107
108 This routine is designed for use in cases where we pay close attention to
109 the size of the data value and are using that (and the exponent) to
110 indicate the accurate part of a result, or similar. For this reason, if
111 there's a carry out we don't store 1 and adjust the exponent, we just
112 leave 100..00. We don't even adjust if there's a carry out of prec+1
113 limbs, but instead give up in that case (which we intend shouldn't arise
114 in normal circumstances). */
115
116 void
refmpf_add_ulp(mpf_ptr f)117 refmpf_add_ulp (mpf_ptr f)
118 {
119 mp_ptr fp = PTR(f);
120 mp_size_t fsize = SIZ(f);
121 mp_size_t abs_fsize = ABSIZ(f);
122 mp_limb_t c;
123
124 if (fsize == 0)
125 {
126 printf ("Oops, refmpf_add_ulp called with f==0\n");
127 abort ();
128 }
129
130 c = refmpn_add_1 (fp, fp, abs_fsize, CNST_LIMB(1));
131 if (c != 0)
132 {
133 if (abs_fsize >= PREC(f) + 1)
134 {
135 printf ("Oops, refmpf_add_ulp carried out of prec+1 limbs\n");
136 abort ();
137 }
138
139 fp[abs_fsize] = c;
140 abs_fsize++;
141 SIZ(f) = (fsize > 0 ? abs_fsize : - abs_fsize);
142 EXP(f)++;
143 }
144 }
145
146 /* Fill f with size limbs of the given value, setup as an integer. */
147 void
refmpf_fill(mpf_ptr f,mp_size_t size,mp_limb_t value)148 refmpf_fill (mpf_ptr f, mp_size_t size, mp_limb_t value)
149 {
150 ASSERT (size >= 0);
151 size = MIN (PREC(f) + 1, size);
152 SIZ(f) = size;
153 EXP(f) = size;
154 refmpn_fill (PTR(f), size, value);
155 }
156
157 /* Strip high zero limbs from the f data, adjusting exponent accordingly. */
158 void
refmpf_normalize(mpf_ptr f)159 refmpf_normalize (mpf_ptr f)
160 {
161 while (SIZ(f) != 0 && PTR(f)[ABSIZ(f)-1] == 0)
162 {
163 SIZ(f) = (SIZ(f) >= 0 ? SIZ(f)-1 : SIZ(f)+1);
164 EXP(f) --;
165 }
166 if (SIZ(f) == 0)
167 EXP(f) = 0;
168 }
169
170 /* refmpf_set_overlap sets up dst as a copy of src, but with PREC(dst)
171 unchanged, in preparation for an overlap test.
172
173 The full value of src is copied, and the space at PTR(dst) is extended as
174 necessary. The way PREC(dst) is unchanged is as per an mpf_set_prec_raw.
175 The return value is the new PTR(dst) space precision, in bits, ready for
176 a restoring mpf_set_prec_raw before mpf_clear. */
177
178 unsigned long
refmpf_set_overlap(mpf_ptr dst,mpf_srcptr src)179 refmpf_set_overlap (mpf_ptr dst, mpf_srcptr src)
180 {
181 mp_size_t dprec = PREC(dst);
182 mp_size_t ssize = ABSIZ(src);
183 unsigned long ret;
184
185 refmpf_set_prec_limbs (dst, (unsigned long) MAX (dprec, ssize));
186 mpf_set (dst, src);
187
188 ret = mpf_get_prec (dst);
189 PREC(dst) = dprec;
190 return ret;
191 }
192
193 /* Like mpf_set_prec, but taking a precision in limbs.
194 PREC(f) ends up as the given "prec" value. */
195 void
refmpf_set_prec_limbs(mpf_ptr f,unsigned long prec)196 refmpf_set_prec_limbs (mpf_ptr f, unsigned long prec)
197 {
198 mpf_set_prec (f, __GMPF_PREC_TO_BITS (prec));
199 }
200
201
202 void
refmpf_sub(mpf_ptr w,mpf_srcptr u,mpf_srcptr v)203 refmpf_sub (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
204 {
205 mp_size_t hi, lo, size;
206 mp_ptr ut, vt, wt;
207 int neg;
208 mp_exp_t exp;
209 TMP_DECL;
210
211 TMP_MARK;
212
213 if (SIZ (u) == 0)
214 {
215 size = ABSIZ (v);
216 wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
217 MPN_COPY (wt, PTR (v), size);
218 exp = EXP (v);
219 neg = SIZ (v) > 0;
220 goto done;
221 }
222 if (SIZ (v) == 0)
223 {
224 size = ABSIZ (u);
225 wt = (mp_ptr) TMP_ALLOC ((size+1) * BYTES_PER_MP_LIMB);
226 MPN_COPY (wt, PTR (u), size);
227 exp = EXP (u);
228 neg = SIZ (u) < 0;
229 goto done;
230 }
231 if ((SIZ (u) ^ SIZ (v)) < 0)
232 {
233 mpf_t tmp;
234 SIZ (tmp) = -SIZ (v);
235 EXP (tmp) = EXP (v);
236 PTR (tmp) = PTR (v);
237 refmpf_add (w, u, tmp);
238 if (SIZ (u) < 0)
239 mpf_neg (w, w);
240 return;
241 }
242 neg = SIZ (u) < 0;
243
244 /* Compute the significance of the hi and lo end of the result. */
245 hi = MAX (EXP (u), EXP (v));
246 lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
247 size = hi - lo;
248 ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
249 vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
250 wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
251 MPN_ZERO (ut, size);
252 MPN_ZERO (vt, size);
253 {int off;
254 off = size + (EXP (u) - hi) - ABSIZ (u);
255 MPN_COPY (ut + off, PTR (u), ABSIZ (u));
256 off = size + (EXP (v) - hi) - ABSIZ (v);
257 MPN_COPY (vt + off, PTR (v), ABSIZ (v));
258 }
259
260 if (mpn_cmp (ut, vt, size) >= 0)
261 mpn_sub_n (wt, ut, vt, size);
262 else
263 {
264 mpn_sub_n (wt, vt, ut, size);
265 neg ^= 1;
266 }
267 exp = hi;
268 while (size != 0 && wt[size - 1] == 0)
269 {
270 size--;
271 exp--;
272 }
273
274 done:
275 if (size > PREC (w))
276 {
277 wt += size - PREC (w);
278 size = PREC (w);
279 }
280 MPN_COPY (PTR (w), wt, size);
281 SIZ (w) = neg == 0 ? size : -size;
282 EXP (w) = exp;
283 TMP_FREE;
284 }
285
286
287 /* Validate got by comparing to want. Return 1 if good, 0 if bad.
288
289 The data in got is compared to that in want, up to either PREC(got) limbs
290 or the size of got, whichever is bigger. Clearly we always demand
291 PREC(got) of accuracy, but we go further and say that if got is bigger
292 then any extra must be correct too.
293
294 want needs to have enough data to allow this comparison. The size in
295 want doesn't have to be that big though, if it's smaller then further low
296 limbs are taken to be zero.
297
298 This validation approach is designed to allow some flexibility in exactly
299 how much data is generated by an mpf function, ie. either prec or prec+1
300 limbs. We don't try to make a reference function that emulates that same
301 size decision, instead the idea is for a validation function to generate
302 at least as much data as the real function, then compare. */
303
304 int
refmpf_validate(const char * name,mpf_srcptr got,mpf_srcptr want)305 refmpf_validate (const char *name, mpf_srcptr got, mpf_srcptr want)
306 {
307 int bad = 0;
308 mp_size_t gsize, wsize, cmpsize, i;
309 mp_srcptr gp, wp;
310 mp_limb_t glimb, wlimb;
311
312 MPF_CHECK_FORMAT (got);
313
314 if (EXP (got) != EXP (want))
315 {
316 printf ("%s: wrong exponent\n", name);
317 bad = 1;
318 }
319
320 gsize = SIZ (got);
321 wsize = SIZ (want);
322 if ((gsize < 0 && wsize > 0) || (gsize > 0 && wsize < 0))
323 {
324 printf ("%s: wrong sign\n", name);
325 bad = 1;
326 }
327
328 gsize = ABS (gsize);
329 wsize = ABS (wsize);
330
331 /* most significant limb of respective data */
332 gp = PTR (got) + gsize - 1;
333 wp = PTR (want) + wsize - 1;
334
335 /* compare limb data */
336 cmpsize = MAX (PREC (got), gsize);
337 for (i = 0; i < cmpsize; i++)
338 {
339 glimb = (i < gsize ? gp[-i] : 0);
340 wlimb = (i < wsize ? wp[-i] : 0);
341
342 if (glimb != wlimb)
343 {
344 printf ("%s: wrong data starting at index %ld from top\n",
345 name, (long) i);
346 bad = 1;
347 break;
348 }
349 }
350
351 if (bad)
352 {
353 printf (" prec %d\n", PREC(got));
354 printf (" exp got %ld\n", (long) EXP(got));
355 printf (" exp want %ld\n", (long) EXP(want));
356 printf (" size got %d\n", SIZ(got));
357 printf (" size want %d\n", SIZ(want));
358 printf (" limbs (high to low)\n");
359 printf (" got ");
360 for (i = ABSIZ(got)-1; i >= 0; i--)
361 {
362 gmp_printf ("%MX", PTR(got)[i]);
363 if (i != 0)
364 printf (",");
365 }
366 printf ("\n");
367 printf (" want ");
368 for (i = ABSIZ(want)-1; i >= 0; i--)
369 {
370 gmp_printf ("%MX", PTR(want)[i]);
371 if (i != 0)
372 printf (",");
373 }
374 printf ("\n");
375 return 0;
376 }
377
378 return 1;
379 }
380
381
382 int
refmpf_validate_division(const char * name,mpf_srcptr got,mpf_srcptr n,mpf_srcptr d)383 refmpf_validate_division (const char *name, mpf_srcptr got,
384 mpf_srcptr n, mpf_srcptr d)
385 {
386 mp_size_t nsize, dsize, sign, prec, qsize, tsize;
387 mp_srcptr np, dp;
388 mp_ptr tp, qp, rp;
389 mpf_t want;
390 int ret;
391
392 nsize = SIZ (n);
393 dsize = SIZ (d);
394 ASSERT_ALWAYS (dsize != 0);
395
396 sign = nsize ^ dsize;
397 nsize = ABS (nsize);
398 dsize = ABS (dsize);
399
400 np = PTR (n);
401 dp = PTR (d);
402 prec = PREC (got);
403
404 EXP (want) = EXP (n) - EXP (d) + 1;
405
406 qsize = prec + 2; /* at least prec+1 limbs, after high zero */
407 tsize = qsize + dsize - 1; /* dividend size to give desired qsize */
408
409 /* dividend n, extended or truncated */
410 tp = refmpn_malloc_limbs (tsize);
411 refmpn_copy_extend (tp, tsize, np, nsize);
412
413 qp = refmpn_malloc_limbs (qsize);
414 rp = refmpn_malloc_limbs (dsize); /* remainder, unused */
415
416 ASSERT_ALWAYS (qsize == tsize - dsize + 1);
417 refmpn_tdiv_qr (qp, rp, (mp_size_t) 0, tp, tsize, dp, dsize);
418
419 PTR (want) = qp;
420 SIZ (want) = (sign >= 0 ? qsize : -qsize);
421 refmpf_normalize (want);
422
423 ret = refmpf_validate (name, got, want);
424
425 free (tp);
426 free (qp);
427 free (rp);
428
429 return ret;
430 }
431