diff options
| -rw-r--r-- | exp.c | 107 | ||||
| -rw-r--r-- | exp2.c | 1 | ||||
| -rw-r--r-- | exp3.c | 271 | ||||
| -rw-r--r-- | extract.c | 60 |
4 files changed, 333 insertions, 106 deletions
@@ -51,9 +51,6 @@ mpfr_exp(y, x, rnd_mode) expx = EXP(x); precy = PREC(y); -#ifdef DEBUG - printf("EXP(x)=%d\n",expx); -#endif /* if x > (2^31-1)*ln(2), then exp(x) > 2^(2^31-1) i.e. gives +infinity */ if (expx > 30) { @@ -69,108 +66,8 @@ mpfr_exp(y, x, rnd_mode) mpfr_sub_one_ulp(y); return 1; } - n = (int) floor(mpfr_get_d(x)/LOG2); - - K = (int) sqrt( (double) precy ); - l = (precy-1)/K + 1; - err = K + (int) ceil(log(2.0*(double)l+18.0)/LOG2); - /* add K extra bits, i.e. failure probability <= 1/2^K = O(1/precy) */ - q = precy + err + K + 3; - mpfr_init2(r, q); mpfr_init2(s, q); mpfr_init2(t, q); - /* the algorithm consists in computing an upper bound of exp(x) using - a precision of q bits, and see if we can round to PREC(y) taking - into account the maximal error. Otherwise we increase q. */ - do { -#ifdef DEBUG - printf("n=%d K=%d l=%d q=%d\n",n,K,l,q); -#endif - - /* if n<0, we have to get an upper bound of log(2) - in order to get an upper bound of r = x-n*log(2) */ - mpfr_const_log2(s, (n>=0) ? GMP_RNDZ : GMP_RNDU); -#ifdef DEBUG - printf("n=%d log(2)=",n); mpfr_print_raw(s); putchar('\n'); -#endif - mpfr_mul_ui(r, s, (n<0) ? -n : n, (n>=0) ? GMP_RNDZ : GMP_RNDU); - if (n<0) mpfr_neg(r, r, GMP_RNDD); - /* r = floor(n*log(2)) */ - -#ifdef DEBUG - printf("x=%1.20e\n",mpfr_get_d(x)); - printf(" ="); mpfr_print_raw(x); putchar('\n'); - printf("r=%1.20e\n",mpfr_get_d(r)); - printf(" ="); mpfr_print_raw(r); putchar('\n'); -#endif - mpfr_sub(r, x, r, GMP_RNDU); - if (MPFR_SIGN(r)<0) { /* initial approximation n was too large */ - n--; - mpfr_mul_ui(r, s, (n<0) ? -n : n, GMP_RNDZ); - if (n<0) mpfr_neg(r, r, GMP_RNDD); - mpfr_sub(r, x, r, GMP_RNDU); - } -#ifdef DEBUG - printf("x-r=%1.20e\n",mpfr_get_d(r)); - printf(" ="); mpfr_print_raw(r); putchar('\n'); - if (MPFR_SIGN(r)<0) { fprintf(stderr,"Error in mpfr_exp: r<0\n"); exit(1); } -#endif - mpfr_div_2exp(r, r, K, GMP_RNDU); /* r = (x-n*log(2))/2^K */ - mpfr_set_ui(s, 1, GMP_RNDU); - mpfr_set_ui(t, 1, GMP_RNDU); - - l = 1; expr = EXP(r); - do { - mpfr_mul(t, t, r, GMP_RNDU); - mpfr_div_ui(t, t, l, GMP_RNDU); - mpfr_add(s, s, t, GMP_RNDU); -#ifdef DEBUG - printf("l=%d t=%1.20e\n",l,mpfr_get_d(t)); - printf("s=%1.20e\n",mpfr_get_d(s)); -#endif - l++; - } while (EXP(t)+expr > -q); -#ifdef DEBUG - printf("l=%d q=%d (K+l)*q^2=%1.3e\n", l, q, (K+l)*(double)q*q); -#endif - - /* add 2 ulp to take into account rest of summation */ - mpfr_add_one_ulp(s); - mpfr_add_one_ulp(s); - - for (k=0;k<K;k++) { - mpfr_mul(s, s, s, GMP_RNDU); -#ifdef DEBUG - printf("k=%d s=%1.20e\n",k,mpfr_get_d(s)); -#endif - } - - if (n>0) mpfr_mul_2exp(s, s, n, GMP_RNDU); - else mpfr_div_2exp(s, s, -n, GMP_RNDU); - - /* error is at most 2^K*(2l+18) ulp */ - l = 2*l+17; k=0; while (l) { k++; l >>= 1; } - /* now k = ceil(log(2l+18)/log(2)) */ - K += k; -#ifdef DEBUG - printf("after mult. by 2^n:\n"); - if (EXP(s)>-1024) printf("s=%1.20e\n",mpfr_get_d(s)); - printf(" ="); mpfr_print_raw(s); putchar('\n'); - printf("err=%d bits\n", K); -#endif - - l = mpfr_can_round(s, q-K, GMP_RNDU, rnd_mode, precy); - if (l==0) { -#ifdef DEBUG - printf("not enough precision, use %d\n", q+BITS_PER_MP_LIMB); - printf("q=%d q-K=%d precy=%d\n",q,q-K,precy); -#endif - q += BITS_PER_MP_LIMB; - mpfr_set_prec(r, q); mpfr_set_prec(s, q); mpfr_set_prec(t, q); - } - } while (l==0); - - mpfr_set(y, s, rnd_mode); - - mpfr_clear(r); mpfr_clear(s); mpfr_clear(t); + if (precy > 7000) mpfr_exp3(y,x,rnd_mode); + else mpfr_exp2(y,x,rnd_mode); return 1; } @@ -155,7 +155,6 @@ mpfr_exp2(y, x, rnd_mode) l = mpfr_can_round(s, q-K, GMP_RNDU, rnd_mode, precy); if (l==0) { -printf("not enough precision\n"); #ifdef DEBUG printf("not enough precision, use %d\n", q+BITS_PER_MP_LIMB); printf("q=%d q-K=%d precy=%d\n",q,q-K,precy); @@ -0,0 +1,271 @@ +/* mpfr_exp -- exponential of a floating-point number + +Copyright (C) 1999 PolKA project, Inria Lorraine and Loria + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <math.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "mpfr.h" + +/*#define DEBUG */ +/*#define TIMING */ + +int mpfr_extract(mpz_ptr , mpfr_srcptr , int ); + + +#define MY_INIT_MPZ(x, s) { \ + (x)->_mp_alloc = (s); \ + (x)->_mp_d = (mp_ptr) TMP_ALLOC((s)*BYTES_PER_MP_LIMB); \ + (x)->_mp_size = 0; } + +int +#if __STDC__ +mylog2(int x) +#else +mylog2(x) +int x; +#endif +{ + int i = 0; + for ( ; x != 1; x >>= 1, i++) ; + return i; +} + +int +#if __STDC__ +mpfr_exp_rational(mpfr_ptr y,mpz_srcptr p,int r,int m) +#else +mpfr_exp_rational(y,p,r,m) +mpfr_ptr y; +mpz_srcptr p; +int r; +int m; +#endif +{ + int n,i,k,j,l; + mpz_t* P,*S; + mpz_t* ptoj; + int diff,expo; + int precy = PREC(y); + int * mult; + int prec_i_have; + int *nb_terms; + int accu; + n = 1 << m; + P = (mpz_t*) malloc((m+1) * sizeof(mpz_t)); + S = (mpz_t*) malloc((m+1) * sizeof(mpz_t)); + + ptoj = (mpz_t*) malloc((m+1) * sizeof(mpz_t)); /* ptoj[i] = mantissa^(2^i) */ + mult = (int*) malloc((m+1) * sizeof(int)); + nb_terms = (int*) malloc((m+1) * sizeof(int)); + mult[0] = 0; + for (i=0;i<=m;i++) { mpz_init(P[i]); mpz_init(S[i]); mpz_init(ptoj[i]); + } + mpz_set(ptoj[0], p); + for (i=1;i<m;i++) mpz_mul(ptoj[i], ptoj[i-1], ptoj[i-1]); + mpz_set_ui(P[0], 1); + mpz_set_ui(S[0], 1); + k = 0; + prec_i_have = 0; + for (i=1;(prec_i_have < precy) && (i < n) ;i++) { + k++; + nb_terms[k] = 1; + mpz_set_ui(P[k], i+1); + /* mpz_set_ui(S[k], i+1);*/ + mpz_set(S[k], P[k]);; + j=i+1; l=0; while ((j & 1) == 0) { + mpz_mul(S[k], S[k], ptoj[l]); + mpz_mul(S[k-1], S[k-1], P[k]); + mpz_mul_2exp(S[k-1], S[k-1], r*(1<<l)); + mpz_add(S[k-1], S[k-1], S[k]); + mpz_mul(P[k-1], P[k-1], P[k]); + nb_terms[k-1] = nb_terms[k-1]+ nb_terms[k]; + mult[k] = mult[k-1] + (1 << l)*(r >> 2) + mpz_sizeinbase(P[k],2) - 1; + prec_i_have = mult[k]; + l++; j>>=1; k--; + } + } + l = 0; + accu = 0; + while (k > 0){ + mpz_mul(S[k], S[k], ptoj[mylog2(nb_terms[k])]); + mpz_mul(S[k-1], S[k-1], P[k]); + accu += nb_terms[k]; + mpz_mul_2exp(S[k-1], S[k-1], r* accu); + mpz_add(S[k-1], S[k-1], S[k]); + mpz_mul(P[k-1], P[k-1], P[k]); + l++; k--; + } + + diff = mpz_sizeinbase(S[0],2) - 2*precy; + expo = diff; + if (diff >=0) + { + mpz_div_2exp(S[0],S[0],diff); + } else + { + mpz_mul_2exp(S[0],S[0],-diff); + } + diff = mpz_sizeinbase(P[0],2) - precy; + expo -= diff; + if (diff >=0) + { + mpz_div_2exp(P[0],P[0],diff); + } else + { + mpz_mul_2exp(P[0],P[0],-diff); + } + mpz_tdiv_q(S[0], S[0], P[0]); + mpfr_set_z(y,S[0], GMP_RNDD); + EXP(y) += expo; + + + mpfr_div_2exp(y, y, r*(i-1),GMP_RNDN); + for (i=0;i<=m;i++) { mpz_clear(P[i]); mpz_clear(S[i]); mpz_clear(ptoj[i]); } + TMP_FREE(marker); + return 0; +} + +#define shift (BITS_PER_MP_LIMB/2) + +int +#if __STDC__ +mpfr_exp3(mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) +#else +mpfr_exp3(y,x,rnd_mode) +mpfr_ptr y; +mpfr_srcptr x; +mp_rnd_t rnd_mode; +#endif +{ + mpfr_t t; + mpfr_t x_copy; + int i,k; + mpz_t uk; + mpfr_t tmp; + int ttt; + int twopoweri; + int Prec; +#ifdef TIMING + int st; +#endif +#ifdef TIMING + int dummy; +#endif + int loop; + int prec_x; + int shift_x = 0; + int good = 0; + int realprec = 0; + int iter; + int logn; + /* commencons par 0 */ + if (FLAG_NAN(x)) { SET_NAN(y); return 1; } + if (!NOTZERO(x)) { mpfr_set_ui(y, 1, GMP_RNDN); return 0; } + + /* Decomposer x */ + /* on commence par ecrire x = 1.xxxxxxxxxxxxx + ----- k bits -- */ + prec_x = (int) ceil(log + ((double) (PREC(x)) / (double) BITS_PER_MP_LIMB) + /log(2.0)); + logn = (int) ceil(log + ((double) prec_x) + /log(2.0)); + if (logn < 2) logn = 2; + ttt = EXP(x); + mpfr_init2(x_copy,PREC(x)); + mpfr_set(x_copy,x,GMP_RNDD); + + /* on fait le shift pour que le nombre soit inferieur a 1 */ + if (ttt > 0) + { + shift_x = ttt; + mpfr_mul_2exp(x_copy,x,-ttt, GMP_RNDN); + ttt = EXP(x_copy); + } + realprec = PREC(y)+logn; + while (!good){ +/* Prec = realprec+shift +shift_x + 2*prec_x;*/ + Prec = realprec+shift+2+shift_x; + k = (int) ceil(log + ((double) (Prec) / (double) BITS_PER_MP_LIMB) + /log(2.0)); + /* Maintenant, il faut extraire : */ + mpfr_init2(t, Prec); + mpfr_init2(tmp, Prec); + mpfr_set_ui(tmp,1,GMP_RNDN); + twopoweri = BITS_PER_MP_LIMB; + if (k <= prec_x) iter = k; else iter= prec_x; + for(i = 0; i <= iter; i++){ + mpfr_extract(uk,x_copy,i); +#ifdef DEBUG + mpz_out_str(stderr,2, uk); + fprintf(stderr, "---\n"); + fprintf(stderr, "---%d\n", twopoweri - ttt); +#endif +#ifdef TIMING + st = cputime(); +#endif + if (i) +#ifdef TIMING + for (dummy = 0; dummy < 30; dummy++) +#endif + mpfr_exp_rational(t,uk,twopoweri - ttt, k - i + 1); + else + { + /* cas particulier : on est oblige de faire les calculs avec x/2^. + puis elever au carre (plus rapide) */ +#ifdef TIMING + for (dummy = 0; dummy < 30; dummy++) +#endif + mpfr_exp_rational(t,uk, shift + twopoweri - ttt, k+1); + for (loop= 0 ; loop < shift; loop++) + mpfr_mul(t,t,t,GMP_RNDD); + } + mpfr_mul(tmp,tmp,t,GMP_RNDD); +#ifdef TIMING + fprintf(stderr, "temps : %d ms \n", cputime() - st); +#endif +#ifdef DEBUG + fprintf(stderr, "fin\n"); + mpfr_out_str(stderr, 2, PREC(y), t, GMP_RNDD); + fprintf(stderr, "\n ii --- ii \n"); +#endif + twopoweri <<= 1; + } + for (loop= 0 ; loop < shift_x; loop++) + mpfr_mul(tmp,tmp,tmp,GMP_RNDD); + if (mpfr_can_round(tmp, realprec, GMP_RNDD, rnd_mode, PREC(y))){ + mpfr_set(y,tmp,rnd_mode); + good = 1; + } else { + realprec += 3*logn; + } + } + return 0; +} + + + + + + diff --git a/extract.c b/extract.c new file mode 100644 index 000000000..29655dba4 --- /dev/null +++ b/extract.c @@ -0,0 +1,60 @@ +#include <stdio.h> +#include <math.h> +#include <stdlib.h> + +#include "gmp.h" +#include "gmp-impl.h" +#include "mpfr.h" +#include "assert.h" + + +int +#if __STDC__ +mpfr_extract(mpz_ptr y, mpfr_srcptr p, int i) +#else +mpfr_extract(y,p,i) +mpz_ptr y; +mpfr_srcptr p; +int i; +#endif +{ + int two_i=1 << i; + int two_i_2 = i ? two_i / 2 : 1; + int size; + int j; + if (ABSSIZE(p) < two_i){ + int j; + y->_mp_alloc=two_i_2 ; + y->_mp_size=two_i_2 ; + y->_mp_d = calloc(two_i_2,sizeof( mp_limb_t)); + assert(y->_mp_d!=NULL); + /* attention : initialiser a 0 si on utilise malloc */ + for(j= 0; j < ABSSIZE(p) - two_i_2 ; j++){ + y->_mp_d[j + two_i - ABSSIZE(p)] = p->_mp_d[j]; + } + } else + { + y->_mp_d = p->_mp_d+ABSSIZE(p) - two_i; + y->_mp_alloc=two_i_2 ; + y->_mp_size=two_i_2 ; + } + /* if ISNEG(p) + CHANGE_SIGN(y);*/ + + size = ABSSIZE(y); + for (j = 0; j < size; j++) + { + if (y->_mp_d[j]) + { + if ISNEG(p) + mpz_neg(y,y); + return 0; + } + } + y->_mp_size=0; + + return 0; +} + + + |