[BACK]Return to ref.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / gmp-2.0.2-ssh-2 / mpf / tests

Annotation of OpenXM/src/kan96xx/gmp-2.0.2-ssh-2/mpf/tests/ref.c, Revision 1.1.1.1

1.1       takayama    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: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>