Annotation of OpenXM_contrib/gmp/mpfr/exp3.c, Revision 1.1.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