contrib/libgmp/mpf/tests/ref.c

   1 /* Reference floating point routines.
   2
   3 Copyright (C) 1996 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 Library General Public License as published by
   9 the Free Software Foundation; either version 2 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 Library General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Library 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., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include "gmp.h"
  23 #include "gmp-impl.h"
  24
  25 #if __STDC__
  26 void ref_mpf_add (mpf_t, const mpf_t, const mpf_t);
  27 void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t);
  28 #else
  29 void ref_mpf_add ();
  30 void ref_mpf_sub ();
  31 #endif
  32
  33 void
  34 #if __STDC__
  35 ref_mpf_add (mpf_t w, const mpf_t u, const mpf_t v)
  36 #else
  37 ref_mpf_add (w, u, v)
  38      mpf_t w;
  39      const mpf_t u;
  40      const mpf_t v;
  41 #endif
  42 {
  43   mp_size_t hi, lo, size;
  44   mp_ptr ut, vt, wt;
  45   int neg;
  46   mp_exp_t exp;
  47   mp_limb_t cy;
  48   TMP_DECL (mark);
  49
  50   TMP_MARK (mark);
  51
  52   if (SIZ (u) == 0)
  53     {
  54       size = ABSIZ (v);
  55       wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB);
  56       MPN_COPY (wt, PTR (v), size);
  57       exp = EXP (v);
  58       neg = SIZ (v) < 0;
  59       goto done;
  60     }
  61   if (SIZ (v) == 0)
  62     {
  63       size = ABSIZ (u);
  64       wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB);
  65       MPN_COPY (wt, PTR (u), size);
  66       exp = EXP (u);
  67       neg = SIZ (u) < 0;
  68       goto done;
  69     }
  70   if ((SIZ (u) ^ SIZ (v)) < 0)
  71     {
  72       mpf_t tmp;
  73       SIZ (tmp) = -SIZ (v);
  74       EXP (tmp) = EXP (v);
  75       PTR (tmp) = PTR (v);
  76       ref_mpf_sub (w, u, tmp);
  77       return;
  78     }
  79   neg = SIZ (u) < 0;
  80
  81   /* Compute the significance of the hi and lo end of the result.  */
  82   hi = MAX (EXP (u), EXP (v));
  83   lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
  84   size = hi - lo;
  85   ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
  86   vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
  87   wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
  88   MPN_ZERO (ut, size);
  89   MPN_ZERO (vt, size);
  90   {int off;
  91   off = size + (EXP (u) - hi) - ABSIZ (u);
  92   MPN_COPY (ut + off, PTR (u), ABSIZ (u));
  93   off = size + (EXP (v) - hi) - ABSIZ (v);
  94   MPN_COPY (vt + off, PTR (v), ABSIZ (v));
  95   }
  96
  97   cy = mpn_add_n (wt, ut, vt, size);
  98   wt[size] = cy;
  99   size += cy;
 100   exp = hi + cy;
 101
 102 done:
 103   if (size > PREC (w))
 104     {
 105       wt += size - PREC (w);
 106       size = PREC (w);
 107     }
 108   MPN_COPY (PTR (w), wt, size);
 109   SIZ (w) = neg == 0 ? size : -size;
 110   EXP (w) = exp;
 111   TMP_FREE (mark);
 112 }
 113
 114 void
 115 #if __STDC__
 116 ref_mpf_sub (mpf_t w, const mpf_t u, const mpf_t v)
 117 #else
 118 ref_mpf_sub (w, u, v)
 119      mpf_t w;
 120      const mpf_t u;
 121      const mpf_t v;
 122 #endif
 123 {
 124   mp_size_t hi, lo, size;
 125   mp_ptr ut, vt, wt;
 126   int neg;
 127   mp_exp_t exp;
 128   TMP_DECL (mark);
 129
 130   TMP_MARK (mark);
 131
 132   if (SIZ (u) == 0)
 133     {
 134       size = ABSIZ (v);
 135       wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB);
 136       MPN_COPY (wt, PTR (v), size);
 137       exp = EXP (v);
 138       neg = SIZ (v) > 0;
 139       goto done;
 140     }
 141   if (SIZ (v) == 0)
 142     {
 143       size = ABSIZ (u);
 144       wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB);
 145       MPN_COPY (wt, PTR (u), size);
 146       exp = EXP (u);
 147       neg = SIZ (u) < 0;
 148       goto done;
 149     }
 150   if ((SIZ (u) ^ SIZ (v)) < 0)
 151     {
 152       mpf_t tmp;
 153       SIZ (tmp) = -SIZ (v);
 154       EXP (tmp) = EXP (v);
 155       PTR (tmp) = PTR (v);
 156       ref_mpf_add (w, u, tmp);
 157       if (SIZ (u) < 0)
 158         mpf_neg (w, w);
 159       return;
 160     }
 161   neg = SIZ (u) < 0;
 162
 163   /* Compute the significance of the hi and lo end of the result.  */
 164   hi = MAX (EXP (u), EXP (v));
 165   lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
 166   size = hi - lo;
 167   ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 168   vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 169   wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
 170   MPN_ZERO (ut, size);
 171   MPN_ZERO (vt, size);
 172   {int off;
 173   off = size + (EXP (u) - hi) - ABSIZ (u);
 174   MPN_COPY (ut + off, PTR (u), ABSIZ (u));
 175   off = size + (EXP (v) - hi) - ABSIZ (v);
 176   MPN_COPY (vt + off, PTR (v), ABSIZ (v));
 177   }
 178
 179   if (mpn_cmp (ut, vt, size) >= 0)
 180     mpn_sub_n (wt, ut, vt, size);
 181   else
 182     {
 183       mpn_sub_n (wt, vt, ut, size);
 184       neg ^= 1;
 185     }
 186   exp = hi;
 187   while (size != 0 && wt[size - 1] == 0)
 188     {
 189       size--;
 190       exp--;
 191     }
 192
 193 done:
 194   if (size > PREC (w))
 195     {
 196       wt += size - PREC (w);
 197       size = PREC (w);
 198     }
 199   MPN_COPY (PTR (w), wt, size);
 200   SIZ (w) = neg == 0 ? size : -size;
 201   EXP (w) = exp;
 202   TMP_FREE (mark);
 203 }