bug in j0f()

[Bruce Evans]

On Tue, 2 Dec 2014, Steve Kargl wrote:

> On Tue, Dec 02, 2014 at 04:29:08PM -0800, Steve Kargl wrote:
>> On Tue, Dec 02, 2014 at 04:09:41PM -0800, Steve Kargl wrote:
>>> On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote:
>>>> Anyone object to the following patch?

OK (with 0x54000000).

>>>> Index: e_j0f.c
>>>> ===================================================================
>>>> --- e_j0f.c	(revision 275211)
>>>> +++ e_j0f.c	(working copy)
>>>> @@ -62,7 +62,7 @@
>>>>  	 * 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)/sqrtf(x);
>>>> +		if(ix>0x4b800000) z = (invsqrtpi*cc)/sqrtf(x);
>>>
>>> Exhaustive testing in the range 0x1p38 to 0x1p100
>>> indicated at the constant should be 0x54000000.

My tests agree.  Tested on amd64 and i386.

>> Note, a similar wrong condition exists within y0f().  I have
>> not tested y0f(), but propose making a similar change in y0f()
>> as well.

Not so exhaustive testing gave 0x54800000 on amd64.

> While I'm monologuing, I might as well point out that the
> rational approximations in j0f (and y0f and most likely
> j1f and y1f) are over-specified.  I suspect that the
> polynomials in the rational approximation can be reduced
> by one or two terms.

Also, the cutoffs of 2**-13 and 2**-27 are the same for both precisions,
thus likely to be wrong for float precision.

Bruce