tests/refmpf.c

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