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