svn commit: r215237 - head/lib/msun/src

[Alexander Best]

On Sat Nov 13 10, Ulrich Spoerlein wrote:
> Author: uqs
> Date: Sat Nov 13 10:54:10 2010
> New Revision: 215237
> URL: http://svn.freebsd.org/changeset/base/215237
> 
> Log:
>   Fix bug in jn(3) and jnf(3) that led to -inf results

thank you very much for fixing this long outstanding issue.
you might want to have a look at [1], where two more issues have been reported.

cheers.
alex

[1] http://mailman.oakapple.net/pipermail/numeric-interest/2010-September/thread.html

>   
>   Explanation by Steve:
>   jn[f](n,x) for certain ranges of x uses downward recursion to compute
>   the value of the function.  The recursion sequence that is generated is
>   proportional to the actual desired value, so a normalization step is
>   taken.  This normalization is j0[f](x) divided by the zeroth sequence
>   member.  As Bruce notes, near the zeros of j0[f](x) the computed value
>   can have giga-ULP inaccuracy. I found for the 1st zero of j0f(x) only
>   the leading decimal digit is correct.  The solution to the issue is
>   fairly straight forward.  The zeros of j0(x) and j1(x) never coincide,
>   so as j0(x) approaches a zero, the normalization constant switches to
>   j1[f](x) divided by the 2nd sequence member.  The expectation is that
>   j1[f](x) is a more accurately computed value.
>   
>   PR:		bin/144306
>   Submitted by:	Steven G. Kargl <kargl at troutmask.apl.washington.edu>
>   Reviewed by:	bde
>   MFC after:	7 days
> 
> Modified:
>   head/lib/msun/src/e_jn.c
>   head/lib/msun/src/e_jnf.c
> 
> Modified: head/lib/msun/src/e_jn.c
> ==============================================================================
> --- head/lib/msun/src/e_jn.c	Sat Nov 13 10:38:06 2010	(r215236)
> +++ head/lib/msun/src/e_jn.c	Sat Nov 13 10:54:10 2010	(r215237)
> @@ -200,7 +200,12 @@ __ieee754_jn(int n, double x)
>  			}
>  	     	    }
>  		}
> -	    	b = (t*__ieee754_j0(x)/b);
> +		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;
> 
> Modified: head/lib/msun/src/e_jnf.c
> ==============================================================================
> --- head/lib/msun/src/e_jnf.c	Sat Nov 13 10:38:06 2010	(r215236)
> +++ head/lib/msun/src/e_jnf.c	Sat Nov 13 10:54:10 2010	(r215237)
> @@ -152,7 +152,12 @@ __ieee754_jnf(int n, float x)
>  			}
>  	     	    }
>  		}
> -	    	b = (t*__ieee754_j0f(x)/b);
> +		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;

-- 
a13x