support of cospi and sinpi

[jaakkor2]

Could cospi and sinpi function be supported within @avx?

function f1()
	z = cumsum(ones(4))
	A = zeros(4)
	for i in 1:size(A)[1]
		A[i] = cospi(z[i])
	end
	return A
end

function f1a()
	z = cumsum(ones(4))
	A = zeros(4)
	@avx for i in 1:size(A)[1]
		A[i] = cospi(z[i])
	end
	return A
end

gives

julia> f1()
4-element Array{Float64,1}:
 -1.0
  1.0
 -1.0
  1.0
julia> f1a()
ERROR: MethodError: no method matching cospi(::VectorizationBase.SVec{4,Float64})
Closest candidates are:
  cospi(::T) where T<:AbstractFloat at special/trig.jl:815
  cospi(::Integer) at special/trig.jl:864
  cospi(::T) where T<:Union{Integer, Rational} at special/trig.jl:841
  ...
Stacktrace:
 [1] macro expansion at C:\Users\jr\.julia\packages\LoopVectorization\9Qgzg\src\reconstruct_loopset.jl:232 [inlined]
 [2] _avx_! at C:\Users\jr\.julia\packages\LoopVectorization\9Qgzg\src\reconstruct_loopset.jl:232 [inlined]
 [3] macro expansion at .\gcutils.jl:91 [inlined]
 [4] f1a() at .\REPL[35]:4
 [5] top-level scope at REPL[37]:1

Sometimes sinpi can be faster, I wonder if this applies to @avx

using BenchmarkTools
a=10e6*rand(10_000_000);
b=π*a;
@btime sinpi.($a)
@btime sin.($b);

gives

  167.324 ms (2 allocations: 76.29 MiB)
  322.987 ms (2 allocations: 76.29 MiB)

[chriselrod]

This is ultimately a problem with SLEEFPirates.jl, which forked and modified from SLEEF.jl, which is itself a port of an old version of SLEEF, which is a library written in C.

This old port does not have a sinpi function, so there currently is no method for sinpi(::SVec).
It would be trivial to add sinpi(v::SVec) = sin(π * v). Of course, this won't have any advantages (speed of accuracy) over calculating sin(π*a) yourself aside from convenience.
It is better than throwing an error, so I'll make that addition.

SLEEF (the C library) does now have a sinpi function (as well as cospi and sincospi). Porting that would be a much more difficult, but better, solution.
Longer term, SLEEFPirates.jl isn't great as far as vectorized math libraries go.

SLEEF is now faster and much more complete than it was then.
Recent versions of the gcc compiler also use a far better library than this old version. xsimd looks good too.
(Re-)porting any of these to Julia make be good.
But those would be major undertakings. I hope someone else will do that someday. For now, it isn't high on my own priority list.

I added the simple sinpi definitions. I'll close this issue after I committed them.

[zsoerenm]

It might not be something you asked for, but I have discovered the following:
If you are really after performance and do not care about accuracy too much you can use a fixed point solution. I found one here: https://www.coranac.com/2009/07/sines/
I used the 4th order implementation:

@inline function sin16384pi(x::Int32)
    qN = 13
    qA = 12
    C = Int32(3692) # floor(Int32, 5 * (1 - 3 / pi) * 1 << 14)
    B = Int32(20076) # C + 1 << 14

    c = x << (30 - qN)
    x -= Int32(1) << qN

    x = x << (31 - qN)
    x = x >> (31 - qN)
    x = (x * x) >> (qN << 1 - 14)

    y = B - (x * C) >> 14
    y = (Int32(1) << qA) - ((x * y) >> 16)
    sign(c) * y
end

@inline function cos16384pi(x::Int32)
    qN = 13
    x += Int32(1) << qN
    sin16384pi(x)
end

function calc_sins2!(sins, delta)
    delta_fp = floor(Int32, delta * 16384 * 1 << 17)
    phase_fp = Int32(0)
    @inbounds @simd for i = 1:length(sins) # Unfortunately I haven't been able to use it with @avx (it throws ERROR: InexactError: trunc(Int64, Inf))
        phase_fp = phase_fp + delta_fp
        phase_int = phase_fp >> 17
        sins[i] = sin16384pi(phase_int)
    end
end

# For comparison
function calc_sins1!(sins, delta)
    @avx for i = 1:length(sins)
        sins[i] = sin(i * delta)
    end
end
sins1 = Vector{Float64}(undef, 1000)
sins2 = Vector{Int32}(undef, 1000)

julia> @btime calc_sins1!($sins1, $(0.01π))
  8.165 μs (0 allocations: 0 bytes)
julia> @btime calc_sins2!($sins2, $0.01)
  408.615 ns (0 allocations: 0 bytes)

If you plot both, there is hardly a difference to be seen:

using PyPlot
plot(sins1)
plot(sins2 ./ 4092)

plot(sins2 ./ 4092 .- sins1)

If you'd like to I can share the Int16 version, which is even faster (but also more inaccurate)

[chriselrod]

I am interested.
It could be added as an option for those to whom accuracy isn't important, eg sin_extremely_fast (vs sin and sin_fast.

But mostly, I am interested in using this with the Box-Muller transform for quickly generating random-normal random variables.

Note that VectorizedRNG produces random unsigned integers.
Currently, I'm transforming these UInts into random floats in [1,2), subtracting from 2 to get (0,1], and then multiplying by 2pi before taking sincos.

Instead, I should be using a fixed point approximation like that. Instead of converting the UInts to floating point values, I should just take some random subset of bits. I think I'd want to use a lot more than just 14 bits -- ideally, something closer to 52.

Unfortunately I haven't been able to use it with @avx (it throws ERROR: InexactError: trunc(Int64, Inf))

Please file issues as they come up! But, I'll look into it.

[zsoerenm]

Here is the code for 64-bit Integer:

@inline function sin1073741824pi(x::Int64)
    qN = 29
    qA = 28
    C = Int64(241969553) # floor(Int64, 5 * (1 - 3 / pi) * 1 << 30)
    B = Int64(1315711377) # C + 1 << 30

    c = x << (62 - qN)
    x -= Int64(1) << qN

    x = x << (63 - qN)
    x = x >> (63 - qN)
    x = (x * x) >> (qN << 1 - 30)

    y = B - (x * C) >> 30
    y = (Int64(1) << qA) - ((x * y) >> 32)
    sign(c) * y
end

@inline function cos1073741824pi(x::Int64)
    qN = 29
    x += Int64(1) << qN
    sin1073741824pi(x)
end

function calc_sins3!(sins, delta)
    delta_fp = floor(Int, delta * 1 << 30 * 1 << 33)
    phase_fp = Int(0)
    @inbounds @simd for i = 1:length(sins)
        phase_fp = phase_fp + delta_fp
        phase_int = phase_fp >> 33
        sins[i] = sin1073741824pi(phase_int)
    end
end
sins3 = Vector{Int64}(undef, 1000)

julia> @btime calc_sins3!($sins3, $0.01)
  1.791 μs (0 allocations: 0 bytes)
julia> plot(sins3 ./ 1 << 28 .- sins1)

[chriselrod]

For now I did the lazy thing, but at least they work.

@zsoerenm Thanks, I'll try updating VectorizedRNG some time. It isn't really any faster than the default mersenne twister we get with the Julia standard library (dSFMT), but I think with this normal random sampling could generally be faster than base (which uses the zigguarat method).

julia> x = collect(0.01 .* (1:1000)); y = similar(x);

julia> sins3 = Vector{Int64}(undef, 1000);

julia> @benchmark @. $y = sin($x)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     6.373 μs (0.00% GC)
  median time:      6.413 μs (0.00% GC)
  mean time:        6.526 μs (0.00% GC)
  maximum time:     12.784 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     5

julia> @benchmark @avx @. $y = sin($x)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.239 μs (0.00% GC)
  median time:      2.251 μs (0.00% GC)
  mean time:        2.255 μs (0.00% GC)
  maximum time:     5.747 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     9

julia> @benchmark @avx @. $y = SLEEFPirates.sin_fast($x)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     952.826 ns (0.00% GC)
  median time:      955.609 ns (0.00% GC)
  mean time:        957.051 ns (0.00% GC)
  maximum time:     1.358 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     23

julia> @benchmark calc_sins3!($sins3, 0.01)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     479.518 ns (0.00% GC)
  median time:      481.594 ns (0.00% GC)
  mean time:        491.488 ns (0.00% GC)
  maximum time:     929.401 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     197

What would be the most efficient sincos implementation here?
Interestingly, these sin and cos implementations are each faster than the vsqrt assembly instruction:

julia> @benchmark @avx @. $y = sqrt($x)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     697.726 ns (0.00% GC)
  median time:      697.753 ns (0.00% GC)
  mean time:        698.740 ns (0.00% GC)
  maximum time:     854.945 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     146

meaning using the identity cos(x) = sqrt(1 - abs2(sin(x))) will not be faster.
Although, 0.001 error still seems fairly high.

But there's got to be a better way to calculate sincos than just doing both individually, no?
I'll definitely read the link before updating VectorizedRNG to see if I can figure anything out.

[zsoerenm]

Although, 0.001 error still seems fairly high.

yes, if you are looking for higher accuracy, you might be interested that Sleef released a fast sin and cos within 350ulp in April last year: shibatch/sleef#229

But there's got to be a better way to calculate sincos than just doing both individually, no?

Yes, you could use the 3rd or 5th order implementation to calculate the sin and the 4th order to calculate the cos. Both implementations need x^2. If both are inlined the compiler should be able to calculate it only once.

Edit:
I found another resource which explains better how to set the bits and minimize the quantization error: https://www.nullhardware.com/blog/fixed-point-sine-and-cosine-for-embedded-systems/
I guess my solution above can be improved. I might come with a solution in the next couple of days.