A different approach to compute besselk in intermediate zones

[heltonmc]

@cgeoga I was thinking some more about how we compute the region between asymptotic expansions (power series, uniform expansions, and large argument). Which corresponds to roughly 0<nu<25 and 2 < x < 30...

Right now, we are using the wronskian relation to besseli. The reason we can do this is that the power series of besseli is convergent in this zone so using that along with the continued fraction and wronskian we can get besselk pretty accurately. Now there are a couple drawbacks to that:

1. We need besseli(nu, x) and besseli(nu + 1, x) and then also need to compute besselk(nu+1,x)/besselk(nu,x).
2. Power series are slower to converge as x gets bigger
3. The power series of besseli(nu, x) suffers from cancellation for complex numbers so can't be as uniformly used.

The first point was taken care of by using SIMD instructions to compute both besseli(nu, x) and besseli(nu + 1, x). So that is pretty optimized (though we pay an extra penalty for complex numbers). But points 2 and 3 are pretty much unsolvable. Power series will always be slow for larger arguments and can't use much for complex numbers.

I was toying with an alternative which links a few ideas. Now the continued fraction approach to compute besselk(nu+1,x)/besselk(nu,x)converges more quickly asxgets bigger so it has an inverse relation to the power series. (i.e., power series converges faster for smallxwhile continued fraction converges faster for bigx`). So what if there was a good way to use the uniform expansions to link with the smaller orders. Now Miller's algorithm links stable backward recurrence to higher orders with a normalization factor calculated from low orders. That is fine but usually we do not know the lower orders unless for integer arguments.

I'm sure this idea has been described before but can't seem to find a description of it. But there appears to be a way to link stable forward recurrence with a normalization factor calculated at higher orders to compute besselk at lower orders. It works by instead of considering trial start values (0 and epsilon) in the Miller backwards scheme, it uses the continued fraction approach to generate start values in the forward scheme. So the algorithm would look something like:

1. Generate start values 1.0 and knup1/knu.
2. Use upward recurrence up until a order you can use asymptotic expansion (nu_ceil).
3. Compute besselk(nu_ceil, x) using asymptotic expansion
4. Normalize starting value

It may seem complicated but instead of having to compute besseli(nu, x), besseli(nu + 1, x) and besselk(nu+1,x)/besselk(nu,x)we just needbesselk(nu_ceil, x)andbesselk(nu+1,x)/besselk(nu,x). There is a few extra steps of course to generate the recurrence but that is pretty small.

So now we have an algorithm that computes besselk(nu_ceil, x) at constant time using uniform asymptotic expansion and algorithm to compute besselk(nu+1,x)/besselk(nu,x)which converges faster asx` increases. A sample implementation of the algorithm would look like..

using Bessels

function test(nu, x)
    fpart, _ = modf(nu)
    nu_ceil = 25.0 + fpart
    knu_ceil = Bessels.besselk_large_orders(nu_ceil, x)
    knup1 = Bessels.besselk_ratio_knu_knup1(nu, x)
    return knu_ceil / Bessels.besselk_up_recurrence(x, knup1, 1.0, nu + 1, nu_ceil)[1]
end

Speed tests to other approach.

## high precision calculation
julia> 2.308550810254649099632E-14
2.308550810254649e-14

julia> test(2.2, 30.0)
2.3085508102546522e-14

julia> Bessels.besselk_continued_fraction(2.2, 30.0)
2.3085508102546462e-14

julia> @benchmark test(2.2, x) setup = (x = 30.0 + rand())
BenchmarkTools.Trial: 10000 samples with 968 evaluations.
 Range (min … max):  79.846 ns … 120.610 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     79.976 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   80.832 ns ±   2.255 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆█▃                    ▂▃           ▃▅▁                      ▁
  ███▆▅▅▄▃▄▄▄▄▆▄▄▃▁▁▄▃▄▁▆██▆▅▃▃▃▄▄▄▅▄▆███▅▄▅▃▃▁▃▁▄▃▁▄▃▃▁▃▃▁▁▃█ █
  79.8 ns       Histogram: log(frequency) by time      87.2 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Bessels.besselk_continued_fraction(2.2, x) setup = (x = 30.0 + rand())
BenchmarkTools.Trial: 10000 samples with 898 evaluations.
 Range (min … max):  127.180 ns … 178.360 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     129.315 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   128.840 ns ±   1.969 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   █▇▂              █▇                                           
  ▃███▃▂▂▂▂▂▂▂▂▁▂▂▂███▄▃▂▂▂▂▂▃▃▃▄▃▂▂▂▂▂▂▁▂▂▂▂▃▃▂▂▂▂▂▁▁▂▂▂▂▂▂▂▂▂ ▃
  127 ns           Histogram: frequency by time          134 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Which is about 40% faster for bigger arguments. For smallish arguments its about 10% faster.

julia> @benchmark Bessels.besselk_continued_fraction(2.2, x) setup = (x = 10.0 + rand())
BenchmarkTools.Trial: 10000 samples with 940 evaluations.
 Range (min … max):  102.305 ns … 145.612 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     104.122 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   104.488 ns ±   2.483 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅▆     ▄█▅▂   ▁    ▅▄▃▁     ▁                                 ▂
  ██▇▁▁▁▁████▆▄▄█▇▅▅▄████▇▅▅▄▇█▇▆▆▇▇▅▄▁▃▄▅▅▆▅▆▅▆▆▅▄▅▄▄▃▄▄▅▅▃▄▄▄ █
  102 ns        Histogram: log(frequency) by time        115 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark test(2.2, x) setup = (x = 10.0 + rand())
BenchmarkTools.Trial: 10000 samples with 954 evaluations.
 Range (min … max):  93.203 ns … 134.565 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     97.615 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   96.402 ns ±   2.687 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇▅              ▁       █▅              ▁         ▁▃         ▂
  ██▇▅▄▅▄▅▅▄▄▁▁▃▄▅█▆▅▄▁▁▁███▇▅▅▅▅▇▆▆▅▁▄▁▃▇█▇▆▅▅▄▃▅▁▅██▆▃▄▁▄▁▄▃ █
  93.2 ns       Histogram: log(frequency) by time       104 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

And for much smaller arguments it's hard to beat the power series so it is slower.

julia> @benchmark test(2.2, x) setup = (x = 4.0 + rand())
BenchmarkTools.Trial: 10000 samples with 918 evaluations.
 Range (min … max):  115.423 ns … 157.588 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     121.642 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   122.518 ns ±   4.119 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆              █             ▂█              ▆              ▃ ▂
  █▇▃▄▄▆▃▃▁▁▃▁▁▁▇██▅▃▆█▆▄▄▄▃▃▄▁███▄▁▅█▆▄▃▄▁▃▃▄▅██▃▄▄▅▅▃▁▃▁▄▁▃▁█ █
  115 ns        Histogram: log(frequency) by time        133 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Bessels.besselk_continued_fraction(2.2, x) setup = (x = 4.0 + rand())
BenchmarkTools.Trial: 10000 samples with 942 evaluations.
 Range (min … max):  100.849 ns … 158.882 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     108.678 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   108.736 ns ±   4.641 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄    ▅      █▂  ▅▂  ▂▆  ▇▄   ▃  ▂▅  ▅▃   ▁   ▃   ▃            ▂
  █▄▁▄▇█▅▁▆▃▃▄██▆▆██▄▄███▅██▆▆██▆▄██▇▆██▅▆▅█▆▅██▇▅▆█▅▅▄▇▄▄▄██▁█ █
  101 ns        Histogram: log(frequency) by time        123 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

So there is of course things to consider but this approach has the advantage of being more applicable to a large parameter set including complex numbers and has the potential for large speed gains over a fairly big domain. My initial implementation could also probably be improved to speed things up a bit but it does serve as a good way to link directly to the asymptotic expansions.

[cgeoga]

@heltonmc, apologies---I somehow completely missed this when you posted it and am just seeing it now for the first time. I like this idea a lot and tried to think about some similar ideas to leverage the large-order expansions for these intermediate argument ranges, but I never really got anywhere. This is very nice. I didn't realize that the power series cutoff was 2---for some reason I thought it was 4. But even by 3 it seems to draw with the power series, and so to me it seems like a pretty obvious win to give up ~10ns for 2 < x < 3 but then gain an increasingly substantial amount when 3 < x < 30.

I think I honestly need a little more time to grok why this works at all, though. I never really understood the Miller recursion in a very deep level and this whole concept of normalization is still a little shaky to me. But I'll definitely play with this and try to understand it, and it seems like a seriously exciting option.

[heltonmc]

Not a problem!! I was working on this before I needed to switch to dissertation stuff so haven't come back to it. Hopefully will work on it a little more in the next couple of weeks.

Yes - the power series for besseli is pretty optimized and rapidly convergent for small arguments so it is hard to beat for x<3. But I was also thinking about the idea of using this approach for bessely because you can't use the besselj power series the same way you can for besseli. So it is is much more applicable to that function but also here when the besseli is slowly converging (larger arguments) or prone to severe cancellation (complex numbers).

I also like the idea of reducing the number of algorithms employed such that if the uniform expansion is the main driver of the speed and accuracy of almost all arguments that is more easily optimized than several different routines. But ya, the Miller algorithm seems to be a pretty common way to generate higher order (integer orders) functions when forward recurrence is not stable. That way you can use an optimized zero order function to normalize for any order. This doesn't work as well for besselk and bessely though when forward recurrence is stable..

Anyway, still need to have the continued fraction approach for bessely..... that could reduce significantly the amount of code and compilation costs (though that seems to be less of an issue now)

[heltonmc]

Was thinking about this approach for bessely again.... it might be possible to derive the continued fraction for bessely but I think its also possible to use the relation between besselk and the Hankel function. Then use the relation of the Hankel function to get bessely from besselk therefore we can just use the same approach described above using besselk_ratio_knu_knup1 with some modifications.

The first challenge is that it requires a rotation of argument by pi/2 so real arguments for bessely will require evaluation of besselk along imaginary axis. That's not too big of a deal for the continued fraction approach other than slowing it down slightly. I made some adjustments to the existing approach here to speed it up slightly for complex arguments....

function besselk_ratio_knu_knup1(v, x::ComplexOrReal{T}) where T
    MaxIter = 1000
    S = eltype(x)
    (hn, Dn, Cn) = (S(1e-50), zero(S), S(1e-50))

    jf = one(T)
    vv = v * v
    @fastmath begin
        for _ in 1:MaxIter
            an = vv - (jf - 0.5)^2
            bn = 2 * (x + jf)
            Cn = an / Cn + bn
            Dn = inv(muladd(an, Dn, bn))
            del = Dn * Cn
            hn *= del
            checktol(del - 1) && break
            jf += one(T)
        end
        xinv = (v + x + hn + 1/2) / x
    end 
    return xinv
end
@inline checktol(err::T) where T = abs(err) < eps(Float64)
@inline checktol(err::Complex{T}) where T = abs2(err) < eps(T)^2

Perhaps, it could be sped up further but it is pretty fast already for the input range we need for bessely (7<x<20)...

julia> @btime besselk_ratio_knu_knup1(v, x) setup=(x= 20.0im; v = rand()*1)
  76.246 ns (0 allocations: 0 bytes)
1.0000915788150249 - 0.045999085366594im

julia> @btime besselk_ratio_knu_knup1(v, x) setup=(x= 7.0im; v = rand()*1)
  141.996 ns (0 allocations: 0 bytes)
1.001602364037009 - 0.11366618431295909im

So a range between 76-142 nanoseconds in the regime we'd utilize this for even with complex arithmetic... Of course we also need recurrence up to higher order expansion (takes around 20 nanoseconds) and then need one large order evaluation...

julia> @btime Bessels.besseljy_debye(y, x) setup=(x=10*rand()+12.0; y=rand()+60.0)
  49.393 ns (0 allocations: 0 bytes)
(1.8133319226273097e-26, -3.0253831328694315e23)

Which takes around 50 nanoseconds. So that puts the total time for this approach around 150-225 nanoseconds in a region where no expansion gives good results. Now that might not seem that fast but the existing approach using Chebyshev approximation is the longest approach we have...

julia> @btime Bessels.bessely_chebyshev(y, x) setup=(x=12*rand()+6; y=rand()*20.0)
  460.577 ns (0 allocations: 0 bytes)

Which takes between 460-500 nanoseconds so the new approach would be about 2-3x faster.

I think it might be worth taking a look at the whole routine to some degree and try to decrease the amount of algorithms to improve branch prediction.

However, there was some discussion on discourse about this but it might be possible to improve the Chebyshev evaluation using SIMD.... I'd assume it would be possible to improve those times by 2-3x

Edit

The relation between besselk and bessely through the Hankel function might require a bit more care. It's true that besselk(v+1, z) / besselk(v, z) = im * besselh1(v + 1, z) / besselh1(v, z) relating that easily to bessely(v+1, z) / bessely(v, z) isn't as straight forward without computing besselj.

So we could compute besselh1 this way as we have the higher order expansion to compute besselh1 with large orders. The challenge though is that forward recurrence isn't stable in the real part of besselh1 whereas the imaginary part is. I believe we would actually need both parts for the normalization condition so this wouldn't necessarily work...

Using the relation bessely(v, x) = -cispi(-v/2) * (conj(besselk(v, x*im)) + cispi(v) * besselk(v, im*x)) / pi might be helpful here if we can calculate besselk. So maybe it is worth starting with computing besselk in the complex plane... The issue is that the along the imaginary axis the debye expansion doesn't converge well. Will need to think more about that.

Edit 2

After some more consideration it might be best to compute Hankel function directly then just take bessely from that.
We can use the relation besselk(v+1, x) / besselk(v, x) = im * hankelh1(v+1, z * im) / hankelh1(v, z * im) and use the continued fraction approach for besselk.

Similar to what I showed above we can construct the solution to the Hankel function with the continued fraction and the solution for higher order term.

# computes Hankel function
# compute bessely from imaginary part
function test(nu, x)
    fpart, _ = modf(nu)
    nu_ceil = 45.0 + fpart
    knu_ceil = complex(Bessels.besseljy_debye(nu_ceil, x)...)
    knup1 = -im * conj(Bessels.besselk_ratio_knu_knup1(nu, x*im))
    return knu_ceil / Bessels.besselj_up_recurrence(x, knup1, one(knup1), nu + 1, nu_ceil)[1]
end

julia> test(2.5, 18.0)
0.11922846888645108 + 0.14657028641952727im

julia> besselj(2.5, 18.0)
0.11922846888645013

julia> bessely(2.5, 18.0)
0.1465702864195262

Even with the code base as is this is pretty fast.

julia> @benchmark test(2.5, x) setup=(x=6.0 + 10.0*rand())
BenchmarkTools.Trial: 10000 samples with 690 evaluations.
 Range (min … max):  181.341 ns … 285.023 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     183.636 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   189.049 ns ±  10.150 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄█▄█  ▃▁▃              ▄▆▁▆  ▂ ▁   ▃▁▃                        ▂
  ████▆▆█████▇▇▄▁▄▄▃▄▁▅▁▄████▇▆█▇█▇▆▆███▇▄▅▆▆▅▄▅▄▃▄▁▄▁▃▁▁█▇▄███ █
  181 ns        Histogram: log(frequency) by time        225 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

So that's much faster than the current approach. I think some more care will help the convergence by only doing continued fraction for small orders and scaling from there.

The question to me is the stability of the recurrence. For bessely i'm not as worried because the forward recurrence is stable but I think this just needs some more testing but this is pretty nice because it again depends on the accuracy and speed of the same expansion that a lot of the routine does. So improving that will improve the whole routine....

[heltonmc]

@cgeoga To continue this discussion a little bit because I think it would be really nice to round out the modified bessel functions completely. I've also been thinking a lot about how to ensure the implementations are differentiable. Enzyme has come a long way on the main branch and can give pretty fast derivatives of the code in Bessels.jl already.

Anyway for computation of modified bessel functions (real arguments for now) there are really three main approaches in order of precedence.

1. large argument expansion
2. power series
3. uniform large argument expansion for large orders

Now the power series for besseli is always numerically stable. And with some of the recent work with SIMD we should probably unroll this loop by a factor of four and use SIMDMath.jl to make sure that vectorizes. That will make computation of the power series much more efficient. But where all of those methods fail (let's call it the intermediate zone) we need a good fourth algorithm. At first it doesn't seem like there is a great method but this thread mentions two:

1. besseli power series coupled with continued fraction for besselk(v+1)/besselk(v) and wronskian.
2. The slightly modified Miller's scheme to use with forward recurrence and continued fraction with normalizing condition with large argument expansion

Those methods advantages and disadvantages were extensively discussed above but essentially method 1 converges faster for smaller arguments while method two converges faster for larger arguments. The slight disadvantage of method 2 is that the accuracy is tied to the uniform expansions which is slightly less accurate (due to roundoff) than the power series. Though approach 1 is not feasible for complex numbers.

I'm going to introduce a different approach that could be a more robust approach meaning that it will work for complex numbers and be accurate. This was inspired by https://arxiv.org/pdf/1112.0072.pdf when I was working on the airy functions but approach three is:

3. Sequence transformation (Levin's approach) of the asymptotic expansions

At first I ignored these approaches because they honestly just didn't seem worth it and not efficient at all. Essentially they are able to sum a diverging asymptotic series and increase the validity of the asymptotic expansions. Computationally they can be computed using recursive schemes but they require significant more computation. Essentially you must compute a Levin table and compute down the counter diagonals of the Levin table. Now just thinking theoretically about this problem this requires at a minimum N*(N+1)/2 muladd operations. I ran some tests and typically you need about 16-24 terms to get this to work. So for 16 terms you need 136 muladd operations and for 24 terms you need 300 muladds. And this is at a minimum. Not to mention the sequence transformation is essentially a more advanced Pade approximant so you must also compute a numerator and denominator so all those operations need to be doubled.

This is probably why I initially avoided this approach was that most of routines are under 100 ns and not only do you need all these arithmetics you also need to store tables and so your cache lines and memory accesses are honestly rubbish my first attempt at this using arrays took like 10 microseconds which is awful. But I want to come back to the performance after discussing where this is accurate.

So here are a few plots of the relative errors as a function of x at a specific nu value. I'll then plot the number N of terms used in each sequence transformation.

So the really cool thing here is we are summing the asymptotic expansion for large x but this sequence transformation allows us to sum a rapidly diverging series even when x << nu.

Ok for the benchmarks it's hard to say as the general form was so slow but I completely unrolled a version (which is probably not what we want) for N=16 but to really see the limits of the performance. Here is the code it is long (@oscardssmith maybe if you have an idea on how to organize the memory most effectively in a loop here 😄 ) ill work on optimizing the general case in a bit. But the benchmarks are pretty good.

julia> p = ntuple(i -> i*2.0, 18)
(2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0)

julia> @benchmark g1($p)
BenchmarkTools.Trial: 10000 samples with 990 evaluations.
 Range (min … max):  45.328 ns … 84.385 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     45.455 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   45.635 ns ±  1.419 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃█▆▄▂                                   ▁                   ▁
  █████▄▁▁▁▁▁▁▃▁▁▁▃▁▄▁▃▃▁▃▅▃▁▁▁▁▃▃▁▄▁▁▁▁▃▇██▇▄▄▁▃▃▁▃▃▄▁▁▄▃▃▅▅ █
  45.3 ns      Histogram: log(frequency) by time      49.5 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

And the advantage is this isn't awful for complex numbers because we have only real/complex multiplies and complex/complex additions.

julia> p = ntuple(i -> i*2.0 + rand()*im, 18)
(2.0 + 0.20671806918378222im, 4.0 + 0.5544791253581327im, 6.0 + 0.8743485225845774im, 8.0 + 0.5842694164427652im, 10.0 + 0.5376242211139655im, 12.0 + 0.24931450428480173im, 14.0 + 0.7889137184169184im, 16.0 + 0.18129161121223591im, 18.0 + 0.9320649223876125im, 20.0 + 0.4467818886051006im, 22.0 + 0.44671130589145225im, 24.0 + 0.09229756644375375im, 26.0 + 0.09649719037936122im, 28.0 + 0.7052070903007531im, 30.0 + 0.8544141179769782im, 32.0 + 0.35520748181052775im, 34.0 + 0.22437122283477695im, 36.0 + 0.9058299917597746im)

julia> @benchmark g1($p)
BenchmarkTools.Trial: 10000 samples with 832 evaluations.
 Range (min … max):  147.886 ns … 208.534 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     148.689 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   149.701 ns ±   3.756 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃▇█▆▂     ▁▃▄▂▁▁ ▁▁                                           ▂
  █████▇▇█▇▇██████████▇▆▅▆▆▅▆▆▇▆▅▅▅▃▆▆▆▄▅▄▅▅▄▅▁▅▆▅▃▅▆▆▆▄▆▅▅▅▆▄▅ █
  148 ns        Histogram: log(frequency) by time        168 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

const CC =  ((1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0),
(0.5, 0.6, 0.6666666666666666, 0.7142857142857143, 0.75, 0.7777777777777778, 0.8, 0.8181818181818182, 0.8333333333333334, 0.8461538461538461, 0.8571428571428571, 0.8666666666666667, 0.875, 0.8823529411764706, 0.8888888888888888),
(0.32, 0.4166666666666667, 0.4897959183673469, 0.546875, 0.5925925925925926, 0.63, 0.6611570247933884, 0.6875, 0.7100591715976331, 0.7295918367346939, 0.7466666666666667, 0.76171875, 0.7750865051903114, 0.7870370370370371),
(0.23148148148148148, 0.31486880466472306, 0.3828125, 0.438957475994513, 0.486, 0.5259203606311045, 0.5601851851851852, 0.589895311788803, 0.6158892128279884, 0.6388148148148148, 0.6591796875, 0.6773865255444739, 0.693758573388203),
(0.17992503123698458, 0.251220703125, 0.31214753848498705, 0.3645, 0.4098080732190424, 0.44931520061728397, 0.48401666608312033, 0.514707413577676, 0.5420246913580247, 0.5664825439453125, 0.5884987009255157, 0.6084152568206066),
(0.14654541015625, 0.20809835898999135, 0.26244, 0.3104606615295776, 0.35303337191358025, 0.3909365379902126, 0.4248378651752246, 0.45530074074074073, 0.48279762268066406, 0.5077243694259351, 0.5304133008179648),
(0.12331754606814303, 0.177147, 0.22578957202151098, 0.2696782702117627, 0.30931242566258577, 0.34518076545487, 0.3777309849108368, 0.4073604941368103, 0.4344165727708536, 0.4592004039488861),
(0.1062882, 0.15394743546921202, 0.1977640648219593, 0.237932635125066, 0.27473571128040675, 0.30848030434385004, 0.3394670784473419, 0.3679763910529583, 0.39426297308742747),
(0.09330147604194668, 0.13596279456509702, 0.17570409978466411, 0.21259310991936234, 0.24678424347508002, 0.27846908778883517, 0.30784952976979524, 0.33512352712431337),
(0.08308837445644818, 0.12164129985092131, 0.15792631022581202, 0.19194330048061778, 0.22376980268745683, 0.25352314216336075, 0.28133826968460873),
(0.07485618452364388, 0.10998439462154766, 0.14331766435886129, 0.17482015834957565, 0.20452287098892968, 0.2324948200865864),
(0.06808557762286284, 0.1003223650512029, 0.13111511876218174, 0.1604100948932782, 0.18821009245104614),
(0.06242280492074847, 0.09219031787965903, 0.12077936556670357, 0.14812831350313815),
(0.057618948674786896, 0.08525602275296724, 0.11191917020237105),
(0.053493975060685324, 0.07927607889334616),
(0.049914568192106844,))


function g1(s)
    b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, b11, b12, b13, b14, b15, b16, b17 = 1 / (s[2] - s[1]), 1 / (s[3] - s[2]), 1 / (s[4] - s[3]), 1 / (s[5] - s[4]), 1 / (s[6] - s[5]), 1 / (s[7] - s[6]), 1 / (s[8] - s[7]), 1 / (s[9] - s[8]), 1 / (s[10] - s[9]), 1 / (s[11] - s[10]), 1 / (s[12] - s[11]), 1 / (s[13] - s[12]), 1 / (s[14] - s[13]), 1 / (s[15] - s[14]), 1 / (s[16] - s[15]), 1 / (s[17] - s[16]), 1 / (s[18] - s[17])
    a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15, a16, a17 = s[1] * b1, s[2] * b2, s[3] * b3, s[4] * b4, s[5] * b5, s[6] * b6, s[7] * b7, s[8] * b8, s[9] * b9, s[10] * b10, s[11] * b11, s[12] * b12, s[13] * b13, s[14] * b14, s[15] * b15, s[16] * b16, s[17] * b17

    index = 1

    a1 = muladd(a1, -CC[index][1], a2)
    b1 = muladd(b1, -CC[index][1], b2)

    a2 = muladd(a2, -CC[index][2], a3)
    b2 = muladd(b2, -CC[index][2], b3)

    a3 = muladd(a3, -CC[index][3], a4)
    b3 = muladd(b3, -CC[index][3], b4)

    a4 = muladd(a4, -CC[index][4], a5)
    b4 = muladd(b4, -CC[index][4], b5)

    a5 = muladd(a5, -CC[index][5], a6)
    b5 = muladd(b5, -CC[index][5], b6)

    a6 = muladd(a6, -CC[index][6], a7)
    b6 = muladd(b6, -CC[index][6], b7)

    a7 = muladd(a7, -CC[index][7], a8)
    b7 = muladd(b7, -CC[index][7], b8)

    a8 = muladd(a8, -CC[index][8], a9)
    b8 = muladd(b8, -CC[index][8], b9)

    a9 = muladd(a9, -CC[index][9], a10)
    b9 = muladd(b9, -CC[index][9], b10)

    a10 = muladd(a10, -CC[index][10], a11)
    b10 = muladd(b10, -CC[index][10], b11)

    a11 = muladd(a11, -CC[index][11], a12)
    b11 = muladd(b11, -CC[index][11], b12)

    a12 = muladd(a12, -CC[index][12], a13)
    b12 = muladd(b12, -CC[index][12], b13)

    a13 = muladd(a13, -CC[index][13], a14)
    b13 = muladd(b13, -CC[index][13], b14)

    a14 = muladd(a14, -CC[index][14], a15)
    b14 = muladd(b14, -CC[index][14], b15)

    a15 = muladd(a15, -CC[index][15], a16)
    b15 = muladd(b15, -CC[index][15], b16)

    a16 = muladd(a16, -CC[index][16], a17)
    b16 = muladd(b16, -CC[index][16], b17)

    index = 2

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]

    a13 = a14 - a13*CC[index][13]
    b13 = b14 - b13*CC[index][13]

    a14 = a15 - a14*CC[index][14]
    b14 = b15 - b14*CC[index][14]

    a15 = a16 - a15*CC[index][15]
    b15 = b16 - b15*CC[index][15]


    index = 3

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]
    
    a13 = a14 - a13*CC[index][13]
    b13 = b14 - b13*CC[index][13]

    a14 = a15 - a14*CC[index][14]
    b14 = b15 - b14*CC[index][14]

    index = 4

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]
    
    a13 = a14 - a13*CC[index][13]
    b13 = b14 - b13*CC[index][13]
    
    index = 5

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]

    a12 = a13 - a12*CC[index][12]
    b12 = b13 - b12*CC[index][12]

    index = 6

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]

    a11 = a12 - a11*CC[index][11]
    b11 = b12 - b11*CC[index][11]


    index = 7

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]

    a10 = a11 - a10*CC[index][10]
    b10 = b11 - b10*CC[index][10]


    index = 8

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    a9 = a10 - a9*CC[index][9]
    b9 = b10 - b9*CC[index][9]


    index = 9

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    a8 = a9 - a8*CC[index][8]
    b8 = b9 - b8*CC[index][8]

    index = 10

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    a7 = a8 - a7*CC[index][7]
    b7 = b8 - b7*CC[index][7]

    index = 11

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    a6 = a7 - a6*CC[index][6]
    b6 = b7 - b6*CC[index][6]

    index = 12

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    a5 = a6 - a5*CC[index][5]
    b5 = b6 - b5*CC[index][5]

    index = 13

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    a4 = a5 - a4*CC[index][4]
    b4 = b5 - b4*CC[index][4]

    index = 14

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    a3 = a4 - a3*CC[index][3]
    b3 = b4 - b3*CC[index][3]

    index = 15

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    a2 = a3 - a2*CC[index][2]
    b2 = b3 - b2*CC[index][2]

    index = 16

    a1 = a2 - a1*CC[index][1]
    b1 = b2 - b1*CC[index][1]

    return a1 / b1

end

[oscardssmith]

Isn't this just 16 pairs of evalpolys of degrees 1 to 16?

[cgeoga]

This is insanity. But I'm into it. Would you mind sharing your un-optimized g1? I would appreciate taking a look at that to get a better idea of what's going on. Will also take a look at this paper---this is super interesting. I haven't actually seen this Levin transformation before at all.

Also, you mention differentiability with Enzyme---that's pretty cool! Would you also mind sharing those snippets/tests? I just tried in a fresh environment pulling the most recent versions and didn't have much luck for any routines that weren't the fixed order variety (besselk0, besselk1, etc).

[heltonmc]

Isn't this just 16 pairs of evalpolys of degrees 1 to 16?

Yes - but I think the difference here is that you have to store the n-1 values after each iteration so it is more of a cumsum in this regard but I'll look into this to see more. The paper that describes this approach is https://arxiv.org/pdf/math/0306302.pdf (section 7) they actually gave a fortran snippet in 7.5 but I didn't look at it much. @cgeoga I'd highly recommend giving that paper a read it has a lot of very good stuff. @oscardssmith the Levin table is described in 7.5-1 computed with a 3 term recurrence. The code I listed above was the fastest way I could come up with doing that recursive tree type problem (though much of my time was spent just getting this to work)

This is insanity. But I'm into it.

Haha ya me too 😄 The code was originally based on https://github.com/MikaelSlevinsky/SequenceTransformations.jl and once you peel away all the generated functions and what have you you get a code that looks like...

function g(s)
    n = length(s) - 1
    l = similar(s)
    N = similar(l)
    D = similar(N)
    
    for k in 1:n
        w = s[k+1] - s[k]
        N[k], D[k] = s[k] / w, 1.0 / w
        kp1 = k + 1.0
        kdkp1 = k/kp1
        kdkp1jm2 = 1/kdkp1
        for j in 1:k-1
            cst = -(kp1-j)*kdkp1jm2/kp1
            N[k-j] = muladd(cst,N[k-j],N[k-j+1])
            D[k-j] = muladd(cst,D[k-j],D[k-j+1])
            kdkp1jm2 *= kdkp1
        end
        l[k] = N[1]/D[1]
    end
    return l
end

Of course I decided to just completely unroll everything and precompute the cst variable.

I just tried in a fresh environment pulling the most recent versions and didn't have much luck for any routines that weren't the fixed order variety (besselk0, besselk1, etc).

Ya I filed a few issues a couple of weeks ago and Billy fixed them on main. They are working on a new rule system I believe so it's been a while since they have released a tagged version so you'll have to try on Enzyme#main.

julia> @benchmark besselj(v, 1.2) setup=(v=rand())
BenchmarkTools.Trial: 10000 samples with 981 evaluations.
 Range (min … max):  62.096 ns … 102.404 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     67.406 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   65.869 ns ±   2.674 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▁▄▆▇▅                             ▅█▆▂                ▁▂  ▂
  ▅████████▄▄▃▁▁▁▄▄▄▃▃▁▄▁▅███▆▆▃▁▁▁▃▃▅█████▅▄▃▁▄▁▁▁▄▄▃▄▁▃▁▇███ █
  62.1 ns       Histogram: log(frequency) by time      70.5 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark autodiff(Forward, v -> besselj(v, 1.2), Duplicated, Duplicated(2.4, 1.0))
BenchmarkTools.Trial: 10000 samples with 952 evaluations.
 Range (min … max):  94.625 ns … 137.693 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     94.889 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   95.189 ns ±   1.758 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▄▇█▇▄▁                             ▂▃▂                      ▂
  ▇██████▅▁▁▃▄▁▄▁▅▁▁▁▁▄▄▁▁▁▃▃▄▁▄▃▃▁▁▅▇████▄▄▁▃▃▃▃▁▁▄▃▁▄▄▅▆▆▆▆▆ █
  94.6 ns       Histogram: log(frequency) by time      99.6 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> 

Here is the code working on besselj with respect to order on the main branch. Honestly that is pretty fast all things considered.

[oscardssmith]

Yeah, my suggestion was the hybrid of this looped verison with your unrolled version (where you precompute cst but otherwise leave it looped). I think that should vectorize well and be fast.

[heltonmc]

Ya I think my initial attempt at this removing that cst steered me away from this approach considering the n(n+1)/2 growth in operations for the sum. Just an example code..

function g(s)
    n = length(s) - 1
    l = similar(s)
    N = similar(l)
    D = similar(N)
    
    for k in 1:n
        w = s[k+1] - s[k]
        N[k], D[k] = s[k] / w, 1.0 / w
        for j in 1:k-1
            N[k-j] = muladd(-2.1,N[k-j],N[k-j+1])
            D[k-j] = muladd(-2.1,D[k-j],D[k-j+1])
        end
        l[k] = N[1]/D[1]
    end
    return l
end

julia> p = rand(16);

julia> @benchmark g($p)
BenchmarkTools.Trial: 10000 samples with 200 evaluations.
 Range (min … max):  406.665 ns …  2.451 μs  ┊ GC (min … max): 0.00% … 82.54%
 Time  (median):     416.250 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   423.480 ns ± 92.645 ns  ┊ GC (mean ± σ):  1.09% ±  4.06%

     ▁▃▅▆████▇▆▅▄▃▂▁▁▁▂▂▂▂▂       ▁▁▁▁                         ▂
  ▅▆█████████████████████████▇█▇▇███████▇▇▇▆▅▆▆▆▆▅▆▆▅▂▄▄▅▃▄▄▆▅ █
  407 ns        Histogram: log(frequency) by time       466 ns <

 Memory estimate: 576 bytes, allocs estimate: 3.

julia> p = rand(20);

julia> @benchmark g($p)
BenchmarkTools.Trial: 10000 samples with 171 evaluations.
 Range (min … max):  632.795 ns …   2.940 μs  ┊ GC (min … max): 0.00% … 77.83%
 Time  (median):     647.421 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   656.425 ns ± 106.955 ns  ┊ GC (mean ± σ):  0.79% ±  3.76%

       ▂▂▄▆▇█▇▇▆▄▃▁▁▁▁▁▃▃▃▂▁▁    ▁▁▁▂                           ▂
  ▂▄▅▆▇███████████████████████████████▇▇▇▇▅▆▆▆▆▆▅▆▅▅▄▅▄▄▂▅▄▃▅▄▄ █
  633 ns        Histogram: log(frequency) by time        713 ns <

 Memory estimate: 672 bytes, allocs estimate: 3.

julia> p = rand(24);

julia> @benchmark g($p)
BenchmarkTools.Trial: 10000 samples with 45 evaluations.
 Range (min … max):  908.333 ns …  10.032 μs  ┊ GC (min … max): 0.00% … 90.52%
 Time  (median):     934.244 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   945.136 ns ± 218.800 ns  ┊ GC (mean ± σ):  0.56% ±  2.21%

        ▁▄▄▅▆▇███▇▆▄▂▁    ▁▃▃▄▃▃▃▁▁                             ▂
  ▄▃▄▄▆▇██████████████▇▅▆▇██████████████▇█▇▇█▇█▇▇█▆▅▅▄▃▄▅▅▅▄▆▅▆ █
  908 ns        Histogram: log(frequency) by time       1.02 μs <

 Memory estimate: 768 bytes, allocs estimate: 3.

So the unrolled tuple version is over 100x faster (of course this is much more memory friendly)

[cgeoga]

Would that vectorize better than a fully unrolled one? Looking at the crazy g1, it looks to me that you could write that function without a ton of repetition using Base.Cartesian.@nexprs for the declaration and the loops and stuff.

[heltonmc]

Both of the versions vectorize in a sense that the numerator and denominator are computed simultaneously. There is an opportunity to further vectorize this if we want to store two different arrays with offset indexes that allow for computation of the whole vector but if we flush the SIMD registers for just a numerator (say we reduce the total computations by 4x) we still have to compute that twice for denominator and numerator since that is not going to be able to vectorize well.

Though I don't think it becomes worth it to go much further than to SIMD the numerator and denominator (which is the natural choice) because no matter what we have to compute 16->15->14 multiplies so if we unroll by a factor of 4 that means we will have 4->4->4->4->3.... multiply structure. As we will have excess multiples in the registers that won't quite fill them all. So if the naive case is (n*(n+1)/2) multiplies (so 136 for n=16) then the perfect SIMD version would have (4 * 4 + 4 * 3 ...) so 40 multiplies. Now the naive case for the 136 non SIMD version will vectorize the numerator and denominotor but if we flush completely the single evaluation then we have 40*2 so the SIMD has a perfect case for 80 multiplies compared to 136.

That's a lot but it starts to become not worth it because we have to store additional vectors and then have them packed efficiently I doubt that it will help much when all said and done. I initially tried it but I then said what happens when we move to the complex domain? It will make more sense to keep the extra SIMD registers to compute those complex additions better.

t looks to me that you could write that function without a ton of repetition using Base.Cartesian.@nexprs

I think you are right here. My initial attempt was not type stable as the tuple length changed in each iteration so need to be aware of that.

[cgeoga]

Wow, well honestly a lot of that is pretty over my head here. I have never really taken SIMD seriously. I'm trying to write the @nexprs equivalent of your g1 and I have one stumper:

function g1_gen(s)
  # declare the a and b scalar variables:
  Base.Cartesian.@nexprs 17 j->(b_j = inv(s[j+1] - s[j]))
  Base.Cartesian.@nexprs 17 j->(a_j = b_j*s[j])
  # first sequence of operations for index=1:
  CC1 = CC[1]
  Base.Cartesian.@nexprs 16 j->begin
    a_j = muladd(a_j, -CC1[j], a_{j+1})
    b_j = muladd(a_j, -CC1[j], b_{j+1})
  end
  # outer loop over "index". 
  # STUMPER: Is there some way to make sure that this gets
  # expanded before the "inner" @nexprs, so that {16-index_m1} expands to the
  # raw int before the inner call to @nexprs?
  Base.Cartesian.@nexprs 15 index_m1->begin
    index = index_m1+1
    # inner loop over the variable index:
    Base.Cartesian.@nexprs {16-index_m1} j->begin
      a_j = a_{j+1}*CC[index][j]
      b_j = b_{j+1}*CC[index][j]
    end
  end
  a_1/b_1
end

But if there is some trick to make that outer @nexprs expand first, then I think this will give what you want. I'm sure it will need to be workshopped plenty, but just playing around for fun. I tried putting a @macroexpand on the outer @nexprs, but that didn't change anything.

[heltonmc]

Thanks for the idea 😄 I think this can go into an example of the what the heck / or amazing examples of julia but

function g1(s)
    @nexprs 17 i -> b_{i} = 1 / (s[i + 1] - s[i])
    @nexprs 17 i -> a_{i} = s[i] * b_{i}
    @nexprs 16 index -> (@nexprs 17-index i -> (a_{i}, b_{i}) = (muladd(a_{i}, -CC[index][i], a_{i+1}), muladd(b_{i}, -CC[index][i], b_{i+1})))
    return a_1 / b_1
end

I believe that is correct (at least it returns correct answer ) and what you were looking for in like three lines of code!

julia> @benchmark g1($(0.11274150564997776, 0.11475646021904119, 0.11370868384312821, 0.11451023277070166, 0.11369906525599734, 0.11472113632452476, 0.11317926911257481, 0.11588910073757686, 0.11045137194340607, 0.12271888810305538, 0.09198318397942494, 0.17666004884002684, -0.07776136204113551, 0.7501986007978472, -2.150973108989948, 8.739300196625988, -34.86022928456346, 150.58310277542904))
BenchmarkTools.Trial: 10000 samples with 991 evaluations.
 Range (min … max):  41.835 ns … 73.957 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     41.961 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   42.132 ns ±  1.291 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▄█▆▃▂                                                       ▁
  █████▅▁▁▁▁▁▁▁▄▁▁▁▃▃▁▃▃▄▄▄▄▃▅▄▁▁▁▄▃▁▁▄▁▅███▅▅▇▄▁▁▁▁▃▄▁▅▃▁▃▃▄ █
  41.8 ns      Histogram: log(frequency) by time        46 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[heltonmc]

I think we could use a generated function as well based on the NTuple length of N.

@oscardssmith I'm curious what you feel like including code like that in the package considering compile times / package load times / cache pressure? The generated llvm code is immense and how does that affect other code who may call besselk amongst other calls within a larger function. It is definitely very very fast in these microbenchmarks...

[oscardssmith]

in terms of speed, the generated function + nexprs approach will be really good. and I don't think the cache pressure/load times will be much worse than the vectorized evalpolys (which also unroll). The one thing you could do that would help is generating the n^2 nexprs as vectorized muladds in the first place (since doing so will reduce the number of instructions by 4x/8x).