[src/jyn_asympt.c] partly solved the slowness of jn(733333,733333)

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

2 files changed, 31 insertions, 7 deletions

diff --git a/src/jyn_asympt.c b/src/jyn_asympt.c
index ff387a264..692327dfb 100644
--- a/src/jyn_asympt.c
+++ b/src/jyn_asympt.c
@@ -49,7 +49,18 @@ FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
   mpfr_exp_t err2, err;
   MPFR_ZIV_DECL (loop);
 
-  mpfr_init (c);
+  mpfr_init2 (c, 64);
+
+  /* The terms of the asymptotic expansion grow like mu^(2k)/(8z)^(2k), where
+     mu = 4n^2, thus we need mu < 8z so that it converges, i.e., n^2/2 < z. */
+  mpfr_set_ui (c, n, MPFR_RNDU);
+  mpfr_mul_ui (c, c, n, MPFR_RNDU);
+  mpfr_div_2ui (c, c, 1, MPFR_RNDU);
+  if (mpfr_cmp (c, z) >= 0)
+    {
+      mpfr_clear (c);
+      return 0; /* asymptotic expansion failed */
+    }
 
   w = MPFR_PREC(res) + MPFR_INT_CEIL_LOG2(MPFR_PREC(res)) + 4;
 
diff --git a/tests/tjn.c b/tests/tjn.c
index d17b81020..40146a207 100644
--- a/tests/tjn.c
+++ b/tests/tjn.c
@@ -22,6 +22,22 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
 
 #include "mpfr-test.h"
 
+/* mpfr_jn doesn't terminate. Bug reported by Alex Coplan on 2020-07-03.
+* https://gcc.gnu.org/bugzilla/show_bug.cgi?id=96044
+*/
+static void
+bug20200703 (void)
+{
+  mpfr_t x, y;
+
+  mpfr_init (x);
+  mpfr_init (y);
+  mpfr_set_si (x, 733333, MPFR_RNDN);
+  mpfr_jn (y, 733333, x, MPFR_RNDN);
+  mpfr_clear (x);
+  mpfr_clear (y);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -41,6 +57,9 @@ main (int argc, char *argv[])
 
   tests_start_mpfr ();
 
+  if (getenv ("MPFR_CHECK_EXPENSIVE") != NULL)
+    bug20200703 ();
+
   mpfr_init (x);
   mpfr_init (y);
 
@@ -294,12 +313,6 @@ main (int argc, char *argv[])
       exit (1);
     }
 
-  /* mpfr_jn doesn't terminate. Bug reported by Alex Coplan on 2020-07-03.
-   * https://gcc.gnu.org/bugzilla/show_bug.cgi?id=96044
-   */
-  mpfr_set_si (x, 733333, MPFR_RNDN);
-  mpfr_jn (y, 733333, x, MPFR_RNDN);
-
   mpfr_clear (x);
   mpfr_clear (y);