First Release

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@653 280ebfd0-de03-0410-8827-d642c229c3f4

4 files changed, 333 insertions, 106 deletions

diff --git a/exp.c b/exp.c
index 67b4e1b7a..a8f939e6e 100644
--- a/exp.c
+++ b/exp.c
@@ -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;
 }
 
diff --git a/exp2.c b/exp2.c
index af61906dc..500fda0fc 100644
--- a/exp2.c
+++ b/exp2.c
@@ -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);
diff --git a/exp3.c b/exp3.c
new file mode 100644
index 000000000..35ad20d26
--- /dev/null
+++ b/exp3.c
@@ -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;
+}
+
+
+