fixed pb found by Damien Fisher

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

2 files changed, 34 insertions, 12 deletions

diff --git a/agm.c b/agm.c
index 9bfde443c..aeab2702e 100644
--- a/agm.c
+++ b/agm.c
@@ -117,14 +117,19 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode)
   while (go_on) {
     int err, can_round;
     mp_prec_t eq;
+    double erraux;
     
-    err=1 + (int) ((3.0/2.0*(double)__gmpfr_ceil_log2((double)p)+1.0)*__gmpfr_ceil_exp2(-(double)p)
-	     +3.0*__gmpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
+    erraux = (vo == uo) ? 0.0 : __gmpfr_ceil_exp2 (-2.0 * (double) p * uo 
+                                                   / (vo - uo));
+    err = 1 + (int) ((3.0 / 2.0 * (double) __gmpfr_ceil_log2 ((double) p) 
+                      + 1.0) * __gmpfr_ceil_exp2 (- (double) p)
+                     + 3.0 * erraux);
+
     if(p-err-3<=q) {
       p=q+err+4;
       err= 1 + 
 	(int) ((3.0/2.0*__gmpfr_ceil_log2((double)p)+1.0)*__gmpfr_ceil_exp2(-(double)p)
-	       +3.0*__gmpfr_ceil_exp2(-2.0*(double)p*uo/(vo-uo)));
+	       +3.0 * erraux);
     }
 
     /* Calculus of un and vn */
diff --git a/tests/tagm.c b/tests/tagm.c
index 57f5c6a37..100317a35 100644
--- a/tests/tagm.c
+++ b/tests/tagm.c
@@ -82,15 +82,32 @@ check_large (void)
 {
   mpfr_t a, b, agm;
 
-  mpfr_init2(a, 82); mpfr_init2(b, 82); mpfr_init2(agm, 82);
-  mpfr_set_ui(a, 1, GMP_RNDN);
-  mpfr_set_str_raw(b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39");
-  mpfr_agm(agm, a, b, GMP_RNDN);
-  mpfr_set_str_raw(a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4");
-  if (mpfr_cmp(agm, a)) {
-    fprintf(stderr, "mpfr_agm failed for precision 82\n"); exit(1);
-  }
-  mpfr_clear(a); mpfr_clear(b); mpfr_clear(agm);
+  mpfr_init2 (a, 82);
+  mpfr_init2 (b, 82);
+  mpfr_init2 (agm, 82);
+
+  mpfr_set_ui (a, 1, GMP_RNDN);
+  mpfr_set_str_raw (b, "0.1111101100001000000001011000110111101000001011111000100001000101010100011111110010E-39");
+  mpfr_agm (agm, a, b, GMP_RNDN);
+  mpfr_set_str_raw (a, "0.1110001000111101101010101010101101001010001001001011100101111011110101111001111100E-4");
+  if (mpfr_cmp (agm, a))
+    {
+      fprintf (stderr, "mpfr_agm failed for precision 82\n");
+      exit (1);
+    }
+  
+  /* problem found by Damien Fischer <damien@maths.usyd.edu.au> 4 Aug 2003:
+     produced a division by zero exception */
+  mpfr_set_prec (a, 268);
+  mpfr_set_prec (b, 268);
+  mpfr_set_prec (agm, 268);
+  mpfr_set_str (a, "703.93543315330225238487276503953366664991725089988315253092140138947103394917006", 10, GMP_RNDN);
+  mpfr_set_str (b, "703.93543315330225238487279020523738740563816490895994499256063816906728642622316", 10, GMP_RNDN);
+  mpfr_agm (agm, a, b, GMP_RNDN);
+
+  mpfr_clear (a);
+  mpfr_clear (b);
+  mpfr_clear (agm);
 }
 
 #if 0