fixed problem when te=ti=1 (i.e. t=0) found by Kevin Ryde

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

1 files changed, 43 insertions, 36 deletions

diff --git a/sinh.c b/sinh.c
index 88b3c458c..dceaf14c7 100644
--- a/sinh.c
+++ b/sinh.c
@@ -25,15 +25,15 @@ MA 02111-1307, USA. */
 #include "mpfr.h"
 #include "mpfr-impl.h"
 
- /* The computation of cosh is done by
+ /* The computation of sinh is done by
 
-    cosh= 1/2[e^(x)-e^(-x)]
+    sinh(x) = 1/2 [e^(x)-e^(-x)]
  */
 
 int
 mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
 {
-    /****** Declaration ******/
+    /****** Declarations ******/
     mpfr_t x;
     mp_prec_t Nxt = MPFR_PREC(xt);
     int flag_neg=0, inexact =0;
@@ -48,7 +48,7 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
     if (MPFR_IS_INF(xt))
       { 
         MPFR_SET_INF(y);
-        MPFR_SET_SAME_SIGN(y,xt);
+        MPFR_SET_SAME_SIGN(y, xt);
         MPFR_RET(0);
       }
 
@@ -57,12 +57,12 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
     if (MPFR_IS_ZERO(xt))
       {
         MPFR_SET_ZERO(y);   /* sinh(0) = 0 */
-        MPFR_SET_SAME_SIGN(y,xt);
+        MPFR_SET_SAME_SIGN(y, xt);
         MPFR_RET(0);
       }
 
-    mpfr_init2(x,Nxt);
-    mpfr_set(x,xt,GMP_RNDN);
+    mpfr_init2 (x, Nxt);
+    mpfr_set (x, xt, GMP_RNDN);
 
     if(MPFR_SIGN(x)<0)
       {
@@ -73,7 +73,7 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
     /* General case */
     {
     /* Declaration of the intermediary variable */
-      mpfr_t t, te,ti;       
+      mpfr_t t, te, ti;
       int d;
 
       /* Declaration of the size variable */
@@ -84,52 +84,59 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
       long int err;  /* Precision of error */
       
       /* compute the precision of intermediary variable */
-      Nt=MAX(Nx,Ny);
+      Nt = MAX(Nx, Ny);
       /* the optimal number of bits : see algorithms.ps */
-      Nt = Nt+_mpfr_ceil_log2(5)+_mpfr_ceil_log2(Nt);
+      Nt = Nt + _mpfr_ceil_log2 (5) + _mpfr_ceil_log2 (Nt);
 
       /* initialise of intermediary	variable */
-      mpfr_init(t);             
-      mpfr_init(te);             
-      mpfr_init(ti);                    
+      mpfr_init (t);
+      mpfr_init (te);
+      mpfr_init (ti);
 
-
-      /* First computation of cosh */
+      /* First computation of sinh */
       do {
 
 	/* reactualisation of the precision */
 
-	mpfr_set_prec(t,Nt);             
-	mpfr_set_prec(te,Nt);             
-	mpfr_set_prec(ti,Nt);             
+	mpfr_set_prec (t, Nt);
+	mpfr_set_prec (te, Nt);
+	mpfr_set_prec (ti, Nt);
 
 	/* compute sinh */
-	mpfr_exp(te,x,GMP_RNDD);         /* exp(x) */
-	mpfr_ui_div(ti,1,te,GMP_RNDU);   /* 1/exp(x) */
-	mpfr_sub(t,te,ti,GMP_RNDN);      /* exp(x) - 1/exp(x)*/
-	mpfr_div_2ui(t,t,1,GMP_RNDN);    /* 1/2(exp(x) - 1/exp(x))*/
-
-        /* calculation of the error*/
-        d = MPFR_EXP(te)-MPFR_EXP(t)+2;
+	mpfr_exp (te, x, GMP_RNDD);        /* exp(x) */
+	mpfr_ui_div (ti, 1, te, GMP_RNDU); /* 1/exp(x) */
+        mpfr_sub (t, te, ti, GMP_RNDN);    /* exp(x) - 1/exp(x) */
+	mpfr_div_2ui (t, t, 1, GMP_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */
+
+        /* it may be that t is zero (in fact, it can only occur when te=1,
+           and thus ti=1 too) */
+
+        if (MPFR_IS_ZERO(t))
+          err = -1;
+        else
+          {
+            /* calculation of the error */
+            d = MPFR_EXP(te) - MPFR_EXP(t) + 2;
 	
-	/* estimation of the error */
-        /* err = Nt-(_mpfr_ceil_log2(1+pow(2,d)));*/
-        err = Nt-(MAX(d,0)+1);
+            /* estimation of the error */
+            /* err = Nt-(_mpfr_ceil_log2(1+pow(2,d)));*/
+            err = Nt - (MAX(d,0) + 1);
+          }
 
 	/* actualisation of the precision */
         Nt += 10; 
 
-      } while ((err < 0) || !mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny));
+      } while ((err < 0) || !mpfr_can_round(t, err, GMP_RNDN, rnd_mode, Ny));
 
-      if (flag_neg==1)
+      if (flag_neg == 1)
           MPFR_CHANGE_SIGN(t);
 
-      inexact = mpfr_set(y,t,rnd_mode);
-      mpfr_clear(t);
-      mpfr_clear(ti);
-      mpfr_clear(te);
+      inexact = mpfr_set (y, t, rnd_mode);
+      mpfr_clear (t);
+      mpfr_clear (ti);
+      mpfr_clear (te);
     }
-    mpfr_clear(x);
-    return inexact;
+    mpfr_clear (x);
 
+    return inexact;
 }