1/* Reference floating point routines.
2
3Copyright 1996, 2001, 2004, 2005 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library test suite.
6
7The GNU MP Library test suite is free software; you can redistribute it
8and/or modify it under the terms of the GNU General Public License as
9published by the Free Software Foundation; either version 3 of the License,
10or (at your option) any later version.
11
12The GNU MP Library test suite is distributed in the hope that it will be
13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15Public License for more details.
16
17You should have received a copy of the GNU General Public License along with
18the 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
27void
28refmpf_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
89done:
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
113void
114refmpf_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. */
144void
145refmpf_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. */
155void
156refmpf_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
175unsigned long
176refmpf_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.  */
192void
193refmpf_set_prec_limbs (mpf_ptr f, unsigned long prec)
194{
195  mpf_set_prec (f, __GMPF_PREC_TO_BITS (prec));
196}
197
198
199void
200refmpf_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
271done:
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
301int
302refmpf_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
379int
380refmpf_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}
428