acosh revue et corrige...+beau

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

1 files changed, 77 insertions, 95 deletions

diff --git a/acosh.c b/acosh.c
index 1a7f88a49..92f31bf85 100644
--- a/acosh.c
+++ b/acosh.c
@@ -20,6 +20,7 @@ the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
 #include <limits.h>
+#include <math.h>
 #include "gmp.h"
 #include "gmp-impl.h"
 #include "mpfr.h"
@@ -43,109 +44,90 @@ mpfr_acosh (y, x, rnd_mode)
 #endif
 {
     
-  /****** Declaration ******/
-
-    /* Variable of Intermediary Calculation*/
-    mpfr_t t;       
-
-    /* Variable of Intermediary Calculation*/
-    mpfr_t ta,tb;
-
-    int round;
-    int boucle;
-    int comp;
-
-    mp_prec_t Nx;   /* Precision of input variable */
-    mp_prec_t Ny;   /* Precision of output variable */
-    mp_prec_t Nt;   /* Precision of Intermediary Calculation variable */
-    mp_prec_t err;  /* Precision of error */
+  int inexact =0;
+  int comp;
 
+  if (MPFR_IS_NAN(x)) 
+    {
+      MPFR_SET_NAN(y); 
+      return 1; 
+    }
     
-    if (MPFR_IS_NAN(x)) {  MPFR_SET_NAN(y); return 1; }
-    MPFR_CLEAR_NAN(y);
+  comp=mpfr_cmp_ui(x,1);
+
+  if(comp < 0)
+    {
+      MPFR_SET_NAN(y); 
+      return(1);
+    }
+  MPFR_CLEAR_NAN(y);
 
-    if (MPFR_IS_INF(x)){ 
+  if(comp == 0)
+    {
+      MPFR_SET_ZERO(y); /* acosh(1) = 0 */
+      return(0);
+    }
+  
+  if (MPFR_IS_INF(x))
+    { 
       MPFR_SET_INF(y);
-      if (MPFR_SIGN(y) < 0) MPFR_CHANGE_SIGN(y);
+      if (MPFR_SIGN(y) < 0) 
+        MPFR_CHANGE_SIGN(y);
       return 1;
     }
 
-    MPFR_CLEAR_INF(y);
+  MPFR_CLEAR_INF(y);
 
-    comp=mpfr_cmp_ui(x,1);
-
-    if(comp <= 0){
-      if(comp == 0){
-	MPFR_SET_ZERO(y); /* acosh(1) = 0 */
-	return(0);
-      }
-      else{
-	/*fprintf(stderr,"Function acosh of MPFR is only defined for x=[1,+Inf]");
-	exit(-1);*/
-
-	/*An other strategy if output is not define for input return NaN*/
-	MPFR_SET_NAN(y); return(-1);
-
-      }
-    }
-    else{
-      
-
-        /* Initialisation of the Precision */
-	Nx=MPFR_PREC(x);
-	Ny=MPFR_PREC(y);
-
-	/* compute the size of intermediary variable */
-	if(Ny>=Nx)
-	  Nt=Ny+2*CHAR_BIT;
-	else
-	  Nt=Nx+2*CHAR_BIT;
-	  
-	  boucle=1;
-
-	    /* initialise of intermediary	variable */
-	    mpfr_init2(t,Nt);                
-	    mpfr_init2(ta,Nt);             
-	    mpfr_init2(tb,Nt);             
-
-	  while(boucle==1){
-
-
-	    /* compute acosh */
-
-	    mpfr_add_ui(ta,x,1,GMP_RNDN);  /* (x+1) */
-	    mpfr_sub_ui(tb,x,1,GMP_RNDN);  /* (x-1) */
-	    mpfr_sqrt(ta,ta,GMP_RNDN);     /* sqrt(x+1) */
-	    mpfr_sqrt(tb,tb,GMP_RNDN);     /* sqrt(x-1) */
-	    mpfr_mul(t,ta,tb,GMP_RNDN);    /* sqrt(x+1)*sqrt(x-1) */
-	    mpfr_add(t,t,x,GMP_RNDN);      /* sqrt(x+1)*sqrt(x-1)+x */
-	    mpfr_log(t,t,GMP_RNDN);        /* ln(sqrt(x+1)*sqrt(x-1)+x)*/
-	    
-	    err=Nt-1-MAX(0,-MPFR_EXP(t));
-
-	    round=mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny);
-
-	    if(round == 1){
-	      mpfr_set(y,t,rnd_mode);
-	      boucle=0;
-	    }
-	    else{
-	      Nt=Nt+10;
-	      /* initialise of intermediary	variable */
-	      mpfr_set_prec(t,Nt);                
-	      mpfr_set_prec(ta,Nt);   
-	      mpfr_set_prec(tb,Nt);             
-
-	      boucle=1;
-	    }
-	    
-	  }
+  /* General case */
+  {
+    /* Declaration of the intermediary variable */
+    mpfr_t t, te,ti;       
     
-	  mpfr_clear(t);
-	  mpfr_clear(ta);
-	  mpfr_clear(tb);
-          return(1);
-      
+    /* Declaration of the size variable */
+    mp_prec_t Nx = MPFR_PREC(x);   /* Precision of input variable */
+    mp_prec_t Ny = MPFR_PREC(y);   /* Precision of input variable */
+    
+    mp_prec_t Nt;   /* Precision of the intermediary variable */
+    mp_prec_t err;  /* Precision of error */
+                
+    /* compute the precision of intermediary variable */
+    Nt=MAX(Nx,Ny);
+    /* the optimal number of bits : see algorithms.ps */
+    Nt=Nt+4+_mpfr_ceil_log2(Nt);
+
+    /* initialise of intermediary	variable */
+    mpfr_init(t);             
+    mpfr_init(te);             
+    mpfr_init(ti);                    
+
+    /* First computation of cosh */
+    do {
+
+      /* reactualisation of the precision */
+      mpfr_set_prec(t,Nt);             
+      mpfr_set_prec(te,Nt);             
+      mpfr_set_prec(ti,Nt);             
+
+      /* compute acosh */
+      mpfr_mul(te,x,x,GMP_RNDD);  /* (x^2) */
+      mpfr_sub_ui(ti,te,1,GMP_RNDD);  /* (x^2-1) */
+      mpfr_sqrt(t,ti,GMP_RNDN);     /* sqrt(x^2-1) */
+      mpfr_add(t,t,x,GMP_RNDN);    /* sqrt(x^2-1)+x */
+      mpfr_log(t,t,GMP_RNDN);        /* ln(sqrt(x^2-1)+x)*/
+
+      /* estimation of the error see- algorithms.ps*/
+      err=Nt-_mpfr_ceil_log2(0.5+pow(2,2-MPFR_EXP(t))+pow(2,1+MPFR_EXP(te)-MPFR_EXP(ti)-MPFR_EXP(t)));
+
+      /* actualisation of the precision */
+      Nt += 10;
+
+    } while (!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny));
  
-    }
+    inexact = mpfr_set(y,t,rnd_mode);
+
+    mpfr_clear(t);
+    mpfr_clear(ti);
+    mpfr_clear(te);
+  }
+  return inexact;
 }