Earlier quoted context omitted.
Where did you find that note? I do not see it in the C standard draft I linked.
It is in n3301, the last 202y draft. In the 201x draft it was note 95.
I see it now. I guess the GCC developers were wrong then.
Reciprocal Approximation with 1 Subtraction
61–70 of 73 posts
Re: Reciprocal Approximation with 1 Subtraction
#62Earlier quoted context omitted.
It is in n3301, the last 202y draft. In the 201x draft it was note 95.
I see it now. I guess the GCC developers were wrong then.
I think it is a bit more complicated. The rule, together with the aliasing rule, if taken a face value, means you could do unrestricted aliasing as long as you cast to an union type on access. I believe that's the interpretation the GCC Devs reject as is makes TBAA ineffective.
Instead they interpret it narrowly to only allow punning through objects that are actual unions (as described in the GCC docs).
Unsurprisingly the standard is kind of a mess.
Re: Reciprocal Approximation with 1 Subtraction
#63A better value could be 0x7EEEEBB3, as in: https://github.com/parallella/pal/blob/bd9389f14fe16db4d9630...
EDIT: after double-checking my work, I realized I have a better bound on maximum error, but not a better average error. So, the magic number depends on the goal or metric, but mean relative error seems reasonable. Leaving my original comment here, but note the big caveat that I’m half wrong.
One can do better still - 0x7EF311BC is a near optimal value at least for inputs in the range of [0.001 .. 1000].
The simple explanation here is:
The post’s number 0x7EFFFFFF results in an approximation that is always equal to or greater than 1/x. The value 0x7EEEEBB3 is better, but it’s less than 1/x around 2/3rds of the time. My number 0x7EF311BC appears to be as well balanced as you can get, half the time greater and half the time less than 1/x.
To find this number, I have a Jupyter notebook that plots the maximum absolute value of relative error over a range of inputs, for a range of magic constants. Once it’s setup, it’s pretty easy to manually binary search and find the minimum. The plot of max error looks like a big “V”. (Edit while the plot of mean error looks like a big “U” near the minimum.
The optimal number does depend on the input range, and using a different range or allowing all finite floats will change where the optimal magic value is. The optimal magic number will also change if you add one or more Newton iterations, like in that github snippet (and also seen in the ‘quake trick’ code).
PPS maybe 0x7EF0F7D0 is a pretty good candidate for minimizing the average relative error…?
Re: Reciprocal Approximation with 1 Subtraction
#64A better value is 0x7eb504f3. You can follow that up with Newton-raphson to refine the approximation, or approximate the Newton-raphson too by multiplying 1.385356f*x. I did this a few days ago in the approximate division thread: https://news.ycombinator.com/item?id=42481612#42489596 In some cases, it can be faster than hardware division.
Not sure I understand your number 0x7eb504f3 - does it require using the 1.385 factor? Is it possible there is a typo in the number? That value doesn’t make sense to me. I’m measuring error here as the absolute value of relative error, e.g., | (v - v_approx) / v |. With that constant alone plugged into the poster’s code, I get a much less accurate approximation, with a minimum error of at least ~0.27 (meaning it’s always far away from the target) and worst case error of ~0.293 over the input range [1/1000…1000]. For comparison, the original 0x7effffff error is min:0 max:0.125. With 0x7eeeebb3, the error I get is min:0 max:0.0667. With 0x7ef311bc, error is min:0 max:0.0505. It’s important that the min error over a large interval is zero, because it means the approximation actually touches the target at least once.
Re: Reciprocal Approximation with 1 Subtraction
#65Whenever I see these tricks(see also: the quake 3 fast inverse sqrt) involving using, not casting, but using integers as floats directly and floats as integers, I wonder if there is a way to do it without the jank. Because what do you really want? some sort of exponent or exponent math right, some variant of the log function should work. is the problem is all the log functions are gated behind the function call inter…
For more general algorithms along these lines, see exercises 25–28 in section 1.2.2 of Knuth[1] (and note that the printed solution to exercise 28 in early printings has an error[2]).
[1] D.E. Knuth, The art of computer programming. Vol. 1, third edition, Addison-Wesley, Reading, MA, 1997.
Re: Reciprocal Approximation with 1 Subtraction
#66Whenever I see these tricks(see also: the quake 3 fast inverse sqrt) involving using, not casting, but using integers as floats directly and floats as integers, I wonder if there is a way to do it without the jank. Because what do you really want? some sort of exponent or exponent math right, some variant of the log function should work. is the problem is all the log functions are gated behind the function call inter…
For more general algorithms along these lines, see exercises 25–28 in section 1.2.2 of Knuth[1] (and note that the printed solution to exercise 28 in early printings has an error[2]). [1] D.E. Knuth, The art of computer programming. Vol. 1 , third edition, Addison-Wesley, Reading, MA, 1997. [2] https://www.werkema.com/2021/10/14/my-knuth-check/
That knuth check story was pretty great. it reminded me of when I had to write pow() from scratch many years ago. However mine was not optimized in any way, and was a very sluggish naive implementation of an algorithm I found in the nist mathematical function library.
I was young and worked night shifts feeding tapes to a ibm mainframe, in the down time I would amuse myself by reading the ibm manuals laying around and writing scripts on the operator console, some useful and some not so much.
One of my more useless scripts was a sort terminal characters drawing routine, I wanted to draw lines at a specific angle which would require sin(), cos(), however, the scripting language used had no included trig functions. So I looked it up. and smuggled in a printout of the nist reference sin function. At which point I found, to my dismay, there also was no floating point power function. So the next night I smuggled in the printout for that one as well. And now, armed with a highly questionable implementation of sin(), cos(), pow() I could finally draw the hands for my stupid analog clock screensaver on the 3270 terminal. Honestly probably the highlight of my time there, but I was too scared to let anyone know about that script, didn't want them to think I had too much downtime.
Re: Reciprocal Approximation with 1 Subtraction
#67A better value could be 0x7EEEEBB3, as in: https://github.com/parallella/pal/blob/bd9389f14fe16db4d9630...
EDIT: after double-checking my work, I realized I have a better bound on maximum error, but not a better average error. So, the magic number depends on the goal or metric, but mean relative error seems reasonable. Leaving my original comment here, but note the big caveat that I’m half wrong. One can do better still - 0x7EF311BC is a near optimal value at least for inputs in the range of [0.001 .. 1000]. The simple ex…
Your suggestion got me intrigued. I have a program that does an exhaustive check for maximum and average error, so I'll give your numbers a spin.
Re: Reciprocal Approximation with 1 Subtraction
#68Earlier quoted context omitted.
On platforms thar require aligned loads and stores (not x86 nor ARM), a direct pointer cast sometimes uses an aligned load/store where a memcpy uses multiple byte loads/stores, even on a good compiler, since memcpy() doesn't require that the pointers are aligned. This can be mitigated by going through a local variable, but it gets pretty verbose.
Sounds like a good place for a macro?
We have memcpy behind a C++ template function that mimics the interface of std::bit_cast.
Re: Reciprocal Approximation with 1 Subtraction
#69Earlier quoted context omitted.
These kinds of tricks are still used today. They're not so useful if you need a reciprocal or square root, since CPUs now have dedicated hardware for that, but it's different if you need a _cube_ root or x^(1/2.4).
I wonder to what extent the dedicated hardware is essentially implementing the same steps but at the transistor level.
The big cores do. They essentially pump division through something like an FMA (fused multiply-add) unit, possibly the same unit that is used for multiplication and addition. That's for the Newton-Raphson steps, or Goldschmidt steps.
In hardware it's much easier to do a LUT-based approximation for the initial estimate rather than the subtraction trick, though.
It's common for CPUs to give 6-8 accurate bits in the approximation. x86 gives 13 accurate bits. Back in 1975, the Cray 1 gave 30 (!) accurate bits in the first approximation, and it didn't even have a division instruction (everything about that machine was big and fast).
Re: Reciprocal Approximation with 1 Subtraction
#70Earlier quoted context omitted.
EDIT: after double-checking my work, I realized I have a better bound on maximum error, but not a better average error. So, the magic number depends on the goal or metric, but mean relative error seems reasonable. Leaving my original comment here, but note the big caveat that I’m half wrong. One can do better still - 0x7EF311BC is a near optimal value at least for inputs in the range of [0.001 .. 1000]. The simple ex…
Your suggestion got me intrigued. I have a program that does an exhaustive check for maximum and average error, so I'll give your numbers a spin.
Given my search criteria, the optimal magic number turns out to be: 0x7ef311c2
Initial approximation:
Good bits min: 4
Good bits avg: 5.242649912834
Error max: 0.0505102872849 (4.30728 bits)
Error avg: 0.0327344845327 (4.93304 bits)
1 NR step:
Good bits min: 8
Good bits avg: 10.642581939697
Error max: 0.00255139507338 (8.61450 bits)
Error avg: 0.00132373889641 (9.56117 bits)
2 NR steps:
Good bits min: 17
Good bits avg: 19.922843217850
Error max: 6.62494557693e-06 (17.20366 bits)
Error avg: 2.62858584054e-06 (18.53728 bits)
3 NR steps:
Good bits min: 23
Good bits avg: 23.674004554749
Error max: 1.19249960972e-07 (22.99951 bits)
Error avg: 3.44158509521e-08 (24.79235 bits)
Also, for the NR steps I used double precision for the inner (2.0 - x * y) part, then rounded to single precision, to simulate FMA, but single precision for the outer multiplication.