Hacker News Reader

Comment by plg | Hacker News Reader

That (last method) is a lot of sines and cosines, cube root and an acos... I bet for 3D the rejection method is not that much slower, in wall clock time.

OK so I coded up the first one and the last one in C and compiled on my mac using whatever gcc links to. I found 1e7 points in the sphere:

first routine: 0.521s

last routine: 0.850s

[edit]: C code here: https://pastebin.com/zH0BehSM

View on HN

Note that his last way does cos(acos(cosϕ))) after generating cosϕ uniformly. Which is a waste, you can just set sinϕ = sqrt(1-cosϕ^2).

After fixing that, this way is faster for me. Some trig, but no branches.

Edit: https://pastebin.com/UExZ5Siw

But you used Julia. I'm almost certain that if you applied your optimization to the C code you'd find that it makes no difference because `gcc -O2` already made it.

Gcc collapses acos(cos(x)) to x on -O2? Surely not since this isn't even accurate for x outside of -pi to pi. That's ignoring floating point errors that make those two different which I thought was a different optimization flag, though I'm not sure about that.

Well, here it's cos(acos(x)) and sin(acos(x)), which are easier. But still that's pretty smart, I did not know this.

Julia is usually ran with LLVM -O3. Is LLVM just missing that optimization then? This sounds like something the optimizer wouldn't due without fastmath flags enabled.

The first routine is about 3 rands per point minus rate of discarded, about half? So 6 rands per point. Now, since the second one has three rands in it, and plenty of trig, it really becomes a question if you can make the trigonometry faster than three rands. But as you saw, straightforward attempt will pretty much certainly end up having worse performance.

But what if one used optimizations like lookup tables and trigonometric approximations instead? Could they end up costing less cycles than three rands? Maybe not LUTs, maybe other trickery someone knows of?

And what about the cost of those three rands? Could one use a fast xorshift with good enough results here? Does one even need thee random vectors? Maybe two + modulo trickery will do?

Typically the generation of the uniform random numbers can be made quite cheaply, especially if the resolution is not critical (e.g. if uniform floats are enough and not uniform doubles etc) or if the quality of the randomness is not critical, and in these cases discarding half of these will win versus other costly operations.

If we are trying to optimize anything, you should probably drop all the unnecessary floating point divisions. I would expect those to be the bottleneck in your code, though you’d probably want to time it to be sure.

(rand()/(double)RAND_MAX) etc.

Instead assign some name to 1.0/(double)RAND_MAX, and multiply.

Reply on news.ycombinator.com