Pythagoras Optimising...

[Xalthorn]

This must have been done countless times, but I'm blowed if I can find out how.

Basically, I'm looking to optimise the following calculation.

A = R / SQR(XX*XX + YY*YY)

If anyone can point me in the direction of a layman's article or simply explain how this is done, that would be perfect

[hellfire]

Can you take any assumptions on XX and/or YY ?
Is the input/output integer or floating-point?
What precision is necessary?

What's expensive here is the division and the square root.
If you're working in floating point, each number is stored in exponential form.
Taking the square root is now very easy, just divide the exponent by 2.
Reciprocal is just as cheap, just negate the exponent.
Notice that you have to renormalize the mantissa when changing the exponent.

edit: wasn't quite right because f= m*2^e and not m^e, so modifying the exponent is a little bit more work.

[Xalthorn]

Can you take any assumptions on XX and/or YY ?
Is the input/output integer or floating-point?
What precision is necessary?

What's expensive here is the division and the suqare root.

If you're working in floating point, each number is stored as:

Code: [Select]

sign * mantissa * 2^exponent
Taking the square root in this exponential form is very easy, just divide the exponent by 2:

Code: [Select]

sign * mantissa * 2^(exponent/2)
Reciprocal is just as cheap:

Code: [Select]

sign * mantissa * 2^(-exponent)
So you just have to extract the exponent (integer) from your floating-point value "x" (as definied in ieee), shift and negate it, put it back into your float and the result is a pretty good approximation of 1/sqrt(x).
Notice that you have to renormalize the mantissa when changing the exponent.

I'm using floating points.

After I posted this, I found the InvSqrt() function that's hugely talked about and that appears to be following the same logic that you're talking about.

However, can I make it work? nope

I'm currently doing this is Freebasic, and I think I'm tying myself in knots.

Code: [Select]

FUNCTION InvSqrt(byval X AS DOUBLE)

' I think X needs to be a positive number
    X=ABS(X)
    DIM AS DOUBLE xhalf
    DIM AS INTEGER i
    xhalf=0.5*X
   
    i=X
   
' assuming the following number is hex, I converted it to decimal in a calculator...
'    i=0x5f3759df-(i SHR 1)
    i=1597463007-(i SHR 1)
   
    x=i
   
    x=x*(1.5-xhalf*x*x)
    return x
   
END FUNCTION


[hellfire]

Code: [Select]

i=X

This does a cast (skip fractional digits of the double).
What you want instead is loading the whole "bit-set" of the double into an integer (but since integers can hold 32bits only you should use "single" as input!), so you can extract the exponent/mantissa parts.
To get it right you have to get a pointer (of type integer) to your single and read (an integer) from the pointer-address.
Do the same (single-pointer to the integer) to store the modified value back into the single.

[Xalthorn]

Code: [Select]

i=X

This does a cast (skip fractional digits of the double).
What you want instead is loading the whole "bit-set" of the double into an integer (but since integers can hold 32bits only you should use "single" as input!), so you can extract the exponent/mantissa parts.
To get it right you have to get a pointer (of type integer) to your single and read (an integer) from the pointer-address.
Do the same (single-pointer to the integer) to store the modified value back into the single.

Hmm,

I'm missing something.  I changed the function to:

Code: [Select]

FUNCTION InvSqrt(byval X AS SINGLE)
    X=ABS(X)
    DIM AS INTEGER PTR P
    DIM AS SINGLE PTR S
    DIM AS SINGLE xhalf
    DIM AS INTEGER i
    xhalf=0.5*X
   
    P=@X
    i=*P
   
'    i=0x5f3759df-(i SHR 1)
    i=1597463007-(i SHR 1)

    S=@i
    X=*S
   
    x=x*(1.5-xhalf*x*x)
    return x
   
END FUNCTION

And calling the function as:

Code: [Select]

            SUM=BR1*InvSqrt(XX*XX+YY*YY)

rather than the working version of

Code: [Select]

            SUM=BR1/SQR(XX*XX+YY*YY)

The variables are defined as:

Code: [Select]

    DIM AS INTEGER XX,YY
    DIM AS SINGLE BR1, SUM


[hellfire]

What's the return-type "InvSqrt" ?

[Xalthorn]

What's the return-type "InvSqrt" ?

I'm asssuming it returns the variable as whatever type the variable is defined as.

[hellfire]

I don't know what your compiler does, but mine (fbc 0.20.0) doesn't even compile:

test.bas(1) error 136: Default types or suffixes are only valid in -lang deprecated

So I specified the return type as single and everything worked fine.

[Xalthorn]

I don't know what your compiler does, but mine (fbc 0.20.0) doesn't even compile:

test.bas(1) error 136: Default types or suffixes are only valid in -lang deprecated

So I specified the return type as single and everything worked fine.

weird, mine compiles and runs but it's just not returning values anywhere near to the sqrt version.

I'll have to dig deeper into the variable storage and start pulling things out of memory more precisely I guess.

[hellfire]

mine compiles and runs but it's just not returning values anywhere near to the sqrt version

Code: [Select]

FUNCTION InvSqrt(byval X AS SINGLE) AS SINGLE
    X=ABS(X)
    DIM AS INTEGER PTR P
    DIM AS SINGLE PTR S
    DIM AS SINGLE xhalf
    DIM AS INTEGER i
    xhalf=0.5*X
   
    P=@X
    i=*P
   
'    i=0x5f3759df-(i SHR 1)
    i=1597463007-(i SHR 1)

    S=@i
    X=*S
   
    x=x*(1.5-xhalf*x*x)
    return x
END FUNCTION
   

Code: [Select]

  x                InvSqrt(x)                 1.0/sqr(x)
 0.0              1.981775e+019              1.#INF       
 0.5              1.41386                    1.41421
 1.0              0.99830                    1.0
 1.5              0.81536                    0.81649
 2.0              0.70693                    0.70710
 2.5              0.63137                    0.63245
 3.0              0.57684                    0.57735
 3.5              0.53442                    0.53452
 4.0              0.49915                    0.5
 4.5              0.47135                    0.47140
 5.0              0.44714                    0.44721
 5.5              0.42605                    0.42640
 6.0              0.40768                    0.40824
 6.5              0.39160                    0.39223
 7.0              0.37744                    0.37796
 7.5              0.36484                    0.36514
 8.0              0.35346                    0.35355
 8.5              0.34276                    0.34299
 9.0              0.33295                    0.33333
 9.5              0.32395                    0.32444
10.0              0.31568                    0.31622

 

[Xalthorn]

Ooh, in that case I'm being dumb

how would I adapt this code...

Code: [Select]

SUM=BR1/SQR(XX*XX+YY*YY)

To work with InvSqrt?  I figured that as 1/SQR is being returned rather than SQR, I would have to flip the calculation.  So I did this (which is clearly not working)

Code: [Select]

SUM=BR1*InvSqrt(XX*XX+YY*YY)


[hellfire]

Your code should work fine. Maybe the result is just not precise enough.
Depending on the range of your input-values you can either tweak the reference value (0x5f3759df) or perform additional approximation steps (see this and this).
Since XX and YY are integers you might want to consider using tables, too.

[stormbringer]

just to satisfy my curiosity... what is the problem you need to solve all this for?

[Xalthorn]

just to satisfy my curiosity... what is the problem you need to solve all this for?

In Shocky's birthday card, I had some meta balls (I think that's what they're called) where they move around and merge into each other.  The calculation is done three times per pixel and as a result it slows the thing down.  I was just looking for a way to optimise things.

I know the optimisation that hellfire is talking about should work, I just can't get it to work.

I think I need to tear that routine out as a stand alone program and post the source, that would make it a lot easier for people to see where I'm missing the obvious

[Voltage]

The code Hellfire posted above works fine for me.  FB version 18.3

[Xalthorn]

Yup, my dumbness knows no bounds.

My *theory* was okay, and my function was *almost* the same as hellfires, with the teeny little difference that the function was declared as:

FUNCTION InvSqrt(byval X AS SINGLE)

When I changed it to

FUNCTION InvSqrt(byval X AS SINGLE) AS SINGLE

I got a Type mismatch compilation error, so I removed it and kept hammering away at the problem.  Major big nooby mistake

It just occurred that I should check my initial function declaration, which oddly enough reads

FUNCTION InvSqrt(byval X AS SINGLE)

I then realised that of course, I had to change them both to show the return type.

Ho hum

It works now, still fairly slow, but better than before.

Thanks for the help and sorry for being so dense

hehe

[Voltage]

I've seen this great code snippet before in C.  It's great to see it working in FB, plus with an explanation I can understand.

Thanks guys.

[Xalthorn]

I figured it would be a cunning plan to put a little test program together that uses this approach.  So here goes.

The zip file contains a compiled executable and the source.

Code: [Select]

' ------------------------------------------------------------------------------
' Blobs - By Xalthorn
'
'  An experiment with meta balls (I think that's what they are called).
'  Basically, we loop through the screen, work out how far each pixel is from
'  each of the balls (well, the centre of them) and colour the pixel according
'  to the sum of these distances.
'
' ------------------------------------------------------------------------------

' ------------------------------------------------------------------------------
' Include any libraries and external code files that we want to use
' ------------------------------------------------------------------------------

    #Include Once "tinyptc_ext.bi"
    #include once "windows.bi"

    ' Get ourselves set up for good coding practice
    OPTION STATIC
    OPTION EXPLICIT

' ------------------------------------------------------------------------------
' Set up any variables that will never change (they will stay constant)
' ------------------------------------------------------------------------------

    ' Our main display resolution
    CONST XRES=640
    CONST YRES=480

    CONST PI AS DOUBLE = 3.1415926535897932
    CONST PIR AS DOUBLE = PI/180

' ------------------------------------------------------------------------------
' We need to set up a bunch of variables that can be used anywhere in the code
' ------------------------------------------------------------------------------

    DIM SHARED AS UINTEGER DRAWBUFFER (XRES*YRES)

    DIM SHARED AS UINTEGER GLOBALROT
    DIM SHARED AS DOUBLE SINES(720), COSINES(720)
    DIM SHARED AS DOUBLE OLDTIMER

' ------------------------------------------------------------------------------
' Now to declare which subroutines we are going to write
' ------------------------------------------------------------------------------
    DECLARE FUNCTION InvSqrt(byval X AS SINGLE) AS SINGLE
    DECLARE SUB BLOBS()
    DECLARE SUB INIT()              :'  Set up our variables and arrays
    INIT()

' ------------------------------------------------------------------------------
' Now we're ready to get going, let's try and open a screen. If it fails, there
' is no point in continuing
' ------------------------------------------------------------------------------

'-------------------------------------------------------------------------------
'1st param. = Enable messagebox (dialog)
'2nd param. = Messagebox text
'3rd param. = Start in fullscreen or windowed mode if 1st param. is zero
'4th param. = Window normal border (0) or borderless (1)
    ptc_allowclose(1) 'enable tinyptc_ext to close on Escape Key press
    ptc_setdialog(1,"Fullscreen Mode?",1,1)
'-------------------------------------------------------------------------------

    IF ( PTC_OPEN ( "Blobby", XRES, YRES ) = 0 ) THEN
        END -1
    END IF

' ------------------------------------------------------------------------------
' We have our variables and arrays set up, we have an open screen display, so
' we can get going.
' ------------------------------------------------------------------------------

DO
    OLDTIMER=TIMER
    GLOBALROT=(GLOBALROT+1) MOD 360

    BLOBS()
    PTC_UPDATE@DRAWBUFFER(0)

    DO
        LOOP UNTIL TIMER>OLDTIMER+0.01
' Keep looping forever (or until the main escape handler catches us)
LOOP UNTIL 1=2

' We're all done now, we pressed escape, time to finish things up
END
PTC_CLOSE()

SUB BLOBS()
    DIM AS INTEGER X,Y,L,XX,YY,COL,YOFF,FULLPURPLE
    DIM AS SINGLE BX1,BX2,BX3,BY1,BY2,BY3,BR1,BR2,BR3
    DIM AS SINGLE SUM
    DIM AS UINTEGER PTR DP

    ' Once we get past a certain threshold, we'll stay purple
    ' so let's prepare a colour value for that
    FULLPURPLE=RGB(255,0,255)
   
    ' Now we'll set the position and radius of the three
    ' blobs according to the global rotation variable
    BX1=SINES(GLOBALROT+60)*60+320
    BX2=SINES(GLOBALROT+120)*90+320
    BX3=SINES(GLOBALROT+180)*120+320
    BY1=COSINES(GLOBALROT+60)*(SINES(GLOBALROT)*120)+280
    BY2=COSINES(GLOBALROT+120)*(SINES(GLOBALROT)*180)+280
    BY3=COSINES(GLOBALROT+180)*(SINES(GLOBALROT)*200)+280
    BR1=200
    BR2=250
    BR3=300

    ' loop through the whole screen, jumping by two.  This
    ' speeds the thing up without any real visible loss of
    ' detail.  To be perfect, you'd step one pixel at a time
    ' and only draw the single pixel in the drawing bit
    FOR Y=0 TO YRES-1 step 2
        YOFF=Y*XRES
        FOR X=0 TO XRES-1 step 2
            DP=@DRAWBUFFER(X+YOFF)
            XX=X-BX1
            YY=Y-BY1
' This is the original raw calculation which we replaced with the
' InvSqrt function.  As the InvSqrt returns 1/SQR, we multiply the
' result by the radius rather than dividing it.
'            SUM=BR1/SQR(XX*XX+YY*YY)
            SUM=BR1*InvSqrt(XX*XX+YY*YY)
            XX=X-BX2
            YY=Y-BY2
'            SUM+=BR2/SQR(XX*XX+YY*YY)
            SUM+=BR2*InvSqrt(XX*XX+YY*YY)
            XX=X-BX3
            YY=Y-BY3
'            SUM+=BR3/SQR(XX*XX+YY*YY)
            SUM+=BR3*InvSqrt(XX*XX+YY*YY)
            COL=SUM*40
            IF COL>255 THEN
                COL-=255
                IF COL>255 THEN COL=255
                COL=RGB(COL,0,255)
            END IF
            IF COL<=0 THEN COL=FULLPURPLE
            ' In this version, we're plotting four pixels as
            ' we're jumping through X and Y by two each time
           
            ' Place the colour in the pixel
            *DP=COL

            ' Move right one and plot the next
            ' pixel
            DP+=1
            *DP=COL
           
            ' Move down one line and left one pixel
            ' and plot the next pixel
            DP+=XRES-1
            *DP=COL
           
            ' Now move right one and plot again
            DP+=1
            *DP=COL
        NEXT X
    NEXT Y
END SUB

' This is the clever function that does a speedy 1/SQR calculation
' Copyright of this function is hard to track down, there are websites
' that are trying to find the source of who came up with the idea and
' the magic number.
'
' Suffice it to say that I didn't work this out, I just (with some
' help from Hellfire) converted it to Freebasic.
'
' It does some freaky stuff and you will get compiler warnings about
' suspicious pointer assignment.  Don't worry, it really *IS* doing
' some suspicious pointer stuff
FUNCTION InvSqrt(byval X AS SINGLE) AS SINGLE
    ' We want X to be positive
    X=ABS(X)
    ' It's important that we're working with integers and singles
    ' because of the clever stuff it does with parts of a computer
    ' stored variable.  You can surf the 'net for an explanation of
    ' what it does as I'm certainly not going to try and explain it
    DIM AS INTEGER PTR P
    DIM AS SINGLE PTR S
    DIM AS SINGLE xhalf
    DIM AS INTEGER i
    xhalf=0.5*X

    ' We want a pointer to the SINGLE variable X
    P=@X
    ' Now we store the appropriate parts of that SINGLE variable
    ' into the INTEGER variable i
    i=*P

    ' This is the wonderful magic number bit
    i=1597463007-(i SHR 1)

    ' We want a pointer to the INTEGER i
    S=@i
    ' So we can store the appropriate parts of that INTEGER variable
    ' into the SINGLE variable X
    X=*S
   
    ' Another quick calculation
    x=x*(1.5-xhalf*x*x)
    ' and throw the result back
    return x
END FUNCTION

' ------------------------------------------------------------------------------
' INIT
'
'   This is where we will initialise all of our variables and arrays.  It is
'   kept at the end of the code to keep the front tidy
' ------------------------------------------------------------------------------
SUB INIT()
    DIM AS INTEGER A

    FOR A=0 TO 719
        SINES(A)=SIN(A*PIR)
        COSINES(A)=COS(A*PIR)
    NEXT A

END SUB


[Jim]

I would question whether you need sqrt at all.  If you remove the sqrt and just use the squares and change the scale factor I suspect the fx will look almost identical.  You're getting quite a cheap divide though (the mul by 1/sqrt).

Jim

[Xalthorn]

I would question whether you need sqrt at all.  If you remove the sqrt and just use the squares and change the scale factor I suspect the fx will look almost identical.  You're getting quite a cheap divide though (the mul by 1/sqrt).

Jim

That's a thought.

I changed the calculations to

SUM=BR1/(XX*XX+YY*YY)

And increased the radius values to 16000, 20000, and 24000 respectively and it kinda works, but not as pretty.  I guess I would just need to find nicer radius values.  It seems to be a much more direct or linear scale than before.

It's damn fast though...