[BACK]Return to exprfra.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / demos / expr

Annotation of OpenXM_contrib/gmp/demos/expr/exprfra.c, Revision 1.1.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>