Annotation of OpenXM_contrib/gmp/mpfr/exp3.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_exp -- exponential of a floating-point number
! 2:
! 3: Copyright 1999, 2001, 2002 Free Software Foundation, Inc.
! 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 Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 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 Lesser General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Lesser 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 "gmp.h"
! 24: #include "gmp-impl.h"
! 25: #include "mpfr.h"
! 26: #include "mpfr-impl.h"
! 27:
! 28: static int mpfr_exp_rational _PROTO ((mpfr_ptr, mpz_srcptr, int, int));
! 29:
! 30: static int
! 31: mpfr_exp_rational (mpfr_ptr y, mpz_srcptr p, int r, int m)
! 32: {
! 33: int n,i,k,j,l;
! 34: mpz_t* P,*S;
! 35: mpz_t* ptoj;
! 36: int diff,expo;
! 37: int precy = MPFR_PREC(y);
! 38: int * mult;
! 39: int prec_i_have;
! 40: int *nb_terms;
! 41: int accu;
! 42: TMP_DECL (marker);
! 43:
! 44: TMP_MARK (marker);
! 45: n = 1 << m;
! 46: P = (mpz_t*) TMP_ALLOC((m+1) * sizeof(mpz_t));
! 47: S = (mpz_t*) TMP_ALLOC((m+1) * sizeof(mpz_t));
! 48: ptoj = (mpz_t*) TMP_ALLOC((m+1) * sizeof(mpz_t)); /* ptoj[i] = mantissa^(2^i) */
! 49: mult = (int*) TMP_ALLOC((m+1) * sizeof(int));
! 50: nb_terms = (int*) TMP_ALLOC((m+1) * sizeof(int));
! 51: mult[0] = 0;
! 52: for (i=0;i<=m;i++) { mpz_init(P[i]); mpz_init(S[i]); mpz_init(ptoj[i]); }
! 53: mpz_set(ptoj[0], p);
! 54: for (i=1;i<m;i++) mpz_mul(ptoj[i], ptoj[i-1], ptoj[i-1]);
! 55: mpz_set_ui(P[0], 1);
! 56: mpz_set_ui(S[0], 1);
! 57: k = 0;
! 58: nb_terms[0] = 1;
! 59: prec_i_have = 0;
! 60: for (i=1;(prec_i_have < precy) && (i < n) ;i++) {
! 61: k++;
! 62: nb_terms[k] = 1;
! 63: mpz_set_ui(P[k], i+1);
! 64: mpz_set(S[k], P[k]);;
! 65: j=i+1; l=0; while ((j & 1) == 0) {
! 66: mpz_mul(S[k], S[k], ptoj[l]);
! 67: mpz_mul(S[k-1], S[k-1], P[k]);
! 68: mpz_mul_2exp(S[k-1], S[k-1], r*(1<<l));
! 69: mpz_add(S[k-1], S[k-1], S[k]);
! 70: mpz_mul(P[k-1], P[k-1], P[k]);
! 71: nb_terms[k-1] = nb_terms[k-1]+ nb_terms[k];
! 72: mult[k] = mult[k-1] + (1 << l)*(r >> 2) + mpz_sizeinbase(P[k],2) - 1;
! 73: prec_i_have = mult[k];
! 74: l++; j>>=1; k--;
! 75: }
! 76: }
! 77: l = 0;
! 78: accu = 0;
! 79: while (k > 0){
! 80: mpz_mul(S[k], S[k], ptoj[_mpfr_ceil_log2((double) nb_terms[k])]);
! 81: mpz_mul(S[k-1], S[k-1], P[k]);
! 82: accu += nb_terms[k];
! 83: mpz_mul_2exp(S[k-1], S[k-1], r* accu);
! 84: mpz_add(S[k-1], S[k-1], S[k]);
! 85: mpz_mul(P[k-1], P[k-1], P[k]);
! 86: l++; k--;
! 87: }
! 88:
! 89: diff = mpz_sizeinbase(S[0],2) - 2*precy;
! 90: expo = diff;
! 91: if (diff >=0)
! 92: {
! 93: mpz_div_2exp(S[0],S[0],diff);
! 94: } else
! 95: {
! 96: mpz_mul_2exp(S[0],S[0],-diff);
! 97: }
! 98: diff = mpz_sizeinbase(P[0],2) - precy;
! 99: expo -= diff;
! 100: if (diff >=0)
! 101: {
! 102: mpz_div_2exp(P[0],P[0],diff);
! 103: } else
! 104: {
! 105: mpz_mul_2exp(P[0],P[0],-diff);
! 106: }
! 107:
! 108: mpz_tdiv_q(S[0], S[0], P[0]);
! 109: mpfr_set_z(y,S[0], GMP_RNDD);
! 110: MPFR_EXP(y) += expo;
! 111:
! 112: mpfr_div_2ui(y, y, r*(i-1),GMP_RNDN);
! 113: for (i=0;i<=m;i++) { mpz_clear(P[i]); mpz_clear(S[i]); mpz_clear(ptoj[i]); }
! 114: TMP_FREE (marker);
! 115: return 0;
! 116: }
! 117:
! 118: #define shift (BITS_PER_MP_LIMB/2)
! 119:
! 120: int
! 121: mpfr_exp3 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
! 122: {
! 123: mpfr_t t;
! 124: mpfr_t x_copy;
! 125: int i,k;
! 126: mpz_t uk;
! 127: mpfr_t tmp;
! 128: int ttt;
! 129: int twopoweri;
! 130: int Prec;
! 131: int loop;
! 132: int prec_x;
! 133: int shift_x = 0;
! 134: int good = 0;
! 135: int realprec = 0;
! 136: int iter;
! 137: int logn, inexact = 0;
! 138:
! 139: /* decompose x */
! 140: /* we first write x = 1.xxxxxxxxxxxxx
! 141: ----- k bits -- */
! 142: prec_x = _mpfr_ceil_log2 ((double) (MPFR_PREC(x)) / BITS_PER_MP_LIMB);
! 143: if (prec_x < 0) prec_x = 0;
! 144: logn = _mpfr_ceil_log2 ((double) prec_x + MPFR_PREC(y));
! 145: if (logn < 2) logn = 2;
! 146: ttt = MPFR_EXP(x);
! 147: mpfr_init2(x_copy,MPFR_PREC(x));
! 148: mpfr_set(x_copy,x,GMP_RNDD);
! 149: /* we shift to get a number less than 1 */
! 150: if (ttt > 0)
! 151: {
! 152: shift_x = ttt;
! 153: mpfr_div_2ui(x_copy, x, ttt, GMP_RNDN);
! 154: ttt = MPFR_EXP(x_copy);
! 155: }
! 156: realprec = MPFR_PREC(y)+logn;
! 157: mpz_init (uk);
! 158: while (!good){
! 159: Prec = realprec+shift+2+shift_x;
! 160: k = _mpfr_ceil_log2 ((double) Prec / BITS_PER_MP_LIMB);
! 161:
! 162: /* now we have to extract */
! 163: mpfr_init2 (t, Prec);
! 164: mpfr_init2 (tmp, Prec);
! 165: mpfr_set_ui(tmp,1,GMP_RNDN);
! 166: twopoweri = BITS_PER_MP_LIMB;
! 167: if (k <= prec_x) iter = k; else iter= prec_x;
! 168: for(i = 0; i <= iter; i++){
! 169: mpfr_extract (uk, x_copy, i);
! 170: if (i)
! 171: mpfr_exp_rational (t, uk, twopoweri - ttt, k - i + 1);
! 172: else
! 173: {
! 174: /* particular case: we have to compute with x/2^., then
! 175: do squarings (this is faster) */
! 176: mpfr_exp_rational (t, uk, shift + twopoweri - ttt, k+1);
! 177: for (loop= 0 ; loop < shift; loop++)
! 178: mpfr_mul(t,t,t,GMP_RNDD);
! 179:
! 180: }
! 181: mpfr_mul(tmp,tmp,t,GMP_RNDD);
! 182: twopoweri <<= 1;
! 183: }
! 184: mpfr_clear (t);
! 185: for (loop= 0 ; loop < shift_x; loop++)
! 186: mpfr_mul (tmp, tmp, tmp, GMP_RNDD);
! 187: if (mpfr_can_round (tmp, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(y)))
! 188: {
! 189: inexact = mpfr_set (y, tmp, rnd_mode);
! 190: good = 1;
! 191: }
! 192: else
! 193: realprec += 3*logn;
! 194: mpfr_clear (tmp);
! 195: }
! 196: mpz_clear (uk);
! 197: mpfr_clear(x_copy);
! 198: return inexact;
! 199: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to exp3.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr