9 ms·

Reciprocal Approximation with 1 Subtraction

99 points by mbitsnbites 2y ago

Today's find: You can get a floating-point approximation of 1/x that's accurate to 3 bits with a single integer subtraction instruction. float fast_reciprocal(float x) { unsigned i = *(unsigned *) &x; i = 0x7effffffU - i; return *(float *) &i; } The magic number 0x7effffff accomplishes two things: 1) The exponent is calculated as 253-e, which effectively negates the exponent and subtracts 1. 2) The mantissa is approximated as a 1st order polynomial in the interval [1, 2). Interesting, but perhaps not very useful (as most CPU:s have more accurate reciprocal approximations these days).

magicalhippo 2y ago
Similar trick to that used in the fast inverse square root routine[1] popularized by Quake 3. [1]: https://en.wikipedia.org/wiki/Fast_inverse_square_root https://en.wikipedia.org/wiki/Fast_inverse_square_root
- mbitsnbites 2y ago
  Yep. Along the same lines. This one is even simpler, though, as it requires only a single integer CPU instruction (and the simplest of all instructions too). If you want full precision, you need to do three Newton-Raphson iterations after the initial approximation. One iteration is: y = y * (2.0F - x * y);
  - magicalhippo 2y ago
    It's a neat trick, and could be very useful on microcontrollers which doesn't have hardware division but does have hardware multiplication.
    - actionfromafar 2y ago
      Hm, like 68000?
      - pechay 2y ago
        The 68k has both signed and unsigned 32/16 bit divide instructions.
- fph 2y ago
  For some more explanation: the main idea behind both tricks is that the IEEE floating point formats are designed so that the most significant bits represent its exponent, that is, floor(log2(x)). Hence reinterpret-casting a float x to an unsigned integer uint(x) approximates a multiple of log2(x). So these kinds of approximations work like logarithm arithmetic: float(C + a*uint(x)) approximates x^a, for a suitable constant C. Quake's invsqrt is a=-1/2, this post is a=-1. More detail on https://en.wikipedia.org/wiki/Fast_inverse_square_root#Aliasing_to_an_integer_as_an_approximate_logarithm https://en.wikipedia.org/wiki/Fast_inverse_square_root#Alias... . IEEE754 floats are a very well-designed binary format, and the fact that these approximations are possible is part of this design; indeed, the first known instance of this trick is for a=1/2 by W. Kahan, the main designer of IEE754.
torusle 2y ago
There are couple of tricks you can do if you fiddle with the bits of a floating point value using integer arithmetic and binary logic. That was a thing back in the 90th.. I wonder how hard the performance hit from moving values between integer and float pipeline is nowadays. Last time I looked into that was the Cortex-A8 (first I-Phone area). Doing that kind of trick costed around 26 cycles (back and forth) due to pipeline stalls back then.
- stephencanon 2y ago
  There are basic integer operations in the FP/SIMD units on most CPUs, so there’s no generally need to “move back and forth” unless you need to branch on the result of a comparison, use a value as an address, or do some more specialized arithmetic.
  - stephencanon 2y ago
    (For that matter, though, most modern FP/SIMD units have a direct approximate-reciprocal instruction that is single-cycle throughput or better and much more accurate--generally around 10-12 bits, so there's no need for this sort of thing. See, e.g. FRECPE on ARM NEON and [V]RCPP[S/D] on x86.)
  - gpderetta 2y ago
    On x86 there is sometimes (depending on the specific microarchitecture) an extra cycle additional latency when using an integer operation on a xmm register last used with a float operation. I have seen it explained as the integer and foat ALUs's being physically distant and the forwarding network needing an extra cycle to transport the operands.
    - stephencanon 2y ago
      This is correct, but it’s happily pretty rare for it to matter in practice (because the domain bypass penalty is small and does not directly impact throughput, only latency).
- ack_complete 2y ago
Y_Y 2y ago
Quick comparison with exact and one Newton-Raphson: Value | True | Fast | Fast+Newton ---------------------------------------- 0.1 | 10.000 | 11.200 | 9.8560 0.5 | 2.0000 | 2.0000 | 2.0000 1.0 | 1.0000 | 1.0000 | 1.0000 2.0 | 0.5000 | 0.5000 | 0.5000 5.0 | 0.2000 | 0.2187 | 0.1982 10.0 | 0.1000 | 0.1094 | 0.0991 (where the extra correction was done with: y *= (2.0f - x*y); )
pkhuong 2y ago
You can do a little better if you tweak the low order bits (https://pvk.ca/Blog/LowLevel/software-reciprocal.html https://pvk.ca/Blog/LowLevel/software-reciprocal.html for the double float version)
jcmeyrignac 2y ago
A better value could be 0x7EEEEBB3, as in: https://github.com/parallella/pal/blob/bd9389f14fe16db4d963088a530a4f4908f42b01/src/math/p_inv.c#L20 https://github.com/parallella/pal/blob/bd9389f14fe16db4d9630...
- dahart 2y ago
  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…?
  - mbitsnbites 2y ago
    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.
    - mbitsnbites 2y ago
      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) Here, "good bits" is 24 minus the number of trailing non-zero-bits in the integer difference between the approximation and the correct value, looking at the IEEE 754 binary representation (if that makes sense). 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.
dahart 2y ago
This code is technically UB in C++, right? [1] Has anyone run into a case where it actually didn’t work? Just curious. I’ve often assumed that if C++ compilers didn’t compile this code, all hell would break loose. It might be nice to start sharing modern/safe versions of this snippet & the Quake thing. Is using memcpy the only option that is safe in both C and C++? That always felt really awkward to me. [1] https://tttapa.github.io/Pages/Programming/Cpp/Practices/type-punning.html https://tttapa.github.io/Pages/Programming/Cpp/Practices/typ...
- tekknolagi 2y ago
  We used memcpy everywhere in our runtime and after the 10th or so time doing it, it becomes less awkward.
  - dahart 2y ago
    And it’s always reliably optimized out in release builds, I assume?
    - LegionMammal978 2y ago
      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.
      - stouset 2y ago
        Sounds like a good place for a macro?
        mbitsnbites 2y ago
        We have memcpy behind a C++ template function that mimics the interface of std::bit_cast.
AlotOfReading 2y ago
A 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 https://news.ycombinator.com/item?id=42481612#42489596 In some cases, it can be faster than hardware division.
- dahart 2y ago
  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.
quanto 2y ago
A quick intuition: magic number 7e ffffff is negating by two's complement both the mantissa and exponent. 1. 7e: the first sig bit has to be zero to preserve the sign of the overall number. 2. ffffff: due to Taylor series 1/(1 + X) = 1 - X ..., negating gives the multiplicative inverse. Although an IEEE float has a sign bit (thus 2's complement does not work on the float itself) but mantissa and the exponent individually work with 2's complement for different reasons. The exponent has a biased range due to being chopped by half; where as the mantissa has a biased range due to its definition and constant offset. The exponent's 1 offset (7e instead of 7f) is a bit more difficult to see at first -- the devil is in the details, but 2's complement is the basic intuition.
somat 2y ago
Whenever 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 interface. where as the subtract function is less heavy being behind the operator interface. or are they trying to use a floating point accelerated log aka floating point subtract?
- brudgers 2y ago
  Engineering mathematical calculations always looks like "jank." Take a look at Plauger's The C Standard Library. It is full of "magic" numbers. But then again, representing irrational numbers in binary is inherently janky.
- dahart 2y ago
  Yes of course there’s a more graceful approach. You can use an actual divide instruction or routine instead of the bit-casting subtraction trick. It’s not about interfaces. The point is that a single subtract is simpler math, and (probably) faster than a divide, and it gets you surprisingly close to the right answer. The reason it works is subtracting two exponents is a subtraction in log space which is equivalent to a divide in linear space. People don’t use this trick in practice that much if at all, especially these days with GPUs that do reciprocals and square roots in hardware, it’s just an interesting thing to know about that helps solidify our understanding of floating point numbers and logarithms.
- mbitsnbites 2y ago
  If you're not constrained to software solutions you have a whole world of opportunities. E.g. if it's a graphics or neural net pipeline you can pour tricks like this (or better) onto it. If it's a CPU then you can add special instructions that do exponent manipulation and the likes.
- jasomill 2y ago
  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/ https://www.werkema.com/2021/10/14/my-knuth-check/
kkkqkqkqkqlqlql 2y ago
Where is the obligatory "WHAT THE FUCK" comment?
ncruces 2y ago
I collected this, and a few others here: https://github.com/ncruces/fastmath/blob/main/fast.go https://github.com/ncruces/fastmath/blob/main/fast.go It's in Go, but ports easily to other languages. See also: https://stackoverflow.com/questions/32042673/optimized-low-accuracy-approximation-to-rootnx-n https://stackoverflow.com/questions/32042673/optimized-low-a...
- mbitsnbites 2y ago
  Excellent! Will have a look.

99 points by mbitsnbites 2y ago

magicalhippo 2y ago

Similar trick to that used in the fast inverse square root routine[1] popularized by Quake 3. [1]: https://en.wikipedia.org/wiki/Fast_inverse_square_root https://en.wikipedia.org/wiki/Fast_inverse_square_root

mbitsnbites 2y ago
Yep. Along the same lines. This one is even simpler, though, as it requires only a single integer CPU instruction (and the simplest of all instructions too). If you want full precision, you need to do three Newton-Raphson iterations after the initial approximation. One iteration is: y = y * (2.0F - x * y);
- magicalhippo 2y ago
  It's a neat trick, and could be very useful on microcontrollers which doesn't have hardware division but does have hardware multiplication.
  - actionfromafar 2y ago
    Hm, like 68000?
    - pechay 2y ago
      The 68k has both signed and unsigned 32/16 bit divide instructions.
fph 2y ago
For some more explanation: the main idea behind both tricks is that the IEEE floating point formats are designed so that the most significant bits represent its exponent, that is, floor(log2(x)). Hence reinterpret-casting a float x to an unsigned integer uint(x) approximates a multiple of log2(x). So these kinds of approximations work like logarithm arithmetic: float(C + a*uint(x)) approximates x^a, for a suitable constant C. Quake's invsqrt is a=-1/2, this post is a=-1. More detail on https://en.wikipedia.org/wiki/Fast_inverse_square_root#Aliasing_to_an_integer_as_an_approximate_logarithm https://en.wikipedia.org/wiki/Fast_inverse_square_root#Alias... . IEEE754 floats are a very well-designed binary format, and the fact that these approximations are possible is part of this design; indeed, the first known instance of this trick is for a=1/2 by W. Kahan, the main designer of IEE754.

torusle 2y ago

There are couple of tricks you can do if you fiddle with the bits of a floating point value using integer arithmetic and binary logic. That was a thing back in the 90th.. I wonder how hard the performance hit from moving values between integer and float pipeline is nowadays. Last time I looked into that was the Cortex-A8 (first I-Phone area). Doing that kind of trick costed around 26 cycles (back and forth) due to pipeline stalls back then.

stephencanon 2y ago
There are basic integer operations in the FP/SIMD units on most CPUs, so there’s no generally need to “move back and forth” unless you need to branch on the result of a comparison, use a value as an address, or do some more specialized arithmetic.
- stephencanon 2y ago
  (For that matter, though, most modern FP/SIMD units have a direct approximate-reciprocal instruction that is single-cycle throughput or better and much more accurate--generally around 10-12 bits, so there's no need for this sort of thing. See, e.g. FRECPE on ARM NEON and [V]RCPP[S/D] on x86.)
- gpderetta 2y ago
  On x86 there is sometimes (depending on the specific microarchitecture) an extra cycle additional latency when using an integer operation on a xmm register last used with a float operation. I have seen it explained as the integer and foat ALUs's being physically distant and the forwarding network needing an extra cycle to transport the operands.
  - stephencanon 2y ago
    This is correct, but it’s happily pretty rare for it to matter in practice (because the domain bypass penalty is small and does not directly impact throughput, only latency).
ack_complete 2y ago

Y_Y 2y ago

Quick comparison with exact and one Newton-Raphson: Value | True | Fast | Fast+Newton ---------------------------------------- 0.1 | 10.000 | 11.200 | 9.8560 0.5 | 2.0000 | 2.0000 | 2.0000 1.0 | 1.0000 | 1.0000 | 1.0000 2.0 | 0.5000 | 0.5000 | 0.5000 5.0 | 0.2000 | 0.2187 | 0.1982 10.0 | 0.1000 | 0.1094 | 0.0991 (where the extra correction was done with: y *= (2.0f - x*y); )

pkhuong 2y ago

You can do a little better if you tweak the low order bits (https://pvk.ca/Blog/LowLevel/software-reciprocal.html https://pvk.ca/Blog/LowLevel/software-reciprocal.html for the double float version)

jcmeyrignac 2y ago

A better value could be 0x7EEEEBB3, as in: https://github.com/parallella/pal/blob/bd9389f14fe16db4d963088a530a4f4908f42b01/src/math/p_inv.c#L20 https://github.com/parallella/pal/blob/bd9389f14fe16db4d9630...

dahart 2y ago
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…?
- mbitsnbites 2y ago
  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.
  - mbitsnbites 2y ago
    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) Here, "good bits" is 24 minus the number of trailing non-zero-bits in the integer difference between the approximation and the correct value, looking at the IEEE 754 binary representation (if that makes sense). 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.

dahart 2y ago

This code is technically UB in C++, right? [1] Has anyone run into a case where it actually didn’t work? Just curious. I’ve often assumed that if C++ compilers didn’t compile this code, all hell would break loose. It might be nice to start sharing modern/safe versions of this snippet & the Quake thing. Is using memcpy the only option that is safe in both C and C++? That always felt really awkward to me. [1] https://tttapa.github.io/Pages/Programming/Cpp/Practices/type-punning.html https://tttapa.github.io/Pages/Programming/Cpp/Practices/typ...

tekknolagi 2y ago
We used memcpy everywhere in our runtime and after the 10th or so time doing it, it becomes less awkward.
- dahart 2y ago
  And it’s always reliably optimized out in release builds, I assume?
  - LegionMammal978 2y ago
    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.
    - stouset 2y ago
      Sounds like a good place for a macro?
      - mbitsnbites 2y ago
        We have memcpy behind a C++ template function that mimics the interface of std::bit_cast.

AlotOfReading 2y ago

A 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 https://news.ycombinator.com/item?id=42481612#42489596 In some cases, it can be faster than hardware division.

dahart 2y ago
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.

quanto 2y ago

A quick intuition: magic number 7e ffffff is negating by two's complement both the mantissa and exponent. 1. 7e: the first sig bit has to be zero to preserve the sign of the overall number. 2. ffffff: due to Taylor series 1/(1 + X) = 1 - X ..., negating gives the multiplicative inverse. Although an IEEE float has a sign bit (thus 2's complement does not work on the float itself) but mantissa and the exponent individually work with 2's complement for different reasons. The exponent has a biased range due to being chopped by half; where as the mantissa has a biased range due to its definition and constant offset. The exponent's 1 offset (7e instead of 7f) is a bit more difficult to see at first -- the devil is in the details, but 2's complement is the basic intuition.

somat 2y ago

Whenever 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 interface. where as the subtract function is less heavy being behind the operator interface. or are they trying to use a floating point accelerated log aka floating point subtract?

brudgers 2y ago
Engineering mathematical calculations always looks like "jank." Take a look at Plauger's The C Standard Library. It is full of "magic" numbers. But then again, representing irrational numbers in binary is inherently janky.
dahart 2y ago
Yes of course there’s a more graceful approach. You can use an actual divide instruction or routine instead of the bit-casting subtraction trick. It’s not about interfaces. The point is that a single subtract is simpler math, and (probably) faster than a divide, and it gets you surprisingly close to the right answer. The reason it works is subtracting two exponents is a subtraction in log space which is equivalent to a divide in linear space. People don’t use this trick in practice that much if at all, especially these days with GPUs that do reciprocals and square roots in hardware, it’s just an interesting thing to know about that helps solidify our understanding of floating point numbers and logarithms.
mbitsnbites 2y ago
If you're not constrained to software solutions you have a whole world of opportunities. E.g. if it's a graphics or neural net pipeline you can pour tricks like this (or better) onto it. If it's a CPU then you can add special instructions that do exponent manipulation and the likes.
jasomill 2y ago
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/ https://www.werkema.com/2021/10/14/my-knuth-check/

kkkqkqkqkqlqlql 2y ago

Where is the obligatory "WHAT THE FUCK" comment?

ncruces 2y ago

I collected this, and a few others here: https://github.com/ncruces/fastmath/blob/main/fast.go https://github.com/ncruces/fastmath/blob/main/fast.go It's in Go, but ports easily to other languages. See also: https://stackoverflow.com/questions/32042673/optimized-low-accuracy-approximation-to-rootnx-n https://stackoverflow.com/questions/32042673/optimized-low-a...

mbitsnbites 2y ago
Excellent! Will have a look.