Optimizing Math Functions

[Jonanin]

Hi all.

A while ago I was browsing the tubes and came across an article on Quake's fast inverse square root function. It was pretty interesting (and very clever...). It calculates 1/sqrt(x) using only bit shift and multiplication. I was thinking about how this could maybe be used in Spring to speed things up. I struck up a conversation in #sy and discussed it a little bit but nothing came out of it other than learning that Spring doesn't use approximations at all.

Today I coded up a fast sqrt functions using the magic number from Quake and newton's method. Note that the sqrt is just an approximation - which is what makes it fast.

Code: Select all

float bsqrt(float number) {
    long i;
    float x, y; 
    const float f = 1.5F; 
    x = number * 0.5F; // half number used for approximation
    y  = number; 
    i  = * ( long * ) &y; 
    i  = 0x5f3759df - ( i >> 1 );    // magic number from Quake 3. this is the initial guess
    y  = * ( float * ) &i; 
    y  = y * ( f - ( x * y * y ) );  // Newton's Method (This is 2 iterations, use more means more accuracy
    y  = y * ( f - ( x * y * y ) );  // but it is slower. Less iterations isn't as accurate but faster.
    return number * y; 
}

As you can see the last two lines are the same. These are the iterations of Newton's method. The more iterations you use, the more accurate the guess. However, more iterations is slower.

So far the tests have been great when comparing with streflop's sqrt. When calculating sqrts of 1-100000, 100 times, the fast sqrt is roughly 25x faster.
See here (Times in milliseconds):

Code: Select all

Streflop Sqrt time: 7563
Fast Sqrt 1 Iteration: 265
Fast Sqrt 2 Iterations: 344
Fast Sqrt 3 Iterations: 453

Now, about how accurate it is. Going from 1 to 10000000:
(Note that 'error' is the difference between the fast sqrt and the real square root. so an error of 1 would mean, say, sqrt(x) was 5 and fastsqrt(x) was 6).

Code: Select all

1 Iteration:
     Average error: 1.9649329910134016863310080225346609950065612792969
     Highest error: 4.830078125 at sqrt(9999917)
2 Iterations:
     Average Error: 0.0036126651906251906927469708108446866390295326709747
     Highest Error: 0.01123046875 at sqrt(9975537)
3 Iterations:
     Average Error: 0.000043552740287780761837405085756813605257775634527206
     Highest Error: 0.000244140625 at sqrt(4196813)

Now, I was wondering if this would be worth implementing, and if so, how I would go about implementing and testing its performance against the old sqrt.

Also thanks to lurker for a lot of help.

Comments, etc, please?

[ILMTitan]

Sounds great. I know we could use it in some of the float3 methods at least.

[Jonanin]

Also, I forgot to add.

Has been tested on both windows and linux, it runs exactly the same on both.

[LordMatt]

So it gives the same error on Linux and windows or runs at the same speed?

[Peet]

So it gives the same error on Linux and windows or runs at the same speed?

Yes.

[Jonanin]

Helpful peet :D

Heh, the code never runs the same speed every time in the SAME platform

I meant it returns the same values on both platforms, so yes, same error.

[KDR_11k]

It might still be useful to keep fastinvsqrt around, things like normalizing need that a lot.

[yuritch]

Absolute error values are not that important, it's relative error that says much about approximate functions usefulness. Can you show relative (i.e. percentage) errors of this compared to 'slow' sqrt? I.e. like that: max relative error - 50%, min relative error - 0.0000001%?
Why is this important: compare an absolute error of 1 when computing sqrt(1) and sqrt(10000). In the first case it's a 100% error (0 or 2 instead of 1), in the second case it's 1% error (99 or 101 instead of 100), but in absolute value it's always a 1.

[koshi]

the maximum relative error over all floating point numbers was 0.00175228

Link: http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf

[imbaczek]

given that sqrt consistently shows up in profiles as one of the most time-consuming functions, I'm all for including and using it.

edit: paper states that 0x5f375a86 is a better initial guess, but it doesn't matter _that_ much.

[Auswaschbar]

I did some investigations in math functions in the past and found out that the system's math functions are about 20 times faster than streflop's, but its more accurate than this approximation, but its unsynced.
How fast is this approximation compared with your systems sqrt()-implementation?

[AF]

We could always put it in a branch and test, that way we wouldn't need to speculate.

[lurker]

The system's? With math.h, it came out at the same speed as with streflop, and I'm reasonable sure there were no auto-includes as not having the "streflop_cond.h" line made it not compile. Do I need to do something else to match the mode Spring's streflop runs in?

Edit: After many errors trying to get streflop_init<streflop::Simple>(); to run, I gave up and put the test code directly into game.cpp.
results in a minute

[Kloot]

@Jonanin

Why long instead of int as in Lomont's paper? Even
approximations should always have equal precision
on all hardware (even though in this case it doesn't
matter).

[lurker]

Okay, finally got it going in spring. A two iteration run is ~8 times faster than streflop.

For ten runs of each number 1 - 10,000,000
Streflop: Sqrt Total Time: 11297ms
Fast Sqrt, 1 iteration: 1156ms
Fast Sqrt, 2 iterations: 1516ms
Fast Sqrt, 3 iterations: 1828ms

[Tobi]

What about __builtin_fsqrt() ?

It's one assembler instruction only, and only one line of code so it's easier to maintain too. (Basically it's what GCC would generate with -ffast-math, I think.)

It may even be (a lot) faster then the Quake bsqrt.

[lurker]

streflop: 11610ms
quake, 2 iterations: 1515ms
__builtin_sqrt: 1344ms
1 / __builtin_sqrt: 2219ms

Why don't we use that again? syncing?

[Tobi]

Yes.

It could be that streflop isn't really needed even for sync since we use only GCC anyway, but it's VERY hard to test.

Basically you need to generate a sync matrix like I once did, ie. test every version of GCC you can get with every optimization flag (0, 1, 2, 3, s) vs each other (well vs one reference version should be enough, especially since the old sync matrix is known: every GCC >= 4.0 with every optimization option syncs against each other), to check whether there are any regressions.

(I think there's a script for that somewhere in tools/, when configured correctly it should be possible to run the script and let it run until it finishes (takes a few days); then you can create a matrix from the log file.)

Of course this isn't needed if you make synced code use streflop and only convert unsynced code, but the problem with this is to figure out the possible context in which each sqrt is executed. (for example the one in float3::Normalize(), which is probably one of the most used sqrt calls, is called from both synced and unsynced contexts, so you can't just change it without extensively testing for sync regressions)

EDIT: Before anyone suggests it, SyncedFloat3 is NOT the solution to the latter issue, since the operations are defined in a way to minimize impact on mixed sync calculations: ie. SyncedFloat3::operator+(SyncedFloat3, SyncedFloat3) returns float3, because it is entirely valid to add two synced floats in unsynced context. Would you then do (syncedFloat1 + syncedFloat2).Normalize() and expect SyncedFloat3::Normalize() to get called, you are wrong

[AF]

we could add a global boolean check specifying wetehr we're in synced or unsynced mode but then I dont know if the overhead of that would negate any benefit in a faster sqrt operation. It would also be useful for other things, perhaps its already done somewhere, bah I dont know enough to comment.

[det]

Today's computers have an SSE function to do an equivalent approximation of inverse square root at greater speeds.