Annotation of OpenXM_contrib/gmp/mpfr/exp.c, Revision 1.1
1.1 ! maekawa 1: /* mpfr_exp -- exponential of a floating-point number
! 2:
! 3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
! 4:
! 5: This file is part of the MPFR Library.
! 6:
! 7: The MPFR 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 MPFR 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 MPFR 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 <stdio.h>
! 23: #include <math.h>
! 24: #include "gmp.h"
! 25: #include "gmp-impl.h"
! 26: #include "mpfr.h"
! 27:
! 28: /* #define DEBUG */
! 29:
! 30: #define LOG2 0.69314718055994528622 /* log(2) rounded to zero on 53 bits */
! 31:
! 32: /* use Brent's formula exp(x) = (1+r+r^2/2!+r^3/3!+...)^(2^K)*2^n
! 33: where x = n*log(2)+(2^K)*r
! 34: number of operations = O(K+prec(r)/K)
! 35: */
! 36: int
! 37: #if __STDC__
! 38: mpfr_exp(mpfr_ptr y, mpfr_srcptr x, unsigned char rnd_mode)
! 39: #else
! 40: mpfr_exp(y, x, rnd_mode)
! 41: mpfr_ptr y;
! 42: mpfr_srcptr x;
! 43: unsigned char rnd_mode;
! 44: #endif
! 45: {
! 46: int n, expx, K, precy, q, k, l, expr, err;
! 47: mpfr_t r, s, t;
! 48:
! 49: if (FLAG_NAN(x)) { SET_NAN(y); return 1; }
! 50: if (!NOTZERO(x)) { mpfr_set_ui(y, 1, GMP_RNDN); return 0; }
! 51:
! 52: expx = EXP(x);
! 53: precy = PREC(y);
! 54: #ifdef DEBUG
! 55: printf("EXP(x)=%d\n",expx);
! 56: #endif
! 57:
! 58: /* if x > (2^31-1)*ln(2), then exp(x) > 2^(2^31-1) i.e. gives +infinity */
! 59: if (expx > 30) {
! 60: if (SIGN(x)>0) { printf("+infinity"); return 1; }
! 61: else { SET_ZERO(y); return 1; }
! 62: }
! 63:
! 64: /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */
! 65: if (expx < -precy) { int signx = SIGN(x);
! 66: mpfr_set_ui(y, 1, rnd_mode);
! 67: if (signx>0 && rnd_mode==GMP_RNDU) mpfr_add_one_ulp(y);
! 68: else if (signx<0 && (rnd_mode==GMP_RNDD || rnd_mode==GMP_RNDZ))
! 69: mpfr_sub_one_ulp(y);
! 70: return 1; }
! 71:
! 72: n = (int) floor(mpfr_get_d(x)/LOG2);
! 73:
! 74: K = (int) sqrt( (double) precy );
! 75: l = (precy-1)/K + 1;
! 76: err = K + (int) ceil(log(2.0*(double)l+18.0)/LOG2);
! 77: /* add K extra bits, i.e. failure probability <= 1/2^K = O(1/precy) */
! 78: q = precy + err + K + 3;
! 79: mpfr_init2(r, q); mpfr_init2(s, q); mpfr_init2(t, q);
! 80: /* the algorithm consists in computing an upper bound of exp(x) using
! 81: a precision of q bits, and see if we can round to PREC(y) taking
! 82: into account the maximal error. Otherwise we increase q. */
! 83: do {
! 84: #ifdef DEBUG
! 85: printf("n=%d K=%d l=%d q=%d\n",n,K,l,q);
! 86: #endif
! 87:
! 88: /* if n<0, we have to get an upper bound of log(2)
! 89: in order to get an upper bound of r = x-n*log(2) */
! 90: mpfr_log2(s, (n>=0) ? GMP_RNDZ : GMP_RNDU);
! 91: #ifdef DEBUG
! 92: printf("n=%d log(2)=",n); mpfr_print_raw(s); putchar('\n');
! 93: #endif
! 94: mpfr_mul_ui(r, s, (n<0) ? -n : n, (n>=0) ? GMP_RNDZ : GMP_RNDU);
! 95: if (n<0) mpfr_neg(r, r, GMP_RNDD);
! 96: /* r = floor(n*log(2)) */
! 97:
! 98: #ifdef DEBUG
! 99: printf("x=%1.20e\n",mpfr_get_d(x));
! 100: printf(" ="); mpfr_print_raw(x); putchar('\n');
! 101: printf("r=%1.20e\n",mpfr_get_d(r));
! 102: printf(" ="); mpfr_print_raw(r); putchar('\n');
! 103: #endif
! 104: mpfr_sub(r, x, r, GMP_RNDU);
! 105: if (SIGN(r)<0) { /* initial approximation n was too large */
! 106: n--;
! 107: mpfr_mul_ui(r, s, (n<0) ? -n : n, GMP_RNDZ);
! 108: if (n<0) mpfr_neg(r, r, GMP_RNDD);
! 109: mpfr_sub(r, x, r, GMP_RNDU);
! 110: }
! 111: #ifdef DEBUG
! 112: printf("x-r=%1.20e\n",mpfr_get_d(r));
! 113: printf(" ="); mpfr_print_raw(r); putchar('\n');
! 114: if (SIGN(r)<0) { fprintf(stderr,"Error in mpfr_exp: r<0\n"); exit(1); }
! 115: #endif
! 116: mpfr_div_2exp(r, r, K, GMP_RNDU); /* r = (x-n*log(2))/2^K */
! 117: mpfr_set_ui(s, 1, GMP_RNDU);
! 118: mpfr_set_ui(t, 1, GMP_RNDU);
! 119:
! 120: l = 1; expr = EXP(r);
! 121: do {
! 122: mpfr_mul(t, t, r, GMP_RNDU);
! 123: mpfr_div_ui(t, t, l, GMP_RNDU);
! 124: mpfr_add(s, s, t, GMP_RNDU);
! 125: #ifdef DEBUG
! 126: printf("l=%d t=%1.20e\n",l,mpfr_get_d(t));
! 127: printf("s=%1.20e\n",mpfr_get_d(s));
! 128: #endif
! 129: l++;
! 130: } while (EXP(t)+expr > -q);
! 131: #ifdef DEBUG
! 132: printf("l=%d q=%d (K+l)*q^2=%1.3e\n", l, q, (K+l)*(double)q*q);
! 133: #endif
! 134:
! 135: /* add 2 ulp to take into account rest of summation */
! 136: mpfr_add_one_ulp(s);
! 137: mpfr_add_one_ulp(s);
! 138:
! 139: for (k=0;k<K;k++) {
! 140: mpfr_mul(s, s, s, GMP_RNDU);
! 141: #ifdef DEBUG
! 142: printf("k=%d s=%1.20e\n",k,mpfr_get_d(s));
! 143: #endif
! 144: }
! 145:
! 146: if (n>0) mpfr_mul_2exp(s, s, n, GMP_RNDU);
! 147: else mpfr_div_2exp(s, s, -n, GMP_RNDU);
! 148:
! 149: /* error is at most 2^K*(2l+18) ulp */
! 150: l = 2*l+17; k=0; while (l) { k++; l >>= 1; }
! 151: /* now k = ceil(log(2l+18)/log(2)) */
! 152: K += k;
! 153: #ifdef DEBUG
! 154: printf("after mult. by 2^n:\n");
! 155: if (EXP(s)>-1024) printf("s=%1.20e\n",mpfr_get_d(s));
! 156: printf(" ="); mpfr_print_raw(s); putchar('\n');
! 157: printf("err=%d bits\n", K);
! 158: #endif
! 159:
! 160: l = mpfr_can_round(s, q-K, GMP_RNDU, rnd_mode, precy);
! 161: if (l==0) {
! 162: #ifdef DEBUG
! 163: printf("not enough precision, use %d\n", q+BITS_PER_MP_LIMB);
! 164: printf("q=%d q-K=%d precy=%d\n",q,q-K,precy);
! 165: #endif
! 166: q += BITS_PER_MP_LIMB;
! 167: mpfr_set_prec(r, q); mpfr_set_prec(s, q); mpfr_set_prec(t, q);
! 168: }
! 169: } while (l==0);
! 170:
! 171: mpfr_set(y, s, rnd_mode);
! 172:
! 173: mpfr_clear(r); mpfr_clear(s); mpfr_clear(t);
! 174: return 1;
! 175: }
! 176:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to exp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr