Cauchy distribution: optimise sampling

[dhardy]

@MaximoB this presumably uses 1 bit less precision but is quite a bit faster

# previously
test distr_cauchy                ... bench:      48,639 ns/iter (+/- 1,411) = 164 MB/s                                                 
# now
test distr_cauchy                ... bench:      36,224 ns/iter (+/- 654) = 220 MB/s

Optimisers are weird, but as I guessed, it seems we get a free subtraction operation when using Open01.

I considered checking that the result is finite with a loop, but I'm fairly sure it must be anyway since π is irrational. Do you agree?

[vks]

I considered checking that the result is finite with a loop, but I'm fairly sure it must be anyway since π is irrational.

I think you are right that tan produces finite results for finite input. The only thing I can think of that could go wrong is that r * pi / 2. rounds such that you end up on the other sight of the singularity, changing the sign.

[dhardy]

I realised that might be possible, but frankly I don't think it would have much effect on the distribution (it's valid to sample from any π-length interval excluding the singularities). Probably though it would round towards zero anyway.

[pitdicker]

Optimisers are weird, but as I guessed, it seems we get a free subtraction operation when using Open01.

75% plus of the time is spend in tan. LLVM can't combine both subtractions, because it doesn't know there will be no rounding error. If you use the range code, or use IntoFloat directly, this one instruction can be removed, but it doesn't matter in practice. The performance difference is because of using Open01, and removing the loop.

[vks]

I don't think it would have much effect on the distribution

I agree, it would just wrap around, which is not a problem because it is a random number anyway.

LLVM can't combine both subtractions, because it doesn't know there will be no rounding error.

More concretely, it cannot optimize x - (1.0 - EPSILON / 2.0) - 0.5 to x - ((1.0 - EPSILON / 2.0) - 0.5).

[MaximoB]

@dhardy
Ok, I am not sure how I got such different results in my benchmarking, but the two methods were very close on my computer so if this way is significantly faster for you then it is probably the better option.

Since you are using Open01 the only way you could hit one of the singularities in tan is if Open01 malfunctions or there is a floating point error subtracting 0.5 so I don't really think getting infinity is something to be concerned about.

Over the weekend I read some papers and I have implemented a stunningly fast numerical approximation method for generating Cauchy numbers.
test distr_cauchy ... bench: 9,945 ns/iter (+/- 2,510) = 804 MB/s
I don't really know what the error bounds on this is though so if you are more comfortable going with a more accurate but slower generation I understand (although the paper claims it is very low error).

[dhardy]

The performance difference is because of using Open01, and removing the loop.

You're right; turns out removing the loop makes most of the difference. Moving everything inside the loop also performs well on my CPU, but I think we can just drop the loop since tan(π/2) rounds to a finite number anyway.

[dhardy]

Amended.

@MaximoB that sounds worth looking into. I don't think this is the only distribution which could be sped up significantly using a different approach; see #257. If you're interested then please open a PR.