Annotation of OpenXM_contrib/gmp/mpfr/generic.c, Revision 1.1.1.1
1.1 ohara 1: /* generic file for evaluation of hypergeometric series using binary splitting
2:
3: Copyright 1999, 2000, 2001 Free Software Foundation.
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 MPdFR 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: #ifndef GENERIC
23: # error You should specify a name
24: #endif
25:
26: #ifdef B
27: # ifndef A
28: # error B cannot be used without A
29: # endif
30: #endif
31:
32: /* Compute the first 2^m terms from the hypergeometric series
33: with x = p / 2^r */
34: static int
35: GENERIC (mpfr_ptr y, mpz_srcptr p, int r, int m)
36: {
37: int n,i,k,j,l;
38: int is_p_one = 0;
39: mpz_t* P,*S;
40: #ifdef A
41: mpz_t *T;
42: #endif
43: mpz_t* ptoj;
44: #ifdef R_IS_RATIONAL
45: mpz_t* qtoj;
46: mpfr_t tmp;
47: #endif
48: int diff, expo;
49: int precy = MPFR_PREC(y);
50: TMP_DECL(marker);
51:
52: TMP_MARK(marker);
53: MPFR_CLEAR_FLAGS(y);
54: n = 1 << m;
55: P = (mpz_t*) TMP_ALLOC ((m+1) * sizeof(mpz_t));
56: S = (mpz_t*) TMP_ALLOC ((m+1) * sizeof(mpz_t));
57: ptoj = (mpz_t*) TMP_ALLOC ((m+1) * sizeof(mpz_t)); /* ptoj[i] = mantissa^(2^i) */
58: #ifdef A
59: T = (mpz_t*) TMP_ALLOC ((m+1) * sizeof(mpz_t));
60: #endif
61: #ifdef R_IS_RATIONAL
62: qtoj = (mpz_t*) TMP_ALLOC ((m+1) * sizeof(mpz_t));
63: #endif
64: for (i=0;i<=m;i++)
65: {
66: mpz_init (P[i]);
67: mpz_init (S[i]);
68: mpz_init (ptoj[i]);
69: #ifdef R_IS_RATIONAL
70: mpz_init (qtoj[i]);
71: #endif
72: #ifdef A
73: mpz_init (T[i]);
74: #endif
75: }
76: mpz_set (ptoj[0], p);
77: #ifdef C
78: # if C2 != 1
79: mpz_mul_ui(ptoj[0], ptoj[0], C2);
80: # endif
81: #endif
82: is_p_one = !mpz_cmp_si(ptoj[0], 1);
83: #ifdef A
84: # ifdef B
85: mpz_set_ui(T[0], A1 * B1);
86: # else
87: mpz_set_ui(T[0], A1);
88: # endif
89: #endif
90: if (!is_p_one)
91: for (i=1;i<m;i++) mpz_mul(ptoj[i], ptoj[i-1], ptoj[i-1]);
92: #ifdef R_IS_RATIONAL
93: mpz_set_si(qtoj[0], r);
94: for (i=1;i<=m;i++)
95: {
96: mpz_mul(qtoj[i], qtoj[i-1], qtoj[i-1]);
97: }
98: #endif
99:
100: mpz_set_ui(P[0], 1);
101: mpz_set_ui(S[0], 1);
102: k = 0;
103: for (i=1;(i < n) ;i++) {
104: k++;
105:
106: #ifdef A
107: # ifdef B
108: mpz_set_ui(T[k], (A1 + A2*i)*(B1+B2*i));
109: # else
110: mpz_set_ui(T[k], A1 + A2*i);
111: # endif
112: #endif
113:
114: #ifdef C
115: # ifdef NO_FACTORIAL
116: mpz_set_ui(P[k], (C1 + C2 * (i-1)));
117: mpz_set_ui(S[k], 1);
118: # else
119: mpz_set_ui(P[k], (i+1) * (C1 + C2 * (i-1)));
120: mpz_set_ui(S[k], i+1);
121: # endif
122: #else
123: # ifdef NO_FACTORIAL
124: mpz_set_ui(P[k], 1);
125: # else
126: mpz_set_ui(P[k], i+1);
127: # endif
128: mpz_set(S[k], P[k]);
129: #endif
130: j=i+1; l=0; while ((j & 1) == 0) {
131: if (!is_p_one)
132: mpz_mul(S[k], S[k], ptoj[l]);
133: #ifdef A
134: # ifdef B
135: # if (A2*B2) != 1
136: mpz_mul_ui(P[k], P[k], A2*B2);
137: # endif
138: # else
139: # if A2 != 1
140: mpz_mul_ui(P[k], P[k], A2);
141: # endif
142: #endif
143: mpz_mul(S[k], S[k], T[k-1]);
144: #endif
145: mpz_mul(S[k-1], S[k-1], P[k]);
146: #ifdef R_IS_RATIONAL
147: mpz_mul(S[k-1], S[k-1], qtoj[l]);
148: #else
149: mpz_mul_2exp(S[k-1], S[k-1], r*(1<<l));
150: #endif
151: mpz_add(S[k-1], S[k-1], S[k]);
152: mpz_mul(P[k-1], P[k-1], P[k]);
153: #ifdef A
154: mpz_mul(T[k-1], T[k-1], T[k]);
155: #endif
156: l++; j>>=1; k--;
157: }
158: }
159:
160: diff = mpz_sizeinbase(S[0],2) - 2*precy;
161: expo = diff;
162: if (diff >=0)
163: {
164: mpz_div_2exp(S[0],S[0],diff);
165: } else
166: {
167: mpz_mul_2exp(S[0],S[0],-diff);
168: }
169: diff = mpz_sizeinbase(P[0],2) - precy;
170: expo -= diff;
171: if (diff >=0)
172: {
173: mpz_div_2exp(P[0],P[0],diff);
174: } else
175: {
176: mpz_mul_2exp(P[0],P[0],-diff);
177: }
178:
179: mpz_tdiv_q(S[0], S[0], P[0]);
180: mpfr_set_z(y, S[0], GMP_RNDD);
181: MPFR_EXP(y) += expo;
182:
183: #ifdef R_IS_RATIONAL
184: /* exact division */
185: mpz_div_ui (qtoj[m], qtoj[m], r);
186: mpfr_init2 (tmp, MPFR_PREC(y));
187: mpfr_set_z (tmp, qtoj[m] , GMP_RNDD);
188: mpfr_div (y, y, tmp, GMP_RNDD);
189: mpfr_clear (tmp);
190: #else
191: mpfr_div_2ui(y, y, r*(i-1), GMP_RNDN);
192: #endif
193: for (i=0;i<=m;i++)
194: {
195: mpz_clear (P[i]);
196: mpz_clear (S[i]);
197: mpz_clear (ptoj[i]);
198: #ifdef R_IS_RATIONAL
199: mpz_clear (qtoj[i]);
200: #endif
201: #ifdef A
202: mpz_clear (T[i]);
203: #endif
204: }
205: TMP_FREE(marker);
206: return 0;
207: }
208:
209:
210:
211:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to generic.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr