Annotation of OpenXM_contrib/gmp/mpfr/exp_2.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_exp_2 -- exponential of a floating-point number
! 2: using Brent's algorithms in O(n^(1/2)*M(n)) and O(n^(1/3)*M(n))
! 3:
! 4: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
! 5:
! 6: This file is part of the MPFR Library.
! 7:
! 8: The MPFR Library is free software; you can redistribute it and/or modify
! 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
! 11: option) any later version.
! 12:
! 13: The MPFR Library is distributed in the hope that it will be useful, but
! 14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 16: License for more details.
! 17:
! 18: You should have received a copy of the GNU Lesser General Public License
! 19: along with the MPFR Library; see the file COPYING.LIB. If not, write to
! 20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 21: MA 02111-1307, USA. */
! 22:
! 23: #include <stdio.h>
! 24: #include "gmp.h"
! 25: #include "gmp-impl.h"
! 26: #include "mpfr.h"
! 27: #include "mpfr-impl.h"
! 28:
! 29: static int mpfr_exp2_aux _PROTO ((mpz_t, mpfr_srcptr, int, int*));
! 30: static int mpfr_exp2_aux2 _PROTO ((mpz_t, mpfr_srcptr, int, int*));
! 31: static mp_exp_t mpz_normalize _PROTO ((mpz_t, mpz_t, int));
! 32: static int mpz_normalize2 _PROTO ((mpz_t, mpz_t, int, int));
! 33:
! 34: /* returns floor(sqrt(n)) */
! 35: unsigned long
! 36: _mpfr_isqrt (unsigned long n)
! 37: {
! 38: unsigned long s;
! 39:
! 40: s = 1;
! 41: do {
! 42: s = (s + n / s) / 2;
! 43: } while (!(s*s <= n && n <= s*(s+2)));
! 44: return s;
! 45: }
! 46:
! 47: /* returns floor(n^(1/3)) */
! 48: unsigned long
! 49: _mpfr_cuberoot (unsigned long n)
! 50: {
! 51: double s, is;
! 52:
! 53: s = 1.0;
! 54: do {
! 55: s = (2*s*s*s + (double) n) / (3*s*s);
! 56: is = (double) ((int) s);
! 57: } while (!(is*is*is <= (double) n && (double) n < (is+1)*(is+1)*(is+1)));
! 58: return (unsigned long) is;
! 59: }
! 60:
! 61: #define SWITCH 100 /* number of bits to switch from O(n^(1/2)*M(n)) method
! 62: to O(n^(1/3)*M(n)) method */
! 63:
! 64: #define MY_INIT_MPZ(x, s) { \
! 65: (x)->_mp_alloc = (s); \
! 66: PTR(x) = (mp_ptr) TMP_ALLOC((s)*BYTES_PER_MP_LIMB); \
! 67: (x)->_mp_size = 0; }
! 68:
! 69: /* #define DEBUG */
! 70:
! 71: /* if k = the number of bits of z > q, divides z by 2^(k-q) and returns k-q.
! 72: Otherwise do nothing and return 0.
! 73: */
! 74: static mp_exp_t
! 75: mpz_normalize (mpz_t rop, mpz_t z, int q)
! 76: {
! 77: int k;
! 78:
! 79: k = mpz_sizeinbase(z, 2);
! 80: if (k > q) {
! 81: mpz_div_2exp(rop, z, k-q);
! 82: return k-q;
! 83: }
! 84: else {
! 85: if (rop != z) mpz_set(rop, z);
! 86: return 0;
! 87: }
! 88: }
! 89:
! 90: /* if expz > target, shift z by (expz-target) bits to the left.
! 91: if expz < target, shift z by (target-expz) bits to the right.
! 92: Returns target.
! 93: */
! 94: static int
! 95: mpz_normalize2 (mpz_t rop, mpz_t z, int expz, int target)
! 96: {
! 97: if (target > expz) mpz_div_2exp(rop, z, target-expz);
! 98: else mpz_mul_2exp(rop, z, expz-target);
! 99: return target;
! 100: }
! 101:
! 102: /* use Brent's formula exp(x) = (1+r+r^2/2!+r^3/3!+...)^(2^K)*2^n
! 103: where x = n*log(2)+(2^K)*r
! 104: together with Brent-Kung O(t^(1/2)) algorithm for the evaluation of
! 105: power series. The resulting complexity is O(n^(1/3)*M(n)).
! 106: */
! 107: int
! 108: mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
! 109: {
! 110: int n, K, precy, q, k, l, err, exps, inexact;
! 111: mpfr_t r, s, t; mpz_t ss;
! 112: TMP_DECL(marker);
! 113:
! 114: precy = MPFR_PREC(y);
! 115:
! 116: n = (int) (mpfr_get_d1 (x) / LOG2);
! 117:
! 118: /* for the O(n^(1/2)*M(n)) method, the Taylor series computation of
! 119: n/K terms costs about n/(2K) multiplications when computed in fixed
! 120: point */
! 121: K = (precy<SWITCH) ? _mpfr_isqrt((precy + 1) / 2) : _mpfr_cuberoot (4*precy);
! 122: l = (precy-1)/K + 1;
! 123: err = K + (int) _mpfr_ceil_log2 (2.0 * (double) l + 18.0);
! 124: /* add K extra bits, i.e. failure probability <= 1/2^K = O(1/precy) */
! 125: q = precy + err + K + 3;
! 126: mpfr_init2 (r, q);
! 127: mpfr_init2 (s, q);
! 128: mpfr_init2 (t, q);
! 129: /* the algorithm consists in computing an upper bound of exp(x) using
! 130: a precision of q bits, and see if we can round to MPFR_PREC(y) taking
! 131: into account the maximal error. Otherwise we increase q. */
! 132: do {
! 133: #ifdef DEBUG
! 134: printf("n=%d K=%d l=%d q=%d\n",n,K,l,q);
! 135: #endif
! 136:
! 137: /* if n<0, we have to get an upper bound of log(2)
! 138: in order to get an upper bound of r = x-n*log(2) */
! 139: mpfr_const_log2 (s, (n>=0) ? GMP_RNDZ : GMP_RNDU);
! 140: #ifdef DEBUG
! 141: printf("n=%d log(2)=",n); mpfr_print_binary(s); putchar('\n');
! 142: #endif
! 143: mpfr_mul_ui (r, s, (n<0) ? -n : n, (n>=0) ? GMP_RNDZ : GMP_RNDU);
! 144: if (n<0) mpfr_neg(r, r, GMP_RNDD);
! 145: /* r = floor(n*log(2)) */
! 146:
! 147: #ifdef DEBUG
! 148: printf("x=%1.20e\n", mpfr_get_d1 (x));
! 149: printf(" ="); mpfr_print_binary(x); putchar('\n');
! 150: printf("r=%1.20e\n", mpfr_get_d1 (r));
! 151: printf(" ="); mpfr_print_binary(r); putchar('\n');
! 152: #endif
! 153: mpfr_sub(r, x, r, GMP_RNDU);
! 154: if (MPFR_SIGN(r)<0) { /* initial approximation n was too large */
! 155: n--;
! 156: mpfr_mul_ui(r, s, (n<0) ? -n : n, GMP_RNDZ);
! 157: if (n<0) mpfr_neg(r, r, GMP_RNDD);
! 158: mpfr_sub(r, x, r, GMP_RNDU);
! 159: }
! 160: #ifdef DEBUG
! 161: printf("x-r=%1.20e\n", mpfr_get_d1 (r));
! 162: printf(" ="); mpfr_print_binary(r); putchar('\n');
! 163: if (MPFR_SIGN(r)<0) { fprintf(stderr,"Error in mpfr_exp: r<0\n"); exit(1); }
! 164: #endif
! 165: mpfr_div_2ui(r, r, K, GMP_RNDU); /* r = (x-n*log(2))/2^K */
! 166:
! 167: TMP_MARK(marker);
! 168: MY_INIT_MPZ(ss, 3 + 2*((q-1)/BITS_PER_MP_LIMB));
! 169: exps = mpfr_get_z_exp(ss, s);
! 170: /* s <- 1 + r/1! + r^2/2! + ... + r^l/l! */
! 171: l = (precy<SWITCH) ? mpfr_exp2_aux(ss, r, q, &exps) /* naive method */
! 172: : mpfr_exp2_aux2(ss, r, q, &exps); /* Brent/Kung method */
! 173:
! 174: #ifdef DEBUG
! 175: printf("l=%d q=%d (K+l)*q^2=%1.3e\n", l, q, (K+l)*(double)q*q);
! 176: #endif
! 177:
! 178: for (k=0;k<K;k++) {
! 179: mpz_mul(ss, ss, ss); exps <<= 1;
! 180: exps += mpz_normalize(ss, ss, q);
! 181: }
! 182: mpfr_set_z(s, ss, GMP_RNDN); MPFR_EXP(s) += exps;
! 183: TMP_FREE(marker); /* don't need ss anymore */
! 184:
! 185: if (n>0) mpfr_mul_2ui(s, s, n, GMP_RNDU);
! 186: else mpfr_div_2ui(s, s, -n, GMP_RNDU);
! 187:
! 188: /* error is at most 2^K*(3l*(l+1)) ulp for mpfr_exp2_aux */
! 189: if (precy<SWITCH) l = 3*l*(l+1);
! 190: else l = l*(l+4);
! 191: k=0; while (l) { k++; l >>= 1; }
! 192: /* now k = ceil(log(error in ulps)/log(2)) */
! 193: K += k;
! 194: #ifdef DEBUG
! 195: printf("after mult. by 2^n:\n");
! 196: if (MPFR_EXP(s) > -1024)
! 197: printf("s=%1.20e\n", mpfr_get_d1 (s));
! 198: printf(" ="); mpfr_print_binary(s); putchar('\n');
! 199: printf("err=%d bits\n", K);
! 200: #endif
! 201:
! 202: l = mpfr_can_round(s, q-K, GMP_RNDN, rnd_mode, precy);
! 203: if (l==0) {
! 204: #ifdef DEBUG
! 205: printf("not enough precision, use %d\n", q+BITS_PER_MP_LIMB);
! 206: printf("q=%d q-K=%d precy=%d\n",q,q-K,precy);
! 207: #endif
! 208: q += BITS_PER_MP_LIMB;
! 209: mpfr_set_prec(r, q); mpfr_set_prec(s, q); mpfr_set_prec(t, q);
! 210: }
! 211: } while (l==0);
! 212:
! 213: inexact = mpfr_set (y, s, rnd_mode);
! 214:
! 215: mpfr_clear(r); mpfr_clear(s); mpfr_clear(t);
! 216: return inexact;
! 217: }
! 218:
! 219: /* s <- 1 + r/1! + r^2/2! + ... + r^l/l! while MPFR_EXP(r^l/l!)+MPFR_EXPR(r)>-q
! 220: using naive method with O(l) multiplications.
! 221: Return the number of iterations l.
! 222: The absolute error on s is less than 3*l*(l+1)*2^(-q).
! 223: Version using fixed-point arithmetic with mpz instead
! 224: of mpfr for internal computations.
! 225: s must have at least qn+1 limbs (qn should be enough, but currently fails
! 226: since mpz_mul_2exp(s, s, q-1) reallocates qn+1 limbs)
! 227: */
! 228: static int
! 229: mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, int q, int *exps)
! 230: {
! 231: int l, dif, qn;
! 232: mpz_t t, rr; mp_exp_t expt, expr;
! 233: TMP_DECL(marker);
! 234:
! 235: TMP_MARK(marker);
! 236: qn = 1 + (q-1)/BITS_PER_MP_LIMB;
! 237: MY_INIT_MPZ(t, 2*qn+1); /* 2*qn+1 is neeeded since mpz_div_2exp may
! 238: need an extra limb */
! 239: MY_INIT_MPZ(rr, qn+1);
! 240: mpz_set_ui(t, 1); expt=0;
! 241: mpz_set_ui(s, 1); mpz_mul_2exp(s, s, q-1); *exps = 1-q; /* s = 2^(q-1) */
! 242: expr = mpfr_get_z_exp(rr, r); /* no error here */
! 243:
! 244: l = 0;
! 245: do {
! 246: l++;
! 247: mpz_mul(t, t, rr); expt=expt+expr;
! 248: dif = *exps + mpz_sizeinbase(s, 2) - expt - mpz_sizeinbase(t, 2);
! 249: /* truncates the bits of t which are < ulp(s) = 2^(1-q) */
! 250: expt += mpz_normalize(t, t, q-dif); /* error at most 2^(1-q) */
! 251: mpz_div_ui(t, t, l); /* error at most 2^(1-q) */
! 252: /* the error wrt t^l/l! is here at most 3*l*ulp(s) */
! 253: #ifdef DEBUG
! 254: if (expt != *exps) {
! 255: fprintf(stderr, "Error: expt != exps %d %d\n", expt, *exps); exit(1);
! 256: }
! 257: #endif
! 258: mpz_add(s, s, t); /* no error here: exact */
! 259: /* ensures rr has the same size as t: after several shifts, the error
! 260: on rr is still at most ulp(t)=ulp(s) */
! 261: expr += mpz_normalize(rr, rr, mpz_sizeinbase(t, 2));
! 262: } while (mpz_cmp_ui(t, 0));
! 263:
! 264: TMP_FREE(marker);
! 265: return l;
! 266: }
! 267:
! 268: /* s <- 1 + r/1! + r^2/2! + ... + r^l/l! while MPFR_EXP(r^l/l!)+MPFR_EXPR(r)>-q
! 269: using Brent/Kung method with O(sqrt(l)) multiplications.
! 270: Return l.
! 271: Uses m multiplications of full size and 2l/m of decreasing size,
! 272: i.e. a total equivalent to about m+l/m full multiplications,
! 273: i.e. 2*sqrt(l) for m=sqrt(l).
! 274: Version using mpz. ss must have at least (sizer+1) limbs.
! 275: The error is bounded by (l^2+4*l) ulps where l is the return value.
! 276: */
! 277: static int
! 278: mpfr_exp2_aux2 (mpz_t s, mpfr_srcptr r, int q, int *exps)
! 279: {
! 280: int expr, l, m, i, sizer, *expR, expt, ql;
! 281: unsigned long int c;
! 282: mpz_t t, *R, rr, tmp;
! 283: TMP_DECL(marker);
! 284:
! 285: /* estimate value of l */
! 286: l = q / (-MPFR_EXP(r));
! 287: m = (int) _mpfr_isqrt (l);
! 288: /* we access R[2], thus we need m >= 2 */
! 289: if (m < 2) m = 2;
! 290: TMP_MARK(marker);
! 291: R = (mpz_t*) TMP_ALLOC((m+1)*sizeof(mpz_t)); /* R[i] stands for r^i */
! 292: expR = (int*) TMP_ALLOC((m+1)*sizeof(int)); /* exponent for R[i] */
! 293: sizer = 1 + (MPFR_PREC(r)-1)/BITS_PER_MP_LIMB;
! 294: mpz_init(tmp);
! 295: MY_INIT_MPZ(rr, sizer+2);
! 296: MY_INIT_MPZ(t, 2*sizer); /* double size for products */
! 297: mpz_set_ui(s, 0); *exps = 1-q; /* initialize s to zero, 1 ulp = 2^(1-q) */
! 298: for (i=0;i<=m;i++) MY_INIT_MPZ(R[i], sizer+2);
! 299: expR[1] = mpfr_get_z_exp(R[1], r); /* exact operation: no error */
! 300: expR[1] = mpz_normalize2(R[1], R[1], expR[1], 1-q); /* error <= 1 ulp */
! 301: mpz_mul(t, R[1], R[1]); /* err(t) <= 2 ulps */
! 302: mpz_div_2exp(R[2], t, q-1); /* err(R[2]) <= 3 ulps */
! 303: expR[2] = 1-q;
! 304: for (i=3;i<=m;i++) {
! 305: mpz_mul(t, R[i-1], R[1]); /* err(t) <= 2*i-2 */
! 306: mpz_div_2exp(R[i], t, q-1); /* err(R[i]) <= 2*i-1 ulps */
! 307: expR[i] = 1-q;
! 308: }
! 309: mpz_set_ui(R[0], 1); mpz_mul_2exp(R[0], R[0], q-1); expR[0]=1-q; /* R[0]=1 */
! 310: mpz_set_ui(rr, 1); expr=0; /* rr contains r^l/l! */
! 311: /* by induction: err(rr) <= 2*l ulps */
! 312:
! 313: l = 0;
! 314: ql = q; /* precision used for current giant step */
! 315: do {
! 316: /* all R[i] must have exponent 1-ql */
! 317: if (l) for (i=0;i<m;i++) {
! 318: expR[i] = mpz_normalize2(R[i], R[i], expR[i], 1-ql);
! 319: }
! 320: /* the absolute error on R[i]*rr is still 2*i-1 ulps */
! 321: expt = mpz_normalize2(t, R[m-1], expR[m-1], 1-ql);
! 322: /* err(t) <= 2*m-1 ulps */
! 323: /* computes t = 1 + r/(l+1) + ... + r^(m-1)*l!/(l+m-1)!
! 324: using Horner's scheme */
! 325: for (i=m-2;i>=0;i--) {
! 326: mpz_div_ui(t, t, l+i+1); /* err(t) += 1 ulp */
! 327: mpz_add(t, t, R[i]);
! 328: }
! 329: /* now err(t) <= (3m-2) ulps */
! 330:
! 331: /* now multiplies t by r^l/l! and adds to s */
! 332: mpz_mul(t, t, rr); expt += expr;
! 333: expt = mpz_normalize2(t, t, expt, *exps);
! 334: /* err(t) <= (3m-1) + err_rr(l) <= (3m-2) + 2*l */
! 335: #ifdef DEBUG
! 336: if (expt != *exps) {
! 337: fprintf(stderr, "Error: expt != exps %d %d\n", expt, *exps); exit(1);
! 338: }
! 339: #endif
! 340: mpz_add(s, s, t); /* no error here */
! 341:
! 342: /* updates rr, the multiplication of the factors l+i could be done
! 343: using binary splitting too, but it is not sure it would save much */
! 344: mpz_mul(t, rr, R[m]); /* err(t) <= err(rr) + 2m-1 */
! 345: expr += expR[m];
! 346: mpz_set_ui (tmp, 1);
! 347: for (i=1, c=1; i<=m; i++) {
! 348: if (l+i > ~((unsigned long int) 0)/c) {
! 349: mpz_mul_ui(tmp, tmp, c);
! 350: c = l+i;
! 351: }
! 352: else c *= (unsigned long int) l+i;
! 353: }
! 354: if (c != 1) mpz_mul_ui (tmp, tmp, c); /* tmp is exact */
! 355: mpz_fdiv_q(t, t, tmp); /* err(t) <= err(rr) + 2m */
! 356: expr += mpz_normalize(rr, t, ql); /* err_rr(l+1) <= err_rr(l) + 2m+1 */
! 357: ql = q - *exps - mpz_sizeinbase(s, 2) + expr + mpz_sizeinbase(rr, 2);
! 358: l+=m;
! 359: } while (expr+mpz_sizeinbase(rr, 2) > -q);
! 360:
! 361: TMP_FREE(marker);
! 362: mpz_clear(tmp);
! 363: return l;
! 364: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to exp_2.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr