[newlib-cygwin] Fix spurious underflow exceptions for Bessel functions for double(from glibc bug 14155)

[newlib-cygwin] Fix spurious underflow exceptions for Bessel functions for double(from glibc bug 14155)

[Corinna Vinschen]

https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;h=5e24839658f6576b68b26c977897b9ad3fc3c23f

commit 5e24839658f6576b68b26c977897b9ad3fc3c23f
Author: Joseph S. Myers <joseph@codesourcery.com>
Date:   Wed Mar 25 11:18:44 2020 -0700

    Fix spurious underflow exceptions for Bessel functions for double(from glibc bug 14155)
    
            This fix comes from glibc, from files which originated from
            the same place as the newlib files. Those files in glibc carry
            the same license as the newlib files.
    
    Bug 14155 is spurious underflow exceptions from Bessel functions for
    large arguments.  (The correct results for large x are roughly
    constant * sin or cos (x + constant) / sqrt (x), so no underflow
    exceptions should occur based on the final result.)
    
    There are various places underflows may occur in the intermediate
    calculations that cause the failures listed in that bug.  This patch
    fixes problems for the double version where underflows occur in
    calculating the intermediate functions P and Q (in particular, x**-12
    gets computed while calculating Q).  Appropriate approximations are
    used for P and Q for arguments at least 0x1p28 and above to avoid the
    underflows.
    
    For sufficiently large x - 0x1p129 and above - the code already has a
    cut-off to avoid calculating P and Q at all, which means the
    approximations -0.125 / x and 0.375 / x can't themselves cause
    underflows calculating Q.  This cut-off is heuristically reasonable
    for the point beyond which Q can be neglected (based on expecting
    around 0x1p-64 to be the least absolute value of sin or cos for large
    arguments representable in double).
    
    The float versions use a cut-off 0x1p17, which is less heuristically
    justifiable but should still only affect values near zeroes of the
    Bessel functions where these implementations are intrinsically
    inaccurate anyway (bugs 14469-14472), and should serve to avoid
    underflows (the float underflow for jn in bug 14155 probably comes
    from the recurrence to compute jn).  ldbl-96 uses 0x1p129, which may
    not really be enough heuristically (0x1p143 or so might be safer - 143
    = 64 + 79, number of mantissa bits plus total number of significant
    bits in representation) but again should avoid underflows and only
    affect values where the code is substantially inaccurate anyway.
    ldbl-128 and ldbl-128ibm share a completely different implementation
    with no such cut-off, which I propose to fix separately.
    
    Signed-off-by: Keith Packard <keithp@keithp.com>

Diff:
---
 newlib/libm/math/e_j0.c  | 6 ++++--
 newlib/libm/math/e_j1.c  | 6 ++++--
 newlib/libm/math/ef_j0.c | 6 +++---
 newlib/libm/math/ef_j1.c | 4 ++--
 4 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/newlib/libm/math/e_j0.c b/newlib/libm/math/e_j0.c
index 13773cbf9..d3af9d32c 100644
--- a/newlib/libm/math/e_j0.c
+++ b/newlib/libm/math/e_j0.c
@@ -338,7 +338,8 @@ static double pS2[5] = {
 	__int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = pR8; q= pS8;}
+	if (ix>=0x41b00000)    {return one;}
+	else if(ix>=0x40200000){p = pR8; q= pS8;}
 	else if(ix>=0x40122E8B){p = pR5; q= pS5;}
 	else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
       else {p = pR2; q= pS2;}
@@ -474,7 +475,8 @@ static double qS2[6] = {
 	__int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qR8; q= qS8;}
+	if (ix>=0x41b00000)    {return -.125/x;}
+	else if(ix>=0x40200000){p = qR8; q= qS8;}
 	else if(ix>=0x40122E8B){p = qR5; q= qS5;}
 	else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
       else {p = qR2; q= qS2;}
diff --git a/newlib/libm/math/e_j1.c b/newlib/libm/math/e_j1.c
index 098eb569e..72855e3fa 100644
--- a/newlib/libm/math/e_j1.c
+++ b/newlib/libm/math/e_j1.c
@@ -336,7 +336,8 @@ static double ps2[5] = {
         __int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-        if(ix>=0x40200000)     {p = pr8; q= ps8;}
+	if (ix>=0x41b00000)    {return one;}
+	else if(ix>=0x40200000){p = pr8; q= ps8;}
         else if(ix>=0x40122E8B){p = pr5; q= ps5;}
         else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
         else {p = pr2; q= ps2;}
@@ -473,7 +474,8 @@ static double qs2[6] = {
 	__int32_t ix;
 	GET_HIGH_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qr8; q= qs8;}
+	if (ix>=0x41b00000)    {return .375/x;}
+	else if(ix>=0x40200000){p = qr8; q= qs8;}
 	else if(ix>=0x40122E8B){p = qr5; q= qs5;}
 	else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
       else {p = qr2; q= qs2;}
diff --git a/newlib/libm/math/ef_j0.c b/newlib/libm/math/ef_j0.c
index 866cfcf96..854801f1d 100644
--- a/newlib/libm/math/ef_j0.c
+++ b/newlib/libm/math/ef_j0.c
@@ -74,7 +74,7 @@ static float zero = 0.0;
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x80000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);
+		if(ix>0x5c000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);
 		else {
 		    u = pzerof(x); v = qzerof(x);
 		    z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(x);
@@ -156,14 +156,14 @@ v04  =  4.4111031494e-10; /* 0x2ff280c2 */
                     if ((s*c)<zero) cc = z/ss;
                     else            ss = z/cc;
                 }
-                if(ix>0x80000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
+                if(ix>0x5c000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
                 else {
                     u = pzerof(x); v = qzerof(x);
                     z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x);
                 }
                 return z;
 	}
-	if(ix<=0x32000000) {	/* x < 2**-27 */
+	if(ix<=0x39800000) {	/* x < 2**-27 */
 	    return(u00 + tpi*__ieee754_logf(x));
 	}
 	z = x*x;
diff --git a/newlib/libm/math/ef_j1.c b/newlib/libm/math/ef_j1.c
index 01bd24cf1..f4c9c9dd3 100644
--- a/newlib/libm/math/ef_j1.c
+++ b/newlib/libm/math/ef_j1.c
@@ -75,7 +75,7 @@ static float zero    = 0.0;
 	 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 	 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x80000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(y);
+		if(ix>0x5c000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(y);
 		else {
 		    u = ponef(y); v = qonef(y);
 		    z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(y);
@@ -153,7 +153,7 @@ static float V0[5] = {
          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
          * to compute the worse one.
          */
-                if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
+                if(ix>0x5c000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
                 else {
                     u = ponef(x); v = qonef(x);
                     z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x);