diff options
Diffstat (limited to 'contrib/libgmp/mpf/tests/ref.c')
| -rw-r--r-- | contrib/libgmp/mpf/tests/ref.c | 203 |
1 files changed, 203 insertions, 0 deletions
diff --git a/contrib/libgmp/mpf/tests/ref.c b/contrib/libgmp/mpf/tests/ref.c new file mode 100644 index 0000000..b2c3e6c --- /dev/null +++ b/contrib/libgmp/mpf/tests/ref.c @@ -0,0 +1,203 @@ +/* Reference floating point routines. + +Copyright (C) 1996 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +The GNU MP Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +#if __STDC__ +void ref_mpf_add (mpf_t, const mpf_t, const mpf_t); +void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t); +#else +void ref_mpf_add (); +void ref_mpf_sub (); +#endif + +void +#if __STDC__ +ref_mpf_add (mpf_t w, const mpf_t u, const mpf_t v) +#else +ref_mpf_add (w, u, v) + mpf_t w; + const mpf_t u; + const mpf_t v; +#endif +{ + mp_size_t hi, lo, size; + mp_ptr ut, vt, wt; + int neg; + mp_exp_t exp; + mp_limb_t cy; + TMP_DECL (mark); + + TMP_MARK (mark); + + if (SIZ (u) == 0) + { + size = ABSIZ (v); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (v), size); + exp = EXP (v); + neg = SIZ (v) < 0; + goto done; + } + if (SIZ (v) == 0) + { + size = ABSIZ (u); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (u), size); + exp = EXP (u); + neg = SIZ (u) < 0; + goto done; + } + if ((SIZ (u) ^ SIZ (v)) < 0) + { + mpf_t tmp; + SIZ (tmp) = -SIZ (v); + EXP (tmp) = EXP (v); + PTR (tmp) = PTR (v); + ref_mpf_sub (w, u, tmp); + return; + } + neg = SIZ (u) < 0; + + /* Compute the significance of the hi and lo end of the result. */ + hi = MAX (EXP (u), EXP (v)); + lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v)); + size = hi - lo; + ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + MPN_ZERO (ut, size); + MPN_ZERO (vt, size); + {int off; + off = size + (EXP (u) - hi) - ABSIZ (u); + MPN_COPY (ut + off, PTR (u), ABSIZ (u)); + off = size + (EXP (v) - hi) - ABSIZ (v); + MPN_COPY (vt + off, PTR (v), ABSIZ (v)); + } + + cy = mpn_add_n (wt, ut, vt, size); + wt[size] = cy; + size += cy; + exp = hi + cy; + +done: + if (size > PREC (w)) + { + wt += size - PREC (w); + size = PREC (w); + } + MPN_COPY (PTR (w), wt, size); + SIZ (w) = neg == 0 ? size : -size; + EXP (w) = exp; + TMP_FREE (mark); +} + +void +#if __STDC__ +ref_mpf_sub (mpf_t w, const mpf_t u, const mpf_t v) +#else +ref_mpf_sub (w, u, v) + mpf_t w; + const mpf_t u; + const mpf_t v; +#endif +{ + mp_size_t hi, lo, size; + mp_ptr ut, vt, wt; + int neg; + mp_exp_t exp; + TMP_DECL (mark); + + TMP_MARK (mark); + + if (SIZ (u) == 0) + { + size = ABSIZ (v); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (v), size); + exp = EXP (v); + neg = SIZ (v) > 0; + goto done; + } + if (SIZ (v) == 0) + { + size = ABSIZ (u); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (u), size); + exp = EXP (u); + neg = SIZ (u) < 0; + goto done; + } + if ((SIZ (u) ^ SIZ (v)) < 0) + { + mpf_t tmp; + SIZ (tmp) = -SIZ (v); + EXP (tmp) = EXP (v); + PTR (tmp) = PTR (v); + ref_mpf_add (w, u, tmp); + if (SIZ (u) < 0) + mpf_neg (w, w); + return; + } + neg = SIZ (u) < 0; + + /* Compute the significance of the hi and lo end of the result. */ + hi = MAX (EXP (u), EXP (v)); + lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v)); + size = hi - lo; + ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + MPN_ZERO (ut, size); + MPN_ZERO (vt, size); + {int off; + off = size + (EXP (u) - hi) - ABSIZ (u); + MPN_COPY (ut + off, PTR (u), ABSIZ (u)); + off = size + (EXP (v) - hi) - ABSIZ (v); + MPN_COPY (vt + off, PTR (v), ABSIZ (v)); + } + + if (mpn_cmp (ut, vt, size) >= 0) + mpn_sub_n (wt, ut, vt, size); + else + { + mpn_sub_n (wt, vt, ut, size); + neg ^= 1; + } + exp = hi; + while (size != 0 && wt[size - 1] == 0) + { + size--; + exp--; + } + +done: + if (size > PREC (w)) + { + wt += size - PREC (w); + size = PREC (w); + } + MPN_COPY (PTR (w), wt, size); + SIZ (w) = neg == 0 ? size : -size; + EXP (w) = exp; + TMP_FREE (mark); +} |
