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