Fixed bug found by Fredrik Johansson in the mpfr_jn family of functions.

(merged changesets r9841-9842,9844 from the trunk)

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

3 files changed, 17 insertions, 5 deletions

diff --git a/doc/algorithms.tex b/doc/algorithms.tex
index 416d3d53d..3b987b0ad 100644
--- a/doc/algorithms.tex
+++ b/doc/algorithms.tex
@@ -3472,10 +3472,10 @@ If we first approximate $\chi = z - (n/2 + 1/4) \pi$ with working precision
 $w$, and then approximate $\cos \chi$ and $\sin \chi$, there will be a huge
 relative error if $z > 2^w$. Instead, we use the fact that for $n$ even,
 \[ P(n,z) \cos \chi - Q(n,z) \sin \chi = \frac{1}{\sqrt{2}} (-1)^{n/2}
-   [P (\sin z + \cos z) + Q (\cos z - \sin z)], \]
+   [P(n,z) (\sin z + \cos z) + Q(n,z) (\cos z - \sin z)], \]
 and for $n$ odd,
 \[ P(n,z) \cos \chi - Q(n,z) \sin \chi = \frac{1}{\sqrt{2}} (-1)^{(n-1)/2}
-   [P (\sin z - \cos z) + Q (\cos z + \sin z)], \]
+   [P(n,z) (\sin z - \cos z) + Q(n,z) (\cos z + \sin z)], \]
 where $\cos z$ and $\sin z$ are computed accurately with
 \texttt{mpfr\_sin\_cos}, which uses in turn \texttt{mpfr\_remainder}.
 
diff --git a/src/jyn_asympt.c b/src/jyn_asympt.c
index 89f684689..a97d487ec 100644
--- a/src/jyn_asympt.c
+++ b/src/jyn_asympt.c
@@ -253,9 +253,9 @@ FUNCTION (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
         break;
       if (diverge != 0)
         {
-          mpfr_set (c, z, r); /* will force inex=0 below, which means the
-                               asymptotic expansion failed */
-          break;
+          MPFR_ZIV_FREE (loop);
+          mpfr_clear (c);
+          return 0; /* means that the asymptotic expansion failed */
         }
       MPFR_ZIV_NEXT (loop, w);
     }
diff --git a/tests/tj0.c b/tests/tj0.c
index 66bb65a91..4e83166b4 100644
--- a/tests/tj0.c
+++ b/tests/tj0.c
@@ -99,6 +99,18 @@ main (int argc, char *argv[])
   mpfr_j0 (y, x, MPFR_RNDN);
   MPFR_ASSERTN (! mpfr_nan_p (y) && mpfr_cmp_ui_2exp (y, 41, -11) == 0);
 
+  /* Bug reported by Fredrik Johansson on 19 Jan 2016 */
+  mpfr_set_prec (x, 53);
+  mpfr_set_str (x, "0x4.3328p+0", 0, MPFR_RNDN);
+  mpfr_set_prec (y, 2);
+  mpfr_j0 (y, x, MPFR_RNDD);
+  /* y should be -0.5 */
+  MPFR_ASSERTN (! mpfr_nan_p (y) && mpfr_cmp_si_2exp (y, -1, -1) == 0);
+  mpfr_set_prec (y, 3);
+  mpfr_j0 (y, x, MPFR_RNDD);
+  /* y should be -0.4375 */
+  MPFR_ASSERTN (! mpfr_nan_p (y) && mpfr_cmp_si_2exp (y, -7, -4) == 0);
+
   /* Case for which s = 0 in mpfr_jn */
   mpfr_set_prec (x, 44);
   mpfr_set_prec (y, 44);