File _service:download_src_package:0005-digamma-huge... of Package mpfr

File _service:download_src_package:0005-digamma-hugemem.patch of Package mpfr

diff -Naurd mpfr-4.1.0-a/PATCHES mpfr-4.1.0-b/PATCHES
--- mpfr-4.1.0-a/PATCHES	2021-02-11 12:50:22.384987438 +0000
+++ mpfr-4.1.0-b/PATCHES	2021-02-11 12:50:22.424987002 +0000
@@ -0,0 +1 @@
+digamma-hugemem
diff -Naurd mpfr-4.1.0-a/VERSION mpfr-4.1.0-b/VERSION
--- mpfr-4.1.0-a/VERSION	2021-02-11 12:48:27.370242746 +0000
+++ mpfr-4.1.0-b/VERSION	2021-02-11 12:50:22.424987002 +0000
@@ -1 +1 @@
-4.1.0-p4
+4.1.0-p5
diff -Naurd mpfr-4.1.0-a/src/digamma.c mpfr-4.1.0-b/src/digamma.c
--- mpfr-4.1.0-a/src/digamma.c	2020-06-18 17:17:18.000000000 +0000
+++ mpfr-4.1.0-b/src/digamma.c	2021-02-11 12:50:22.412987133 +0000
@@ -214,19 +214,27 @@
     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
      ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex));
 
-  /* compute a precision q such that x+1 is exact */
-  if (MPFR_PREC(x) < MPFR_GET_EXP(x))
-    q = MPFR_EXP(x);
-  else
-    q = MPFR_PREC(x) + 1;
-
-  /* for very large x, use |digamma(x) - log(x)| < 1/x < 2^(1-EXP(x)) */
-  if (MPFR_PREC(y) + 10 < MPFR_EXP(x))
+  /* For very large x, use |digamma(x) - log(x)| < 1/x < 2^(1-EXP(x)).
+     However, for a fixed value of GUARD, MPFR_CAN_ROUND() might fail
+     with probability 1/2^GUARD, in which case the default code will
+     fail since it requires x+1 to be exact, thus a huge precision if
+     x is huge. There are two workarounds:
+     * either perform a Ziv's loop, by increasing GUARD at each step.
+       However, this might fail if x is moderately large, in which case
+       more terms of the asymptotic expansion would be needed.
+     * implement a full asymptotic expansion (with Ziv's loop). */
+#define GUARD 30
+  if (MPFR_PREC(y) + GUARD < MPFR_EXP(x))
     {
       /* this ensures EXP(x) >= 3, thus x >= 4, thus log(x) > 1 */
-      mpfr_init2 (t, MPFR_PREC(y) + 10);
-      mpfr_log (t, x, MPFR_RNDZ);
-      if (MPFR_CAN_ROUND (t, MPFR_PREC(y) + 10, MPFR_PREC(y), rnd_mode))
+      mpfr_init2 (t, MPFR_PREC(y) + GUARD);
+      mpfr_log (t, x, MPFR_RNDN);
+      /* |t - digamma(x)| <= 1/2*ulp(t) + |digamma(x) - log(x)|
+                          <= 1/2*ulp(t) + 2^(1-EXP(x))
+                          <= 1/2*ulp(t) + 2^(-PREC(y)-GUARD)
+                          <= ulp(t)
+         since |t| >= 1 thus ulp(t) >= 2^(1-PREC(y)-GUARD) */
+      if (MPFR_CAN_ROUND (t, MPFR_PREC(y) + GUARD, MPFR_PREC(y), rnd_mode))
         {
           inex = mpfr_set (y, t, rnd_mode);
           mpfr_clear (t);
@@ -235,6 +243,21 @@
       mpfr_clear (t);
     }
 
+  /* compute a precision q such that x+1 is exact */
+  if (MPFR_PREC(x) < MPFR_GET_EXP(x))
+    {
+      /* The goal of the first assertion is to let the compiler ignore
+         the second one when MPFR_EMAX_MAX <= MPFR_PREC_MAX. */
+      MPFR_ASSERTD (MPFR_EXP(x) <= MPFR_EMAX_MAX);
+      MPFR_ASSERTN (MPFR_EXP(x) <= MPFR_PREC_MAX);
+      q = MPFR_EXP(x);
+    }
+  else
+    q = MPFR_PREC(x) + 1;
+
+  /* FIXME: q can be much too large, e.g. equal to the maximum exponent! */
+  MPFR_LOG_MSG (("q=%Pu\n", q));
+
   mpfr_init2 (x_plus_j, q);
 
   mpfr_init2 (t, p);
diff -Naurd mpfr-4.1.0-a/src/mpfr.h mpfr-4.1.0-b/src/mpfr.h
--- mpfr-4.1.0-a/src/mpfr.h	2021-02-11 12:48:27.366242791 +0000
+++ mpfr-4.1.0-b/src/mpfr.h	2021-02-11 12:50:22.424987002 +0000
@@ -27,7 +27,7 @@
 #define MPFR_VERSION_MAJOR 4
 #define MPFR_VERSION_MINOR 1
 #define MPFR_VERSION_PATCHLEVEL 0
-#define MPFR_VERSION_STRING "4.1.0-p4"
+#define MPFR_VERSION_STRING "4.1.0-p5"
 
 /* User macros:
    MPFR_USE_FILE:        Define it to make MPFR define functions dealing
diff -Naurd mpfr-4.1.0-a/src/version.c mpfr-4.1.0-b/src/version.c
--- mpfr-4.1.0-a/src/version.c	2021-02-11 12:48:27.370242746 +0000
+++ mpfr-4.1.0-b/src/version.c	2021-02-11 12:50:22.424987002 +0000
@@ -25,5 +25,5 @@
 const char *
 mpfr_get_version (void)
 {
-  return "4.1.0-p4";
+  return "4.1.0-p5";
 }
diff -Naurd mpfr-4.1.0-a/tests/tdigamma.c mpfr-4.1.0-b/tests/tdigamma.c
--- mpfr-4.1.0-a/tests/tdigamma.c	2020-06-18 17:17:18.000000000 +0000
+++ mpfr-4.1.0-b/tests/tdigamma.c	2021-02-11 12:50:22.412987133 +0000
@@ -49,12 +49,54 @@
   mpfr_clear (y);
 }
 
+/* With some GMP_CHECK_RANDOMIZE values, test_generic triggers an error
+     tests_addsize(): too much memory (576460752303432776 bytes)
+  Each time on prec = 200, n = 3, xprec = 140.
+  The following test is a more general testcase.
+*/
+static void
+bug20210206 (void)
+{
+#define NPREC 4
+  mpfr_t x, y[NPREC], z;
+  mpfr_exp_t emin, emax;
+  int i, precx, precy[NPREC] = { 200, 400, 520, 1416 };
+
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+  set_emin (MPFR_EMIN_MIN);
+  set_emax (MPFR_EMAX_MAX);
+
+  for (i = 0; i < NPREC; i++)
+    mpfr_init2 (y[i], precy[i]);
+  mpfr_init2 (z, precy[0]);
+
+  for (precx = MPFR_PREC_MIN; precx < 150; precx++)
+    {
+      mpfr_init2 (x, precx);
+      mpfr_setmax (x, __gmpfr_emax);
+      for (i = 0; i < NPREC; i++)
+        mpfr_digamma (y[i], x, MPFR_RNDA);
+      mpfr_set (z, y[1], MPFR_RNDA);
+      MPFR_ASSERTN (mpfr_equal_p (y[0], z));
+      mpfr_clear (x);
+    }
+
+  for (i = 0; i < NPREC; i++)
+    mpfr_clear (y[i]);
+  mpfr_clear (z);
+
+  set_emin (emin);
+  set_emax (emax);
+}
+
 int
 main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
   special ();
+  bug20210206 ();
 
   test_generic (MPFR_PREC_MIN, 200, 20);