Annotation of OpenXM_contrib/gmp/demos/expr/exprfra.c, Revision 1.1
1.1 ! ohara 1: /* mpfr expression evaluation */
! 2:
! 3: /*
! 4: Copyright 2000, 2001 Free Software Foundation, Inc.
! 5:
! 6: This file is part of the GNU MP Library.
! 7:
! 8: The GNU MP Library is free software; you can redistribute it and/or modify
! 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
! 11: option) any later version.
! 12:
! 13: The GNU MP Library is distributed in the hope that it will be useful, but
! 14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Lesser General Public License
! 19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 21: MA 02111-1307, USA.
! 22: */
! 23:
! 24: #include <ctype.h>
! 25: #include <stdio.h>
! 26: #include "gmp.h"
! 27: #include "expr-impl.h"
! 28:
! 29: #ifndef ULONG_MAX
! 30: #define ULONG_MAX (~ (unsigned long) 0)
! 31: #endif
! 32:
! 33:
! 34: /* Change this to "#define TRACE(x) x" to get some traces. */
! 35: #define TRACE(x)
! 36:
! 37:
! 38: static void
! 39: e_mpfr_init2 (mpfr_ptr f, unsigned long prec)
! 40: {
! 41: mpfr_init2 (f, prec);
! 42: }
! 43:
! 44: static void
! 45: e_mpfr_set (mpfr_ptr dst, mpfr_srcptr src)
! 46: {
! 47: mpfr_set (dst, src, ROUND);
! 48: }
! 49:
! 50: /* Test whether fits and is an integer. FIXME: Use mpfr_integer_p and
! 51: mpfr_fits_ulong_p (or just mpfr_cmp_ui), if/when these exist. */
! 52: static int
! 53: e_mpfr_ulong_p (mpfr_srcptr f)
! 54: {
! 55: mp_size_t size = f->_mpfr_size;
! 56: mp_srcptr ptr = f->_mpfr_d;
! 57: mp_exp_t exp = f->_mpfr_exp;
! 58: mp_prec_t high_index = (f->_mpfr_prec-1) / mp_bits_per_limb;
! 59: mp_prec_t i;
! 60: mp_limb_t high = ptr[high_index];
! 61:
! 62: if ((size & 0x60000000) != 0) /* nan or inf don't fit */
! 63: return 0;
! 64:
! 65: if (high == 0) /* zero fits */
! 66: return 1;
! 67:
! 68: for (i = 0; i < high_index; i++) /* low limbs must be zero */
! 69: if (ptr[i] != 0)
! 70: return 0;
! 71:
! 72: if (exp <= 0) /* fractions don't fit */
! 73: return 0;
! 74:
! 75: if (exp > mp_bits_per_limb) /* bigger than a limb doesn't fit */
! 76: return 0;
! 77:
! 78: /* any fraction bits in the high limb must be zero */
! 79: if (exp < mp_bits_per_limb && (high << exp) != 0)
! 80: return 0;
! 81:
! 82: /* value must fit a ulong */
! 83: return (high >> (mp_bits_per_limb - exp)) <= ULONG_MAX;
! 84: }
! 85:
! 86: /* FIXME: Use mpfr_get_ui if/when it exists. */
! 87: static unsigned long
! 88: e_mpfr_get_ui_fits (mpfr_srcptr f)
! 89: {
! 90: mp_srcptr ptr = f->_mpfr_d;
! 91: mp_exp_t exp = f->_mpfr_exp;
! 92: mp_prec_t high_index = (f->_mpfr_prec-1) / mp_bits_per_limb;
! 93: mp_limb_t high = ptr[high_index];
! 94:
! 95: return high >> (mp_bits_per_limb - exp);
! 96: }
! 97:
! 98: static size_t
! 99: e_mpfr_number (mpfr_ptr res, __gmp_const char *e, size_t elen, int base)
! 100: {
! 101: char *edup;
! 102: size_t i, j, ret, extra=0;
! 103:
! 104: TRACE (printf ("mpfr_number prec=%lu, base=%d, \"%.*s\"\n",
! 105: mpfr_get_prec (res), base, (int) elen, e));
! 106:
! 107: /* mpfr_set_str doesn't currently accept 0x for hex in base==0, so do it
! 108: here instead. FIXME: Would prefer to let mpfr_set_str handle this. */
! 109: if (base == 0)
! 110: {
! 111: if (elen >= 2 && e[0] == '0' && (e[1] == 'x' || e[1] == 'X'))
! 112: {
! 113: base = 16;
! 114: extra = 2;
! 115: e += extra;
! 116: elen -= extra;
! 117: }
! 118: else
! 119: base = 10;
! 120: }
! 121:
! 122: i = 0;
! 123: for (;;)
! 124: {
! 125: if (i >= elen)
! 126: goto parsed;
! 127: if (e[i] == '.')
! 128: break;
! 129: if (e[i] == '@' || (base <= 10 && (e[i] == 'e' || e[i] == 'E')))
! 130: goto exponent;
! 131: if (! isasciidigit_in_base (e[i], base == 0 ? 10 : base))
! 132: goto parsed;
! 133: i++;
! 134: }
! 135:
! 136: /* fraction */
! 137: i++;
! 138: for (;;)
! 139: {
! 140: if (i >= elen)
! 141: goto parsed;
! 142: if (e[i] == '@' || (base <= 10 && (e[i] == 'e' || e[i] == 'E')))
! 143: break;
! 144: if (! isasciidigit_in_base (e[i], base == 0 ? 10 : base))
! 145: goto parsed;
! 146: i++;
! 147: }
! 148:
! 149: exponent:
! 150: i++;
! 151: if (i >= elen)
! 152: goto parsed;
! 153: if (e[i] == '-')
! 154: i++;
! 155: for (;;)
! 156: {
! 157: if (i >= elen)
! 158: goto parsed;
! 159: if (! isasciidigit_in_base (e[i], base == 0 ? 10 : base))
! 160: break;
! 161: i++;
! 162: }
! 163:
! 164: parsed:
! 165: TRACE (printf (" parsed i=%d \"%.*s\"\n", i, i, e));
! 166:
! 167: /* mpfr_set_str doesn't currently accept upper case for hex, so convert to
! 168: lower here instead. FIXME: Would prefer to let mpfr_set_str handle
! 169: this. */
! 170: edup = (*__gmp_allocate_func) (i+1);
! 171: for (j = 0; j < i; j++)
! 172: edup[j] = tolower (e[j]);
! 173: edup[i] = '\0';
! 174:
! 175: TRACE (printf (" attempt base=%d, len=%d, \"%s\"\n", base, i, edup));
! 176:
! 177: if (mpfr_set_str (res, edup, base, ROUND) == 0)
! 178: ret = i + extra;
! 179: else
! 180: ret = 0;
! 181:
! 182: (*__gmp_free_func) (edup, i+1);
! 183: return ret;
! 184: }
! 185:
! 186: /* Don't want to change the precision of w, can only do an actual swap when
! 187: w and x have the same precision. */
! 188: static void
! 189: e_mpfr_set_or_swap (mpfr_ptr w, mpfr_ptr x)
! 190: {
! 191: if (mpfr_get_prec (w) == mpfr_get_prec (x))
! 192: mpfr_swap (w, x);
! 193: else
! 194: mpfr_set (w, x, ROUND);
! 195: }
! 196:
! 197: int
! 198: mpfr_expr_a (__gmp_const struct mpexpr_operator_t *table,
! 199: mpfr_ptr res, int base, unsigned long prec,
! 200: __gmp_const char *e, size_t elen,
! 201: mpfr_srcptr var[26])
! 202: {
! 203: struct mpexpr_parse_t p;
! 204:
! 205: p.table = table;
! 206: p.res = (mpX_ptr) res;
! 207: p.base = base;
! 208: p.prec = prec;
! 209: p.e = e;
! 210: p.elen = elen;
! 211: p.var = (mpX_srcptr *) var;
! 212:
! 213: p.mpX_clear = (mpexpr_fun_one_t) mpfr_clear;
! 214: p.mpX_ulong_p = (mpexpr_fun_i_unary_t) e_mpfr_ulong_p;
! 215: p.mpX_get_ui = (mpexpr_fun_get_ui_t) e_mpfr_get_ui_fits;
! 216: p.mpX_init = (mpexpr_fun_unary_ui_t) e_mpfr_init2;
! 217: p.mpX_number = (mpexpr_fun_number_t) e_mpfr_number;
! 218: p.mpX_set = (mpexpr_fun_unary_t) e_mpfr_set;
! 219: p.mpX_set_or_swap = (mpexpr_fun_unary_t) e_mpfr_set_or_swap;
! 220: p.mpX_set_si = (mpexpr_fun_set_si_t) mpfr_set_si;
! 221: p.mpX_swap = (mpexpr_fun_swap_t) mpfr_swap;
! 222:
! 223: return mpexpr_evaluate (&p);
! 224: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>