[src/zeta_ui.c]

  * Support reduced exponent range for the generic case.
  * Added logging.
[tests/tzeta_ui.c] Added tests in reduced exponent range.
(merged changesets r9518-9523 from the trunk, with mpfr_flags_t
replaced by unsigned int)

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/3.1@9524 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 131 insertions, 29 deletions

diff --git a/src/zeta_ui.c b/src/zeta_ui.c
index 32350fdfe..02ba337fb 100644
--- a/src/zeta_ui.c
+++ b/src/zeta_ui.c
@@ -28,12 +28,13 @@ mpfr_zeta_ui (mpfr_ptr z, unsigned long m, mpfr_rnd_t r)
 {
   MPFR_ZIV_DECL (loop);
 
+  MPFR_LOG_FUNC
+    (("m=%lu rnd=%d prec=%Pu", m, r, mpfr_get_prec (z)),
+     ("z[%Pu]=%.*Rg", mpfr_get_prec (z), mpfr_log_prec, z));
+
   if (m == 0)
     {
-      mpfr_set_ui (z, 1, r);
-      mpfr_div_2ui (z, z, 1, r);
-      MPFR_CHANGE_SIGN (z);
-      MPFR_RET (0);
+      return mpfr_set_si_2exp (z, -1, -1, r);
     }
   else if (m == 1)
     {
@@ -49,28 +50,32 @@ mpfr_zeta_ui (mpfr_ptr z, unsigned long m, mpfr_rnd_t r)
       mpz_t d, t, s, q;
       mpfr_t y;
       int inex;
+      MPFR_SAVE_EXPO_DECL (expo);
 
       if (r == MPFR_RNDA)
         r = MPFR_RNDU; /* since the result is always positive */
 
+      MPFR_SAVE_EXPO_MARK (expo);
+
       if (m >= p) /* 2^(-m) < ulp(1) = 2^(1-p). This means that
                      2^(-m) <= 1/2*ulp(1). We have 3^(-m)+4^(-m)+... < 2^(-m)
                      i.e. zeta(m) < 1+2*2^(-m) for m >= 3 */
-
         {
           if (m == 2) /* necessarily p=2 */
-            return mpfr_set_ui_2exp (z, 13, -3, r);
-          else if (r == MPFR_RNDZ || r == MPFR_RNDD || (r == MPFR_RNDN && m > p))
+            inex = mpfr_set_ui_2exp (z, 13, -3, r);
+          else if (r == MPFR_RNDZ || r == MPFR_RNDD ||
+                   (r == MPFR_RNDN && m > p))
             {
               mpfr_set_ui (z, 1, r);
-              return -1;
+              inex = -1;
             }
           else
             {
               mpfr_set_ui (z, 1, r);
               mpfr_nextabove (z);
-              return 1;
+              inex = 1;
             }
+          goto end;
         }
 
       /* now treat also the case where zeta(m) - (1+1/2^m) < 1/2*ulp(1),
@@ -89,9 +94,13 @@ mpfr_zeta_ui (mpfr_ptr z, unsigned long m, mpfr_rnd_t r)
               mpfr_div_2ui (z, z, m, MPFR_RNDZ);
               mpfr_add_ui (z, z, 1, MPFR_RNDZ);
               if (r != MPFR_RNDU)
-                return -1;
-              mpfr_nextabove (z);
-              return 1;
+                inex = -1;
+              else
+                {
+                  mpfr_nextabove (z);
+                  inex = 1;
+                }
+              goto end;
             }
         }
 
@@ -225,6 +234,11 @@ mpfr_zeta_ui (mpfr_ptr z, unsigned long m, mpfr_rnd_t r)
       mpz_clear (s);
       inex = mpfr_set (z, y, r);
       mpfr_clear (y);
-      return inex;
+
+    end:
+      MPFR_LOG_VAR (z);
+      MPFR_LOG_MSG (("inex = %d before mpfr_check_range\n", inex));
+      MPFR_SAVE_EXPO_FREE (expo);
+      return mpfr_check_range (z, inex, r);
     }
 }
diff --git a/tests/tzeta_ui.c b/tests/tzeta_ui.c
index d86de2aed..8d8c2274b 100644
--- a/tests/tzeta_ui.c
+++ b/tests/tzeta_ui.c
@@ -36,9 +36,15 @@ main (int argc, char *argv[])
   mpfr_t x, y, z, t;
   unsigned long n;
   int inex;
+  mpfr_exp_t emin, emax;
+  unsigned int flags, ex_flags;
+  int i;
 
   tests_start_mpfr ();
 
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+
   mpfr_init (x);
   mpfr_init (y);
   mpfr_init (z);
@@ -66,10 +72,73 @@ main (int argc, char *argv[])
       exit (1);
     }
 
-  mpfr_clear_divby0 ();
+  mpfr_clear_flags ();
   inex = mpfr_zeta_ui (x, 0, MPFR_RNDN);
-  MPFR_ASSERTN (inex == 0 && mpfr_cmp_si_2exp (x, -1, -1) == 0
-                && !mpfr_divby0_p ());
+  flags = __gmpfr_flags;
+  MPFR_ASSERTN (inex == 0 && mpfr_cmp_si_2exp (x, -1, -1) == 0 && flags == 0);
+
+  for (i = -2; i <= 2; i += 2)
+    RND_LOOP (rnd)
+      {
+        int ex_inex;
+
+        set_emin (i);
+        set_emax (i);
+        mpfr_clear_flags ();
+        inex = mpfr_zeta_ui (x, 0, (mpfr_rnd_t) rnd);
+        flags = __gmpfr_flags;
+        if (i < 0)
+          {
+            mpfr_set_inf (y, -1);
+            if (rnd == MPFR_RNDU || rnd == MPFR_RNDZ)
+              {
+                mpfr_nextabove (y);
+                ex_inex = 1;
+              }
+            else
+              {
+                ex_inex = -1;
+              }
+            ex_flags = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
+          }
+        else if (i > 0)
+          {
+            mpfr_set_zero (y, -1);
+            if (rnd == MPFR_RNDD || rnd == MPFR_RNDA)
+              {
+                mpfr_nextbelow (y);
+                ex_inex = -1;
+              }
+            else
+              {
+                ex_inex = 1;
+              }
+            ex_flags = MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT;
+          }
+        else
+          {
+            mpfr_set_str_binary (y, "-1e-1");
+            ex_inex = 0;
+            ex_flags = 0;
+          }
+        set_emin (emin);
+        set_emax (emax);
+        if (! (mpfr_equal_p (x, y) && MPFR_IS_NEG (x) &&
+               SAME_SIGN (inex, ex_inex) && flags == ex_flags))
+          {
+            printf ("Failure for zeta(0) in %s, exponent range [%d,%d]\n",
+                    mpfr_print_rnd_mode ((mpfr_rnd_t) rnd), i, i);
+            printf ("Expected ");
+            mpfr_dump (y);
+            printf ("  with inex ~ %d, flags =", ex_inex);
+            flags_out (ex_flags);
+            printf ("Got      ");
+            mpfr_dump (x);
+            printf ("  with inex = %d, flags =", inex);
+            flags_out (flags);
+            exit (1);
+          }
+      }
 
   mpfr_clear_divby0 ();
   inex = mpfr_zeta_ui (x, 1, MPFR_RNDN);
@@ -85,26 +154,45 @@ main (int argc, char *argv[])
       mpfr_set_prec (y, yprec);
 
       for (n = 0; n < 50; n++)
-        for (rnd = 0; rnd < MPFR_RND_MAX; rnd++)
+        RND_LOOP (rnd)
           {
             mpfr_zeta_ui (y, n, MPFR_RNDN);
             if (mpfr_can_round (y, yprec, MPFR_RNDN, MPFR_RNDZ, prec
                                 + (rnd == MPFR_RNDN)))
               {
                 mpfr_set (t, y, (mpfr_rnd_t) rnd);
-                mpfr_zeta_ui (z, n, (mpfr_rnd_t) rnd);
-                if (mpfr_cmp (t, z))
+                for (i = 0; i <= 1; i++)
                   {
-                    printf ("results differ for n=%lu", n);
-                    printf (" prec=%u rnd_mode=%s\n", prec,
-                            mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
-                    printf ("   got      ");
-                    mpfr_dump (z);
-                    printf ("   expected ");
-                    mpfr_dump (t);
-                    printf ("   approx   ");
-                    mpfr_dump (y);
-                    exit (1);
+                    if (i)
+                      {
+                        mpfr_exp_t e;
+
+                        if (MPFR_IS_SINGULAR (t))
+                          break;
+                        e = mpfr_get_exp (t);
+                        set_emin (e);
+                        set_emax (e);
+                      }
+                    mpfr_zeta_ui (z, n, (mpfr_rnd_t) rnd);
+                    if (i)
+                      {
+                        set_emin (emin);
+                        set_emax (emax);
+                      }
+                    if (mpfr_cmp (t, z))
+                      {
+                        printf ("results differ for n = %lu, prec = %u,"
+                                " %s%s\n", n, prec,
+                                mpfr_print_rnd_mode ((mpfr_rnd_t) rnd),
+                                i ? ", reduced exponent range" : "");
+                        printf ("  got      ");
+                        mpfr_dump (z);
+                        printf ("  expected ");
+                        mpfr_dump (t);
+                        printf ("  approx   ");
+                        mpfr_dump (y);
+                        exit (1);
+                      }
                   }
               }
           }