5 ms·

Evolutionary algorithm to find Quake III's fast inverse square hack with Python

173 points by soravux 12y ago

Houshalter 12y ago
There is an amazing, easy to use peice of software for doing stuff like this called Eureqa (http://www.nutonian.com/products/eureqa/ http://www.nutonian.com/products/eureqa/). I've used it to find approximations for a lot of different functions.
- danbruc 12y ago
  Am I correct when I think of this as a curve/surface fitting tool that tries to automatically find a low complexity equation together with parameters that fits some provided data?
  - Houshalter 12y ago
    Yes, that's exactly it.
  - jmpeax 12y ago
    The term for this is called "Symbolic Regression".
- soravux 12y ago
  Eureqa seems a wonderful software for this specific application. If only I could get 2,500$ for a one-year license...
  - Houshalter 12y ago
    It has a 30 day free trial and you can use it for free if you limit to a few hundred datapoints (not ideal but workable.) It's not perfect for this application though because it doesn't bitwise operations or use execution time as a fitness function.
  - deleted 12y ago
    [deleted]
  - cfn 12y ago
    There are cheaper ones such as GeneXproTools.
- chengsun 12y ago
  On the flip side of the coin, there's ries, which takes a number and finds a mathematical expression that approximates it: http://mrob.com/pub/ries/ http://mrob.com/pub/ries/
  - Houshalter 12y ago
    That's really cool, I've been looking for something like that.
  - chinpokomon 12y ago
    I had a program for my HP-48GX that use to do something similar. It was especially useful when I knew that my decimal answer had some numerical constant, such as sqrt(2), and I wanted to convert it to something more reasonable. It would only return a single value, so I had to be careful that I didn't just use it every time, since just because an answer could approximate a multiple of Pi, it might only be a coincidence and not a good answer.
    - wtallis 12y ago
      Some of that was actually built-in to the calculator. The "→Q" and "→Q𝛑" functions took a floating-point number and return a rational number or rational expression in terms of pi. The tolerance for the approximation was determined by how many significant digits the calculator was configured to display. It was pretty nice for a calculator that didn't include a computer algebra system.
      - chinpokomon 12y ago
        It was specifically an improved library of those functions. I'll have to dig around on hpcalc.org and see if I can find it again. You could do algebraic calculations if you wrapped the statement in quotes, but once you got used to RPN that made little sense to do.
      - chinpokomon 12y ago
        It was specifically an improved library of those functions. I'll have to dig around on hpcalc.org and see if I can find it again. You could do algebraic calculations if you wrapped the statement in quotes, but once you got used to RPN that made little sense to do.
lelandbatey 12y ago
Is anyone else seeing a perpetual "loading" string where the example code is supposed to be?
- soravux 12y ago
  It's javascript that loads the Gist into the page. You may have NoScript or similar that prevents the loader to fetch the code. The first code is the one from the Wikipedia page (fast inverse square root) and the second is this one: https://gist.github.com/soravux/9673839 https://gist.github.com/soravux/9673839
s-macke 12y ago
This reminds of another trick on the C64 to multiply two 8 Bit numbers a*b = [ ((a+b)/2)^2 - ((a-b)/2)^2 ] more or less: one 8 bit additon, one 8 bit subtraction, two table lookups for the squares and one 16 bit subtraction.
- bane 12y ago
  Were the divide by twos just bit shift operations? (I have no idea if the 6510 had bit shift operators.)
  - s-macke 12y ago
    Well, I think it is part of the table. http://codebase64.org/doku.php?id=base:fast_8bit_multiplication_16bit_product http://codebase64.org/doku.php?id=base:fast_8bit_multiplicat... But the C64 had also a shifting operation.
  - muyuu 12y ago
    Yep, it does. ROL rotate left ROR rotate right ASL arithmetic shift left LSR logic shift right But for it to work for all 2-complement coded values we need to ensure that the MSB is copied back into itself to keep the value the same sign, which can be done with another shift. ; Divide a signed 16-bit value by 2 LDA MEM+1 ;Load the MSB ASL A ;Copy the sign bit into C ROR MEM+1 ;restore MSB ROR MEM+0 ;Rotate the LSB It's a bunch of cycles but nothing compared to a generic division.
Orangeair 12y ago
Why on earth would they create a named constant for 1.5F, but not for 0x5f3759df?
- jws 12y ago
  What would you name it? #define AND_THEN_A_MIRACLE_OCCURS (0x5f3759df)
  - officialjunk 12y ago
    MAGIC_NUMBER?
  - Orangeair 12y ago
    Fair enough. Still, it seems silly to name a constant THREE_HALVES. What information does that add over 1.5? What happens if you change it and it becomes wrong?
alexholehouse 12y ago
This is probably one of my favourite HN posts of 2014. Outstanding.
chengsun 12y ago
Řrřola (a famous demoscene coder [1]) found a far better constant than the original (a 2.7-fold accuracy improvement) by using an optimisation algorithm over all three constants present in the original inverse square hack. His algorithm and the result (along with a plot of relative error) can be found on his blog [2]. [1]: See, for example, this 256-byte raycaster https://www.youtube.com/watch?v=R35UuntQQF8 https://www.youtube.com/watch?v=R35UuntQQF8 [2]: http://rrrola.wz.cz/inv_sqrt.html http://rrrola.wz.cz/inv_sqrt.html
- soravux 12y ago
  This is really interesting. I'll be sure to include these references in my next blog post on the subject. The goal of the post was to find the equation from scratch, though. As it can be seen, the constant optimization was less the spotlight of the post.
  - zo1 12y ago
    Not to mention your blog engine/site is broken. Please fix the fact that we can't click on the right scrollbar due to your fancy shmancy javascript side menu. Put it on the left hand side, or don't overlay it over the scrollbar. Damn #hipstercoders
    - nitrogen 12y ago
      The placement of any site elements on top of the scrollbar is indeed annoying. Is that a Blogspot thing, or is it due to a customization?
      - masklinn 12y ago
        It's blogspots "Dynamic Views". Complete dreck. I can only assume it's the default when you create a blog nowadays as I have no idea why anybody (least of all developers) would opt into that crap.
        coldpie 12y ago
hcarvalhoalves 12y ago
Now the question is: could you arrive at this algorithm without previous knowledge of what it should converge to?
- pdenya 12y ago
  Why is that a relevant question now? How would doing this be more helpful without the results of the algorithm you're attempting to cut corners to approximate?
- soravux 12y ago
  That is the exact goal of the post. I arrived to the algorithm through evolution with random initialization. I enforced absolutely no heuristic to make it converge to this equation.
  - hcarvalhoalves 12y ago
    I know, but can we prove there's no heuristic on that algorithm? Because if the answer is "yes", this should be general enough to find other optimizations. It's fun to think about.
    - soravux 12y ago
      The best answer I can provide is: The code is there. There are no heuristics in it. Everytime you run it, you get different results (because of random initialization seeds and the stochastic nature of EA). It may find the (a - (x >> 1)) equation on a specific execution, or not. Over the runs I made, this equation (or similar) was the most popular and nothing come close to it. In fact, it finds a lot of other optimizations; either less accurate or way more complex. I remember getting tens and tens of operations in equations with "bof" accuracies.
acchow 12y ago
Oh, I thought this gonna be like regex golf: http://xkcd.com/1313/ http://xkcd.com/1313/ Misleading title.
pixelcort 12y ago
For a moment I thought this was a way to workaround any patents with this hack, by rediscovering the hack dynamically; but I could be remembering a different patented technique.
- Houshalter 12y ago
  That's actually a cool idea. I wonder if that would hold up.
  - personalcompute 12y ago
    Under US law it would not, as even genuine independent (re-)discovery is not a defense to infringement.
    - Houshalter 12y ago
      That's pretty ridiculous. The purpose of the law is to keep people from unfairly copying ideas, not rediscovering them. However this is different. It's not a person rediscovering it, but an algorithm which happens to do exactly the same thing but is provably not unique or patented. It's a different algorithm that just happens to converge on the same process. This is perfectly analogous to patents on real world inventions - there are plenty of cases of different machines that do the same thing. You can't patent the output of an algorithm, just the process by which it gets to that output.
      - polymatter 12y ago
        You can be liable for triple damages if the court decides you knew about the patent and infringed anyway. Normal damages if you did not. I believe independent invention is meant to be so rare as to be suspicious (despite many documented cases where this happened). The (supposed) purpose is to give an incentive to inventors to freely broadcasting their inventions, by giving them a monopoly on it so they can recoup their investment. Without that incentive, inventors would keep their inventions trade secrets.
    - 12y ago
anonymousUser 12y ago
Hacker's Delight [1][2] is an amazing book dedicated to discussing about all kinds of bit hacking, including fast integer square root, cube root calculation. [1]: http://www.hackersdelight.org/ http://www.hackersdelight.org/ [2]: http://www.amazon.com/Hackers-Delight-Edition-Henry-Warren/dp/0321842685 http://www.amazon.com/Hackers-Delight-Edition-Henry-Warren/d...
EGreg 12y ago
Shameless plug some of you might like to read: http://www.flipcode.com/tpractice/ http://www.flipcode.com/tpractice/ (I used to be really into game programming when I was in my teens.)
jheriko 12y ago
With fast native code its more than practical to search the entire space of floating point numbers and guarantee that you find the optimal magic number: http://jheriko-rtw.blogspot.co.uk/2009/04/understanding-and-improving-fast.html http://jheriko-rtw.blogspot.co.uk/2009/04/understanding-and-... As already mentioned in another comment it has been taken further by Řrřola who has optimised a version with the newton step including the constants involved in that. http://rrrola.wz.cz/inv_sqrt.html http://rrrola.wz.cz/inv_sqrt.html

173 points by soravux 12y ago

There is an amazing, easy to use peice of software for doing stuff like this called Eureqa (http://www.nutonian.com/products/eureqa/ http://www.nutonian.com/products/eureqa/). I've used it to find approximations for a lot of different functions.

danbruc 12y ago
Am I correct when I think of this as a curve/surface fitting tool that tries to automatically find a low complexity equation together with parameters that fits some provided data?
- Houshalter 12y ago
  Yes, that's exactly it.
- jmpeax 12y ago
  The term for this is called "Symbolic Regression".
soravux 12y ago
Eureqa seems a wonderful software for this specific application. If only I could get 2,500$ for a one-year license...
- Houshalter 12y ago
  It has a 30 day free trial and you can use it for free if you limit to a few hundred datapoints (not ideal but workable.) It's not perfect for this application though because it doesn't bitwise operations or use execution time as a fitness function.
- deleted 12y ago
  [deleted]
- cfn 12y ago
  There are cheaper ones such as GeneXproTools.
chengsun 12y ago
On the flip side of the coin, there's ries, which takes a number and finds a mathematical expression that approximates it: http://mrob.com/pub/ries/ http://mrob.com/pub/ries/
- Houshalter 12y ago
  That's really cool, I've been looking for something like that.
- chinpokomon 12y ago
  I had a program for my HP-48GX that use to do something similar. It was especially useful when I knew that my decimal answer had some numerical constant, such as sqrt(2), and I wanted to convert it to something more reasonable. It would only return a single value, so I had to be careful that I didn't just use it every time, since just because an answer could approximate a multiple of Pi, it might only be a coincidence and not a good answer.
  - wtallis 12y ago
    Some of that was actually built-in to the calculator. The "→Q" and "→Q𝛑" functions took a floating-point number and return a rational number or rational expression in terms of pi. The tolerance for the approximation was determined by how many significant digits the calculator was configured to display. It was pretty nice for a calculator that didn't include a computer algebra system.
    - chinpokomon 12y ago
      It was specifically an improved library of those functions. I'll have to dig around on hpcalc.org and see if I can find it again. You could do algebraic calculations if you wrapped the statement in quotes, but once you got used to RPN that made little sense to do.
    - chinpokomon 12y ago
      It was specifically an improved library of those functions. I'll have to dig around on hpcalc.org and see if I can find it again. You could do algebraic calculations if you wrapped the statement in quotes, but once you got used to RPN that made little sense to do.

lelandbatey 12y ago

Is anyone else seeing a perpetual "loading" string where the example code is supposed to be?

soravux 12y ago
It's javascript that loads the Gist into the page. You may have NoScript or similar that prevents the loader to fetch the code. The first code is the one from the Wikipedia page (fast inverse square root) and the second is this one: https://gist.github.com/soravux/9673839 https://gist.github.com/soravux/9673839

s-macke 12y ago

This reminds of another trick on the C64 to multiply two 8 Bit numbers a*b = [ ((a+b)/2)^2 - ((a-b)/2)^2 ] more or less: one 8 bit additon, one 8 bit subtraction, two table lookups for the squares and one 16 bit subtraction.

bane 12y ago
Were the divide by twos just bit shift operations? (I have no idea if the 6510 had bit shift operators.)
- s-macke 12y ago
  Well, I think it is part of the table. http://codebase64.org/doku.php?id=base:fast_8bit_multiplication_16bit_product http://codebase64.org/doku.php?id=base:fast_8bit_multiplicat... But the C64 had also a shifting operation.
- muyuu 12y ago
  Yep, it does. ROL rotate left ROR rotate right ASL arithmetic shift left LSR logic shift right But for it to work for all 2-complement coded values we need to ensure that the MSB is copied back into itself to keep the value the same sign, which can be done with another shift. ; Divide a signed 16-bit value by 2 LDA MEM+1 ;Load the MSB ASL A ;Copy the sign bit into C ROR MEM+1 ;restore MSB ROR MEM+0 ;Rotate the LSB It's a bunch of cycles but nothing compared to a generic division.

Orangeair 12y ago

Why on earth would they create a named constant for 1.5F, but not for 0x5f3759df?

jws 12y ago
What would you name it? #define AND_THEN_A_MIRACLE_OCCURS (0x5f3759df)
- officialjunk 12y ago
  MAGIC_NUMBER?
- Orangeair 12y ago
  Fair enough. Still, it seems silly to name a constant THREE_HALVES. What information does that add over 1.5? What happens if you change it and it becomes wrong?

alexholehouse 12y ago

This is probably one of my favourite HN posts of 2014. Outstanding.

chengsun 12y ago

Řrřola (a famous demoscene coder [1]) found a far better constant than the original (a 2.7-fold accuracy improvement) by using an optimisation algorithm over all three constants present in the original inverse square hack. His algorithm and the result (along with a plot of relative error) can be found on his blog [2]. [1]: See, for example, this 256-byte raycaster https://www.youtube.com/watch?v=R35UuntQQF8 https://www.youtube.com/watch?v=R35UuntQQF8 [2]: http://rrrola.wz.cz/inv_sqrt.html http://rrrola.wz.cz/inv_sqrt.html

soravux 12y ago
This is really interesting. I'll be sure to include these references in my next blog post on the subject. The goal of the post was to find the equation from scratch, though. As it can be seen, the constant optimization was less the spotlight of the post.
- zo1 12y ago
  Not to mention your blog engine/site is broken. Please fix the fact that we can't click on the right scrollbar due to your fancy shmancy javascript side menu. Put it on the left hand side, or don't overlay it over the scrollbar. Damn #hipstercoders
  - nitrogen 12y ago
    The placement of any site elements on top of the scrollbar is indeed annoying. Is that a Blogspot thing, or is it due to a customization?
    - masklinn 12y ago
      It's blogspots "Dynamic Views". Complete dreck. I can only assume it's the default when you create a blog nowadays as I have no idea why anybody (least of all developers) would opt into that crap.
      - coldpie 12y ago

hcarvalhoalves 12y ago

Now the question is: could you arrive at this algorithm without previous knowledge of what it should converge to?

pdenya 12y ago
Why is that a relevant question now? How would doing this be more helpful without the results of the algorithm you're attempting to cut corners to approximate?
soravux 12y ago
That is the exact goal of the post. I arrived to the algorithm through evolution with random initialization. I enforced absolutely no heuristic to make it converge to this equation.
- hcarvalhoalves 12y ago
  I know, but can we prove there's no heuristic on that algorithm? Because if the answer is "yes", this should be general enough to find other optimizations. It's fun to think about.
  - soravux 12y ago
    The best answer I can provide is: The code is there. There are no heuristics in it. Everytime you run it, you get different results (because of random initialization seeds and the stochastic nature of EA). It may find the (a - (x >> 1)) equation on a specific execution, or not. Over the runs I made, this equation (or similar) was the most popular and nothing come close to it. In fact, it finds a lot of other optimizations; either less accurate or way more complex. I remember getting tens and tens of operations in equations with "bof" accuracies.

acchow 12y ago

Oh, I thought this gonna be like regex golf: http://xkcd.com/1313/ http://xkcd.com/1313/ Misleading title.

pixelcort 12y ago

For a moment I thought this was a way to workaround any patents with this hack, by rediscovering the hack dynamically; but I could be remembering a different patented technique.

Houshalter 12y ago
That's actually a cool idea. I wonder if that would hold up.
- personalcompute 12y ago
  Under US law it would not, as even genuine independent (re-)discovery is not a defense to infringement.
  - Houshalter 12y ago
    That's pretty ridiculous. The purpose of the law is to keep people from unfairly copying ideas, not rediscovering them. However this is different. It's not a person rediscovering it, but an algorithm which happens to do exactly the same thing but is provably not unique or patented. It's a different algorithm that just happens to converge on the same process. This is perfectly analogous to patents on real world inventions - there are plenty of cases of different machines that do the same thing. You can't patent the output of an algorithm, just the process by which it gets to that output.
    - polymatter 12y ago
      You can be liable for triple damages if the court decides you knew about the patent and infringed anyway. Normal damages if you did not. I believe independent invention is meant to be so rare as to be suspicious (despite many documented cases where this happened). The (supposed) purpose is to give an incentive to inventors to freely broadcasting their inventions, by giving them a monopoly on it so they can recoup their investment. Without that incentive, inventors would keep their inventions trade secrets.
  - 12y ago

anonymousUser 12y ago

Hacker's Delight [1][2] is an amazing book dedicated to discussing about all kinds of bit hacking, including fast integer square root, cube root calculation. [1]: http://www.hackersdelight.org/ http://www.hackersdelight.org/ [2]: http://www.amazon.com/Hackers-Delight-Edition-Henry-Warren/dp/0321842685 http://www.amazon.com/Hackers-Delight-Edition-Henry-Warren/d...

EGreg 12y ago

Shameless plug some of you might like to read: http://www.flipcode.com/tpractice/ http://www.flipcode.com/tpractice/ (I used to be really into game programming when I was in my teens.)

jheriko 12y ago

With fast native code its more than practical to search the entire space of floating point numbers and guarantee that you find the optimal magic number: http://jheriko-rtw.blogspot.co.uk/2009/04/understanding-and-improving-fast.html http://jheriko-rtw.blogspot.co.uk/2009/04/understanding-and-... As already mentioned in another comment it has been taken further by Řrřola who has optimised a version with the newton step including the constants involved in that. http://rrrola.wz.cz/inv_sqrt.html http://rrrola.wz.cz/inv_sqrt.html