github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/tests/refmpf.c (about)

     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.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  }