Fix jn() precision problems around zero points of j0()

There appears to be a really nasty bug in jn() from fdlibm, which is
the foundation for most libm implementations (including glibc libm).
The zeroth-order j0() and first-order j1() cylindrical Bessel functions
are used to recursively generate the jn() value, but only the zeroth-order
Bessel function is used to normalize it; however, each of the functions
gets highly imprecise (approaching "bogus") near its zero point, making
the jn() value itself bogus.

But in fact, the zero points of j0() and j1() never coincide, thus j1()
should be used in case it is more precise than j0(). (That is, simply
when its value is further from zero.)

As an example, 2.4048255576957729_8 is the first zero of j0().
The proper value as calculated by Mathematica is 0.19899990535769...
However, jn() returns -inf on 64-bit arch, or 0.185007 on 32-bit arch.
With the proposed patch below, the returned value is 0.199000.

The fix is based on work by Steve Kargl and Tobias Burnus.

diff --git a/ChangeLog b/ChangeLog
index 601f6f0f433b9c1b1996aa25108dede8975a7658..122a2b874343b6cf1fd0f756903961100d09fcbf 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+2010-05-12  Petr Baudis <pasky@suse.cz>
+
+       [BZ #11589]
+       * sysdeps/ieee754/dbl-64/e_jn.c: Compensate major precision loss
+       around j0() zero points by switching to j1().
+       * sysdeps/ieee754/flt-32/e_jnf.c: Likewise.
+       * sysdeps/ieee754/ldbl-128/e_jnl.c: Likewise.
+       * sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Likewise.
+       * sysdeps/ieee754/ldbl-96/e_jnl.c: Likewise.
+
 2010-06-01  Jan Sembera  <jsembera@suse.cz>
 
        [BZ #6812]

diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c
index bf4a13d974c353c46384aad3bfa4e3c200114634..d9d6f917626347fc85b0b97448ffadbaad870c0e 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_jn.c
+++ b/sysdeps/ieee754/dbl-64/e_jn.c
@@ -215,7 +215,16 @@ static double zero  =  0.00000000000000000000e+00;
                        }
                    }
                }
-               b = (t*__ieee754_j0(x)/b);
+               /* j0() and j1() suffer enormous loss of precision at and
+                * near zero; however, we know that their zero points never
+                * coincide, so just choose the one further away from zero.
+                */
+               z = __ieee754_j0 (x);
+               w = __ieee754_j1 (x);
+               if (fabs (z) >= fabs (w))
+                 b = (t * z / b);
+               else
+                 b = (t * w / a);
            }
        }
        if(sgn==1) return -b; else return b;

diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c
index de2e53de83834d371cc6d0917e10d6f71bf0ede6..dd3d551a39ef6cdd55e4b8ad83978acc711a52d4 100644 (file)
--- a/sysdeps/ieee754/flt-32/e_jnf.c
+++ b/sysdeps/ieee754/flt-32/e_jnf.c
@@ -165,7 +165,16 @@ static float zero  =  0.0000000000e+00;
                        }
                    }
                }
-               b = (t*__ieee754_j0f(x)/b);
+               /* j0() and j1() suffer enormous loss of precision at and
+                * near zero; however, we know that their zero points never
+                * coincide, so just choose the one further away from zero.
+                */
+               z = __ieee754_j0f (x);
+               w = __ieee754_j1f (x);
+               if (fabsf (z) >= fabsf (w))
+                 b = (t * z / b);
+               else
+                 b = (t * w / a);
            }
        }
        if(sgn==1) return -b; else return b;

diff --git a/sysdeps/ieee754/ldbl-128/e_jnl.c b/sysdeps/ieee754/ldbl-128/e_jnl.c
index a4a4e24cf5a3afdc41dc3ec550b4f88fb76d4470..a7f6772c08f27fda43f866871356b8334cad0ea4 100644 (file)
--- a/sysdeps/ieee754/ldbl-128/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-128/e_jnl.c
@@ -285,7 +285,16 @@ __ieee754_jnl (n, x)
                    }
                }
            }
-         b = (t * __ieee754_j0l (x) / b);
+         /* j0() and j1() suffer enormous loss of precision at and
+          * near zero; however, we know that their zero points never
+          * coincide, so just choose the one further away from zero.
+          */
+         z = __ieee754_j0l (x);
+         w = __ieee754_j1l (x);
+         if (fabsl (z) >= fabsl (w))
+           b = (t * z / b);
+         else
+           b = (t * w / a);
        }
     }
   if (sgn == 1)

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
index 0eea7456761265b3911bbf95f28963fea0aed487..372f942bfc489eebee442b3d98ecfb079d6f288b 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
@@ -286,7 +286,16 @@ __ieee754_jnl (n, x)
                    }
                }
            }
-         b = (t * __ieee754_j0l (x) / b);
+         /* j0() and j1() suffer enormous loss of precision at and
+          * near zero; however, we know that their zero points never
+          * coincide, so just choose the one further away from zero.
+          */
+         z = __ieee754_j0l (x);
+         w = __ieee754_j1l (x);
+         if (fabsl (z) >= fabsl (w))
+           b = (t * z / b);
+         else
+           b = (t * w / a);
        }
     }
   if (sgn == 1)

diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c
index 3d715d36aada88f3e0feebf79da0f21e3eb95ffa..bedff7d566b161c7f2de795f158d7f692f43021a 100644 (file)
--- a/sysdeps/ieee754/ldbl-96/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-96/e_jnl.c
@@ -281,7 +281,16 @@ __ieee754_jnl (n, x)
                    }
                }
            }
-         b = (t * __ieee754_j0l (x) / b);
+         /* j0() and j1() suffer enormous loss of precision at and
+          * near zero; however, we know that their zero points never
+          * coincide, so just choose the one further away from zero.
+          */
+         z = __ieee754_j0l (x);
+         w = __ieee754_j1l (x);
+         if (fabsl (z) >= fabsl (w))
+           b = (t * z / b);
+         else
+           b = (t * w / a);
        }
     }
   if (sgn == 1)