Calculating square roots.

[T3sl4co1l]

Nice.  Anyone got an equally simple ln function or e^x?

Yes.  Without using a table (potentially a large one, unwieldy for 8-bit MCUs), you can calculate the remaining bits by left-shifting the operand (N bits) until you have a 1.(N-1) fixed point number (the number of shifts is the integral part of the result N - Lg(operand)), then iterating the subtraction and squaring of the value.  Namely: each time you square the result, if the squaring generates a carry, then the result is greater than 2.  So we subtract 2 from the result (i.e., shift the carry into our output register as the next fractional bit of the log result) and continue.  After N passes, we could continue, but presumably, no more information can be had from it, so that the output has the same number of bits as the input.

Tim

[coppice]

That makes me strongly suspect someone will have done something like the software trick as a paper shortcut long ago, including their use to choose sensible starting approximations for interative algorithms. From the first floating point computers in the 1960s people were documenting all sorts of tricks to prime approximations based on the log nature that can be exposed by the interplay of integer and floating formats.

I'll still have to disagree somewhat. The main use for computers at that time was scientific/engineering/economics where a correct result generally trumped everything else. In a graphics renderer, you may not care about a 1%, 2% or even 5% error in extreme cases, because no one will complain. Especially if it means that the program can run at 60 FPS on someone's computer instead of 30 FPS. The idea of throwing accuracy out the window for raw speed, I argue, is younger than the '60s. In scientific calculations, you can never guarantee that the number won't be used in further calculations. In multimedia programming, you know that the display or speaker output is the final destination of the data, and the result is also not entirely critical. (Ie, the integrity of lives, buildings or economies don't depend on it.)

As for paper shortcuts, I doubt anyone would have used floating point tricks on paper, since binary floating is not a very natural notation system for humans. Of course, there are various other paper shortcuts.

There is no accuracy limitation. The graphics code is limited in accuracy, because they only run the iteration once or twice. The trick is to get a really good priming value for the loop, so it converges really quickly. That benefit is there for the crude graphics case, or with a few more iterations for the precision case.

You might want to read a little more about the kind of calculations computers did a century ago, when computer was a job title, and not a machine.

[ralphrmartin]

Here's my awesome inverse square root function

function invsqroot(x)
return x*x
end

Did I miss something?

hi school?

Inverse  square root is square. Maybe you missed something...

[miguelvp]

the negative value of the root degree makes that happen, not the inverse.

Edit: the term inverse in math is not always 1/ something in any event

[ralphrmartin]

the negative value of the root degree makes that happen, not the inverse.

Edit: the term inverse in math is not always 1/ something in any event

Inverse in advanced mathematics does not (always) mean 1/ something.

The inverse of a function f is a function, written f^-1, that undoes the effect of function f.

Trust me. I am a Chartered Mathematician as well as a Chartered Engineer. 

[miguelvp]

True 1/something as the inverse only applies to multiplication but all this is besides the point.

The function in question computes the (multiplication)inverse of the square root of a given number, so it's not x^2

Bringing semantics won't change that.

Edit: But when I first read your answer I thought it was obvious you said it in tongue in cheek.

[mikerj]

You might want to read a little more about the kind of calculations computers did a century ago, when computer was a job title, and not a machine.

You are still not seeing/understanding the principal operation that makes this function novel and fast, which was well described by Miguel above.  This is fundamentally a bit twiddling hack to reduce the number of iterations required, and is certainly not something you would do if you were manually using Newton Raphson or other methods for approximating roots.

[miguelvp]

I tracked down the original code (but using a different magic number) from:
Jim Blinn, Floating-point tricks, IEEE Comp. Graphics and Applications 17, no 4, 1997

Someone had a scan of the article but it's no longer available so Wayback Machine to the rescue:
https://web.archive.org/web/20130309073539/http://rufus.hackish.org/~rufus/FPtricks.pdf

Or get it from IEEE Xplore, but I haven't been an IEEE member for at least the last 15 years (I do have at least one paper published in there from 1994, however not in computer graphics, And I'm not sure I want to link it because it's really not relevant.) Anyways I'm not willing to pay $31 or even the $13 for members so I got it from the previous link.
http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=595279&queryText=jim+blinn&newsearch=true&searchField=Search_All

So we have to thanks Microsoft's Research for that Gem, even if there was an earlier but more obscure use
http://minnie.tuhs.org/cgi-bin/utree.pl?file=V5/usr/source/s3/sqrt.s (Dennis Ritchie 1974 as a support C library for the PDP-11)

Binn does attribute some ties to Graphic Gems 4 (which I do have but I donated all of them to work but I have access to them), but the actual Quake code as it's coded seems to have derived from that IEEE publication, and the magic number tweaked by Gary Tarolli (3DFX) probably to make it work better for the range they needed.

I will check on Monday to see if there is any similar code on Graphic Gems 4, I did look at Graphic Gems 1 on the past two days but didn't find anything relevant, but since Blinn cites volume 4 I bet there is something in there.

The constant on the IEEE paper is:
(OneAsInteger + (OneAsInteger>>1))

OneAsInteger is 0x3F800000 making the constant 0x5F400000 but as mentioned earlier that constant probably wasn't good enough for the range so it was tweaked.

For the meaning of: (OneAsInteger + (OneAsInteger>>1)) you have to read the paper and see what:
int i = (AsInteger(x)>>1)+(OneAsInteger>>1);

That gives you an approximation of the square of x, so it seems that adding 0x3F800000 as an interger performs the 1/x operation.

I guess I'll read that article in more detail because he (Jim Blinn) goes on to explain a lot of what is going on.

Edit: The article starts with this:

Well, I recently learned about some interesting properties
of the IEEE floating-point representation that sort
of mixes the two. To see (and believe) how this works,
we must first review this representation.

So the knowledge comes from somewhere else? Monday I'll check all of the Graphics Gems 1-4 (I have more but it has to be on the ones pre 1997)

[ralphrmartin]

OK, I'll stop joking around, and just note that the proper name for what someone called inverse square root is actually reciprocal square root.

[mikerj]

I tracked down the original code (but using a different magic number) from:
Jim Blinn, Floating-point tricks, IEEE Comp. Graphics and Applications 17, no 4, 1997

Someone had a scan of the article but it's no longer available so Wayback Machine to the rescue:
https://web.archive.org/web/20130309073539/http://rufus.hackish.org/~rufus/FPtricks.pdf

Thanks for posting that, it's a really interesting article.

[dom0]

If you want the history of the algorithm you need to go back much farther. I really don't know just how far, but Gauss is usually a good guess for who originally came up with things like that. :-)

Nope! Look closer. The real trick here is not the Newton-Raphson successive approximation which, as you mention, is centuries old. The real trick is the dirty casting from float to int to float, and doing an int operation on the float data. This requires a deeper understanding on how floats are stored bitwise, and saves one iteration in the successive approximation, which cuts the number of float multiplications in half.

I think the first publicly available occurrence of this specific method is somewhere in the original UNIX (for PDP 11) sources. Mid 70s.

[miguelvp]

I think the first publicly available occurrence of this specific method is somewhere in the original UNIX (for PDP 11) sources. Mid 70s.

As stated earlier but that used an unknown magic number stored at 0x4E84 and probably was done by Dennis Ritchie

...
So we have to thanks Microsoft's Research for that Gem, even if there was an earlier but more obscure use
http://minnie.tuhs.org/cgi-bin/utree.pl?file=V5/usr/source/s3/sqrt.s (Dennis Ritchie 1974 as a support C library for the PDP-11)
...