Better precision for cispi(::Complex)

[antoine-levitt]

See https://discourse.julialang.org/t/cispi-is-not-accurate/83811/.

[giordano]

Can add any tests?

[antoine-levitt]

Done, thanks for nudging me. I also fixed thanks to the existing tests a very fun bug: exp(-pi*v) is wrong because it's really exp((-pi)*v), which converts pi to Float64 (which is bad for higher precisions)

[giordano]

Sounds like you have to update a doctest as well: https://buildkite.com/julialang/julia-master/builds/13663#0181e2d3-4b6d-4359-8341-a387a685cc86/1260-1268.

Just to confirm the goodness of the change, the new result seems to be closer to what multiple precision gives:

julia> ComplexF64(cispi(big(0.25 + 1im)))
0.030556854645954555 + 0.030556854645954555im

[oscardssmith]

technically you could do this faster by combing the reduction steps, but that is probably overkill here.

[giordano]

The new implementation is already faster than the current one, in addition to being more accurate:

julia> using BenchmarkTools

julia> function cispi_new(z::Complex)
           v = exp(-(pi*imag(z)))
           s, c = sincospi(real(z))
           Complex(v * c, v * s)
       end
cispi_new (generic function with 1 method)

julia> @benchmark cispi(z) setup=(z=rand(ComplexF64))
BenchmarkTools.Trial: 10000 samples with 992 evaluations.
 Range (min … max):  35.504 ns … 112.049 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     42.331 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   42.003 ns ±   4.903 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

    ▇              █   ▄▂                                       
  █▁█▂▃█▂▂▂▂▂▂▁▁▃▇▁█▂▂▇██▇▃▂▃█▂▂▃▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ ▃
  35.5 ns         Histogram: frequency by time         58.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark cispi_new(z) setup=(z=rand(ComplexF64))
BenchmarkTools.Trial: 10000 samples with 996 evaluations.
 Range (min … max):  22.240 ns … 72.815 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     23.563 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   24.159 ns ±  3.092 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark w .= cispi.(z) setup=(z=randn(ComplexF64, 1_000_000); w=similar(z))
BenchmarkTools.Trial: 75 samples with 1 evaluation.
 Range (min … max):  49.548 ms … 71.129 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     50.281 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   52.622 ms ±  3.839 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                        ▂                                   
  █▇▄▄▃▁▃▃▁▁▁▁▃▁▁▁▁▁▁▃▁▁▁▁▄█▄▃▁▁▁▃▁▁▁▁▃▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
  49.5 ms         Histogram: frequency by time        63.4 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark w .= cispi_new.(z) setup=(z=randn(ComplexF64, 1_000_000); w=similar(z))
BenchmarkTools.Trial: 120 samples with 1 evaluation.
 Range (min … max):  27.201 ms … 34.093 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     27.489 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   29.053 ms ±  2.356 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

   █▂                                                          
  ▅██▇▃▂▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆█▆▁▂▁▁▂ ▂
  27.2 ms         Histogram: frequency by time          33 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

[antoine-levitt]

Good catch, thanks! (I did look at CI, but it appeared to be broken for some reason so I didn't get very far).