faster besselj0 for x<25 #15

oscardssmith · 2022-02-28T03:41:01Z

this is about 2x faster on my computer where applicable. This also should be done for besselj1 and bessely0, bessely1

oscardssmith · 2022-02-28T03:42:56Z

Here's a plot of relative errors.

heltonmc · 2022-02-28T03:47:17Z

Though it's another branch. It might be best to keep the power series for very small arguments to ensure that besselj0(0.0)= 1.0 (test that fails currently).

oscardssmith · 2022-02-28T03:51:08Z

Fixed. I used 1 to few terms, and didn't force the 1st coefficient to be 1.0.

codecov-commenter · 2022-02-28T03:52:37Z

Codecov Report

Merging #15 (927c1a5) into master (f772ccf) will decrease coverage by 0.00%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #15      +/-   ##
==========================================
- Coverage   99.56%   99.55%   -0.01%     
==========================================
  Files          12       13       +1     
  Lines         683      676       -7     
==========================================
- Hits          680      673       -7     
  Misses          3        3

Impacted Files	Coverage Δ
src/Bessels.jl	`100.00% <ø> (ø)`
src/Polynomials/besselj_polys.jl	`100.00% <100.00%> (ø)`
src/besselj.jl	`100.00% <100.00%> (ø)`
src/bessely.jl	`100.00% <100.00%> (ø)`
src/constants.jl	`100.00% <100.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update f772ccf...927c1a5. Read the comment docs.

oscardssmith · 2022-02-28T03:56:34Z

Also, for future reference, here is the code used to generate the polynomials

using ArbNumerics
function r(x,root)
    abs(x)<eps(Float64) && return -besselj(1,ArbFloat(x)+root[1]+root[2])
    return besselj(0,ArbFloat(x)+root[1]+root[2])/x
end
function r2(x,root)
    return besselj(0,ArbFloat(x)+root[1]+root[2])
end
polys = []
for i in 2:2:16
    poly = Float64.(Tuple(ratfn_minimax(x->r(x,pts[i]), ((pts[i-1][1]-pts[i][1]-.1)/2, (pts[i+1][1]-pts[i][1]+.1)/2), 12, 0)[1]))
    push!(polys, (0.0, poly...))
    poly = Float64.(Tuple(ratfn_minimax(x->r2(x,pts[i+1]), ((pts[i][1]-pts[i+1][1]-.1)/2, (pts[i+2][1]-pts[i+1][1]+.1)/2), 13, 0)[1]))
    push!(polys, poly)
end
Tuple(polys)

for the 0-pi/2 range, the invocation was (1.0, Float64.(Tuple(ratfn_minimax(x->(besselj(0,sqrt(x))-1)/x, (1e-20,nextfloat(pi/2)^2), 7, 0)[1]))...)

heltonmc · 2022-02-28T04:00:24Z

We may also want to extend this range just a bit more to blend with the asymptotic expansion better. I had chosen that cutoff to blend well with the previous method but we may want to use this for x < 26.0 now ?

oscardssmith · 2022-02-28T04:02:04Z

Makes sense. Do you want to do the honors, or should I?

heltonmc · 2022-02-28T04:07:50Z

If you have it set up that might be easiest! Otherwise, I can look more at this tomorrow and also for J1, Y0, and Y1. But this is awesome. Hoping to finally finish these special cases! (at least for Float64)

oscardssmith · 2022-02-28T04:34:13Z

updated (the post with computation details and error plots are also updated).

heltonmc · 2022-02-28T16:28:15Z

To clarify the above code. You will need pts defined. This code runs for me and produces similar values. Though, I think the last point used is different? Also, pts contains the zeros and extrema of the bessel functions. What is the first point you used?

using ArbNumerics, Remez
pts = (
    (0.73676709589402, 0.0),
    (2.404825557695773, -1.176691651530894e-16),
    (3.8317059702075125, -1.5269184090088067e-16),
    (5.520078110286311,  8.088597146146722e-17),
    (7.015586669815619,  -9.414165653410389e-17),
    (8.653727912911013, -2.92812607320779e-16 ),
    (10.173468135062722,   4.482162274768888e-16),
    (11.791534439014281,  2.812956912778735e-16),
    (13.323691936314223,   2.600408064718813e-16),
    (14.930917708487787, -7.070514505983074e-16),
    (16.470630050877634,  -1.619019544798128e-15),
    (18.071063967910924, -9.658048089426209e-16),
    (19.615858510468243,  -1.004445634526616e-15),
    (21.21163662987926 ,  4.947077428784068e-16),
    (22.760084380592772,  -4.925749373614922e-16),
    (24.352471530749302,  9.169067133951066e-16),
    (25.903672087618382, 4.894530726419825e-16),
    (26.1,0.0)
)
function r(x,root)
    abs(x)<eps(Float64) && return -besselj(1,ArbFloat(x)+root[1]+root[2])
    return besselj(0,ArbFloat(x)+root[1]+root[2])/x
end
function r2(x,root)
    return besselj(0,ArbFloat(x)+root[1]+root[2])
end
polys = []
for i in 2:2:16
    poly = Float64.(Tuple(ratfn_minimax(x->r(x,pts[i]), ((pts[i-1][1]-pts[i][1]-.1)/2, (pts[i+1][1]-pts[i][1]+.1)/2), 12, 0)[1]))
    push!(polys, (0.0, poly...))
    poly = Float64.(Tuple(ratfn_minimax(x->r2(x,pts[i+1]), ((pts[i][1]-pts[i+1][1]-.1)/2, (pts[i+2][1]-pts[i+1][1]+.1)/2), 13, 0)[1]))
    push!(polys, poly)
end
Tuple(polys)

Edited with correct values.

oscardssmith · 2022-02-28T16:56:50Z

The last 2 points are wrong. They should be
(25.903672087618382, 4.894530726419825e-16) and (26.1,0.0). (you want the last point to be the reflection of the last root around your cutoff value).

heltonmc · 2022-02-28T17:16:17Z

It looks like I also misspoke were pts contains the roots of besselj0 and besselj1? Unsure what the first point is though.

oscardssmith · 2022-02-28T17:41:44Z

The derivative of besselj0 is -besselj1, so extrema of besselj0 are roots of besselj1. The first point is a reflection of the second point around x=pi/2 (which is the lower cutoff)

oscardssmith · 2022-03-01T05:20:50Z

This now is implemented for besselj1 although something with my polynomials is slightly broken, which is causing spikes in error at the extrema of the function which really shouldn't be happening.

oscardssmith · 2022-03-01T16:18:57Z

This is now failing tests because SpecialFunctions is too inaccurate near 3.8317059702075125. The tests probably need to use Arb.

heltonmc · 2022-03-01T18:15:01Z

I agree that we should probably switch the tests over to Arb for all this. I actually had a draft of this early on but ran into JeffreySarnoff/ArbNumerics.jl#48 causing some issues. Now that it is fixed we should move over.
Though, for this test besselj1 is computed in higher precision I believe using SpecialFunctions.jl as they provide BigFloat implementations. I'll have to check this some more.

src/besselj.jl

heltonmc · 2022-03-01T21:28:17Z

src/constants.jl

-    2.56987256757748830383E5, 6.20836478118054335476E2,
-    1.00000000000000000000E0
-)
-const PP_j1(::Type{Float64}) = (


Unfortunately several of these constants are used for the bessely implementations.

it will work out once we replace those implementations as well :)

O most definitely!! Didn't know if you wanted to do that in separate PR or all at once 😃

I'm probably going to save those for a separate PR. They look worse.

src/bessely.jl

heltonmc · 2022-03-02T04:45:49Z

src/besselj.jl

+        (24.352471530749302, 9.169067133951066e-16),
+        (25.903672087618382, 4.894530726419825e-16)
+)
+@inline J0_POLYS(::Type{Float64}) = (


Is there a reason we sometimes inline these, sometimes we don't, and sometimes we use constants for the polynomials? It might be good to have a defined way going forward that is consistent. Though this one is slightly different as we have to determine which polynomial at runtime. Does this affect compile time?

Also, I would be a fan of moving the polynomials to their own file and slightly different file structure to make the core code easier to read. I think this approach could also be used for Double64 implementations that would also take up a lot of space in main file.

heltonmc · 2022-03-02T04:49:04Z

src/besselj.jl

@@ -114,34 +217,29 @@ function besselj1(x::Float64)
    x = abs(x)
    isinf(x) && return zero(x)


We should probably be consistent with the structure between besselj0 and besselj1. Either have the isinf check in the large argument branch (where I check if xinv is zero) or just have it at beginning.

heltonmc · 2022-03-02T04:52:54Z

src/besselj.jl

-        p = p * sc[2] - w * q * sc[1]
-        return p * SQ2OPI(T) / sqrt(x)
+    if x <= 26.0
+        x <= pi/2 && return x*evalpoly(x*x, J1_POLY_PIO2(T))


Also, do we want to set a code style for this repo? I think in general I've been trying to use blue but sometimes I obsess over spaces too much 😅

I tend to write with relatively little spacing consistency. If you have a style guide you would prefer, I can adhere to it, but my default tends to be somewhat arbitrary.

I usually try to follow https://github.com/jrevels/YASGuide. I can make some spacing changes if that's ok.

Otherwise I think this is ready to merge?

sounds good. (make any changes you want).

heltonmc · 2022-03-02T04:55:09Z

Tests look good now. Do we want to also update the Float32 implementations here or later?

oscardssmith · 2022-03-02T05:06:23Z

I think we should do Float32 separately because I think this method might not be the best for lower precision.

heltonmc · 2022-03-04T00:20:38Z

I added a fix to the sign error for besselj1 and updated spacing, moved polys to own file. The sign error could be done in a different way perhaps. This is ready to merge on my end.

oscardssmith · 2022-03-04T00:21:43Z

src/besselj.jl

    else
        xinv = inv(x)
-        x2 = xinv*xinv
+        iszero(xinv) && return zero(T)


technically this should probably return xinv so it returns -0.0 when applicable.

I believe since we have x=abs(x) this would still always return 0.0 so might need to handle that separately. I'm actually starting to not appreciate that syntax of overwriting that variable and prefering xabs = abs(x). Thoughts?

I don't think it really matters. I would probably just merge this as is.

faster besselj0 for x<25

4fe2d3b

fix behavior at 0

73179d2

oscardssmith added 2 commits February 27, 2022 23:24

extend to 26 (gives better matching)

51cb088

update docs

809e51b

oscardssmith and others added 2 commits February 28, 2022 09:04

remove unused coefficients.

a084a62

fix type stability and spacing

c9c81d9

update besselj1

f349610

oscardssmith added 7 commits March 1, 2022 00:58

comment out generation

2cc0361

fix coeficients

022db33

fix typo

4d10243

remove unused j1 polys.

f06cd02

typo

9c09a19

Update besselj.jl

3bab7d7

Update besselj.jl

1213631

heltonmc reviewed Mar 1, 2022

View reviewed changes

src/besselj.jl Outdated Show resolved Hide resolved

oops

6a5ccc5

heltonmc reviewed Mar 1, 2022

View reviewed changes

re-add constants

0a8f77a

heltonmc reviewed Mar 2, 2022

View reviewed changes

src/bessely.jl Outdated Show resolved Hide resolved

fix typo

d3e7deb

heltonmc reviewed Mar 2, 2022

View reviewed changes

add neg tests, fix spacing

927c1a5

oscardssmith commented Mar 4, 2022

View reviewed changes

heltonmc merged commit ded7c2d into JuliaMath:master Mar 4, 2022

oscardssmith deleted the faster-besselj0s branch March 4, 2022 00:47

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

faster besselj0 for x<25 #15

faster besselj0 for x<25 #15

oscardssmith commented Feb 28, 2022

oscardssmith commented Feb 28, 2022 •

edited

Loading

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

codecov-commenter commented Feb 28, 2022 •

edited

Loading

oscardssmith commented Feb 28, 2022 •

edited

Loading

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

heltonmc commented Feb 28, 2022 •

edited

Loading

oscardssmith commented Feb 28, 2022

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

oscardssmith commented Mar 1, 2022

oscardssmith commented Mar 1, 2022

heltonmc commented Mar 1, 2022

heltonmc Mar 1, 2022

oscardssmith Mar 1, 2022

heltonmc Mar 1, 2022

oscardssmith Mar 2, 2022

heltonmc Mar 2, 2022

heltonmc Mar 2, 2022

heltonmc Mar 2, 2022

heltonmc Mar 2, 2022

oscardssmith Mar 2, 2022

heltonmc Mar 3, 2022

oscardssmith Mar 3, 2022

heltonmc commented Mar 2, 2022

oscardssmith commented Mar 2, 2022

heltonmc commented Mar 4, 2022 •

edited

Loading

oscardssmith Mar 4, 2022

heltonmc Mar 4, 2022

oscardssmith Mar 4, 2022

		@@ -114,34 +217,29 @@ function besselj1(x::Float64)
		x = abs(x)
		isinf(x) && return zero(x)

faster besselj0 for x<25 #15

faster besselj0 for x<25 #15

Conversation

oscardssmith commented Feb 28, 2022

oscardssmith commented Feb 28, 2022 • edited Loading

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

codecov-commenter commented Feb 28, 2022 • edited Loading

Codecov Report

oscardssmith commented Feb 28, 2022 • edited Loading

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

heltonmc commented Feb 28, 2022 • edited Loading

oscardssmith commented Feb 28, 2022

heltonmc commented Feb 28, 2022

oscardssmith commented Feb 28, 2022

oscardssmith commented Mar 1, 2022

oscardssmith commented Mar 1, 2022

heltonmc commented Mar 1, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

heltonmc commented Mar 2, 2022

oscardssmith commented Mar 2, 2022

heltonmc commented Mar 4, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

oscardssmith commented Feb 28, 2022 •

edited

Loading

codecov-commenter commented Feb 28, 2022 •

edited

Loading

oscardssmith commented Feb 28, 2022 •

edited

Loading

heltonmc commented Feb 28, 2022 •

edited

Loading

heltonmc commented Mar 4, 2022 •

edited

Loading