8 ms·

Faster asin() was hiding in plain sight

263 points by def-pri-pub 6mo ago

deleted 6mo ago
[deleted]
erichocean 6mo ago
Ideal HN content, thanks!
orangepanda 6mo ago
> Nobody likes throwing away work they've done I like throwing away work I've done. Frees up my mental capacity for other work to throw away.
patchnull 6mo ago
[flagged]
- stephencanon 6mo ago
  These sorts of approximations (and more sophisticated methods) are fairly widely used in systems programming, as seen by the fact that Apple's asin is only a couple percent slower and sub-ulp accurate (https://members.loria.fr/PZimmermann/papers/accuracy.pdf https://members.loria.fr/PZimmermann/papers/accuracy.pdf). I would expect to get similar performance on non-Apple x86 using Intel's math library, which does not seem to have been measured, and significantly better performance while preserving accuracy using a vectorized library call. The approximation reported here is slightly faster but only accurate to about 2.7e11 ulp. That's totally appropriate for the graphics use in question, but no one would ever use it for a system library; less than half the bits are good. Also worth noting that it's possible to go faster without further loss of accuracy--the approximation uses a correctly rounded square root, which is much more accurate than the rest of the approximation deserves. An approximate square root will deliver the same overall accuracy and much better vectorized performance.
  - Pannoniae 6mo ago
    Yeah, the only big problem with approx. sqrt is that it's not consistent across systems, for example Intel and AMD implement RSQRT differently... Fine for graphics, but if you need consistency, that messes things up.
    - stephencanon 6mo ago
      Newer rsqrt approximations (ARM NEON and SVE, and the AVX512F approximations on x86) make the behavior architectural so this is somewhat less of a problem (it still varies between _architectures_, however).
    - def-pri-pub 6mo ago
      Wait, what? Do you have a resource I could read up on about that? That is moderately concerning if your math isn't portable across chips.
      - stephencanon 6mo ago
        When Intel specced the rsqrt[ps]s and rcp[ps]s instructions ~30 years ago, they didn't fully specify their behavior. They just said their relative error is "smaller than 1.5 * 2⁻¹²," which someone thought was very clever because it gave them leeway to use tables or piecewise linear approximations or digit-by-digit computation or whatever was best suited to future processors. Since these are not IEEE 754 correctly-rounded operations, and there was (by definition) no software that currently used them, this was "fine". And mostly it has been OK, except for some cases like games or simulations that want to get bitwise identical results across HW, which (if they're lucky) just don't use these operations or (if they're unlucky) use them and have to handle mismatches somehow. Compilers never generate these operations implicitly unless you're compiling with some sort of fast-math flag, so you mostly only get to them by explicitly using an intrinsic, and in theory you know what you're signing up for if you do that. However, this did make them unusable for some scenarios where you would otherwise like to use them, so a bunch of graphics and scientific computing and math library developers said "please fully specify these operations next time" and now NEON/SVE and AVX512 have fully-specified reciprocal estimates,¹ which solves the problem unless you have to interoperate between x86 and ARM. ¹ e.g. Intel "specifies" theirs here: https://www.intel.com/content/www/us/en/developer/articles/code-sample/reference-implementations-for-ia-approximation-instructions-vrcp14-vrsqrt14-vrcp28-vrsqrt28-vexp2.html https://www.intel.com/content/www/us/en/developer/articles/c... ARM's is a little more readable: https://developer.arm.com/documentation/ddi0596/2021-03/Shared-Pseudocode/Shared-Functions?lang=en#impl-shared.RecipSqrtEstimate.2 https://developer.arm.com/documentation/ddi0596/2021-03/Shar...
- def-pri-pub 6mo ago
  I did scan some (major) open source games and graphics related project and found a few of them using `std::asin()`. I plan on submitting some patches.
adampunk 6mo ago
We love to leave faster functions languishing in library code. The basis for Q3A’s fast inverse square root had been sitting in fdlibm since 1986, on the net since 1993: https://www.netlib.org/fdlibm/e_sqrt.c https://www.netlib.org/fdlibm/e_sqrt.c
- def-pri-pub 6mo ago
  Funny enough that fdlimb implementation of asin() did come up in my research. I believe it might have been more performant in the past. But taking a quick scan of `e_asin.c`, I see it doing something similar to the Cg asin() implementation (but with more terms and more multiplications, which my guess is that it's slower). I think I see it also taking more branches (which could also lead to more of a slowdown).
  - adampunk 6mo ago
    Yeah Ng’s work in fdlibm is cool and really clever in parts but a lot of branching. Some of the ways they reach correct rounding are…so cool.
drsopp 6mo ago
Did some quick calculations, and at this precision, it seems a table lookup might be able to fit in the L1 cache depending on the CPU model.
- Pannoniae 6mo ago
  Microbenchmarks. A LUT will win many of them but you pessimise the rest of the code. So unless a significant (read: 20+%) portion of your code goes into the LUT, there isn't that much point to bother. For almost any pure calculation without I/O, it's better to do the arithmetic than to do memory access.
  - jcalvinowens 6mo ago
    Locality within the LUT matters too: if you know you're looking up identical or nearby-enough values to benefit from caching, an LUT can be more of a win. You only pay the cache cost for the portion you actually touch at runtime. I could imagine some graphics workloads tend compute asin() repeatedly with nearby input values. But I'd guess the locality isn't local enough to matter, only eight double precision floats fit in a cache line.
  - hrmtst93837 6mo ago
    [flagged]
- groundzeros2015 6mo ago
  I don’t want to fill up L1 for sin.
- jcalvinowens 6mo ago
  Surely the loss in precision of a 32KB LUT for double precision asin() would be unacceptable?
  - Jyaif 6mo ago
    By interpolating between values you can get excellent results with LUTs much smaller than 32KB. Will it be faster than the computation from op, that I don't know.
scottlamb 6mo ago
Isn't the faster approach SIMD [edit: or GPU]? A 1.05x to 1.90x speedup is great. A 16x speedup is better! They could be orthogonal improvements, but if I were prioritizing, I'd go for SIMD first. I searched for asin on Intel's intrinsics guide. They have a AVX-512 instrinsic `_mm512_asin_ps` but it says "sequence" rather than single-instruction. Presumably the actual sequence they use is in some header file somewhere, but I don't know off-hand where to look, so I don't know how it compares to a SIMDified version of `fast_asin_cg`. https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#text=asin&ig_expand=359 https://www.intel.com/content/www/us/en/docs/intrinsics-guid...
- deleted 6mo ago
  [deleted]
- TimorousBestie 6mo ago
  I don’t know much about raytracing but it’s probably tricky to orchestrate all those asin calls so that the input and output memory is aligned and contiguous. My uneducated intuition is that there’s little regularity as to which pixels will take which branches and will end up requiring which asin calls, but I might be wrong.
  - scottlamb 6mo ago
    I'd expect it to come down to data-oriented design: SoA (structure of arrays) rather than AoS (array of structures). I skimmed the author's source code, and this is where I'd start: https://github.com/define-private-public/PSRayTracing/blob/8dea5113f6b00e1ef6bb5c0c117562e971280271/render_library/Objects/HittableList.cpp#L39 https://github.com/define-private-public/PSRayTracing/blob/8... Instead of an `_objects`, I might try for a `_spheres`, `_boxes`, etc. (Or just `_lists` still using the virtual dispatch but for each list, rather than each object.) The `asin` seems to be used just for spheres. Within my `Spheres::closest_hit` (note plural), I'd work to SIMDify it. (I'd try to SIMDify the others too of course but apparently not with `asin`.) I think it's doable: https://github.com/define-private-public/PSRayTracing/blob/8dea5113f6b00e1ef6bb5c0c117562e971280271/render_library/Objects/Sphere.cpp#L34 https://github.com/define-private-public/PSRayTracing/blob/8... I don't know much about ray tracers either (having only written a super-naive one back in college) but this is the general technique used to speed up games, I believe. Besides enabling SIMD, it's more cache-efficient and minimizes dispatch overhead. edit: there's also stuff that you can hoist in this impl. Restructuring as SoA isn't strictly necessary to do that, but it might make it more obvious and natural. As an example, this `ray_dir.length_squared()` is the same for the whole list. You'd notice that when iterating over the spheres. https://github.com/define-private-public/PSRayTracing/blob/8dea5113f6b00e1ef6bb5c0c117562e971280271/render_library/Objects/Sphere.cpp#L43 https://github.com/define-private-public/PSRayTracing/blob/8...
AlotOfReading 6mo ago
I'm pretty sure it's not faster, but it was fun to write: float asin(float x) { float x2 = 1.0f-fabs(x); u32 i = bitcast(x2); i = 0x5f3759df - (i>>1); float inv = bitcast(i); return copysign(pi/2-pi/2*(x2*inv),x); } Courtesy of evil floating point bithacks.
- adampunk 6mo ago
  // what the fuck
- def-pri-pub 6mo ago
  > floating point bithacks The forbidden magic
- chuckadams 6mo ago
  You brought Zalgo. I blame this decade on you.
- moffkalast 6mo ago
  > float asinine(float x) { FTFY :P
- jacquesm 6mo ago
  That could do with some subtitles.
  - irishcoffee 6mo ago
    https://en.wikipedia.org/wiki/Fast_inverse_square_root https://en.wikipedia.org/wiki/Fast_inverse_square_root
- teo_zero 6mo ago
  The bad thing about this method is that it's slower than native CPU instructions. The good thing is that the result is very precise for at least 2 values of x, namely 1.0 and -1.0 JK
  - 6mo ago
LegionMammal978 6mo ago
In general, I find that minimax approximation is an underappreciated tool, especially the quite simple Remez algorithm to generate an optimal polynomial approximation [0]. With some modifications, you can adapt it to optimize for either absolute or relative error within an interval, or even come up with rational-function approximations. (Though unfortunately, many presentations of the algorithm use overly-simple forms of sample point selection that can break down on nontrivial input curves, especially if they contain small oscillations.) [0] https://en.wikipedia.org/wiki/Remez_algorithm https://en.wikipedia.org/wiki/Remez_algorithm
- jason_s 6mo ago
  Not sure I would call Remez "simple"... it's all relative; I prefer Chebyshev approximation which is simpler than Remez.
  - stephencanon 6mo ago
    Ideally either one is just a library call to generate the coefficients. Remez can get into trouble near the endpoints of the interval for asin and require a little bit of manual intervention, however.
  - LegionMammal978 6mo ago
    Perhaps, but at least I find it very simple for the optimality properties it gives: there is no inherent need to say, "I know that better parameters likely exist, but the algorithm to find them would be hopelessly expensive," as is the case in many forms of minimax optimization.
- herf 6mo ago
  They teach a lot of Taylor/Maclaurin series in Math classes (and trig functions are sometimes called "CORDIC" which is an old method too) but these are not used much in actual FPUs and libraries. Maybe we should update the curricula so people know better ways.
  - bee_rider 6mo ago
deleted 6mo ago
[deleted]
stephc_int13 6mo ago
My favorite tool to experiment with math approximation is lolremez. And you can easily ask your llm to do it for you.
glitchc 6mo ago
The 4% improvement doesn't seem like it's worth the effort. On a general note, instructions like division and square root are roughly equal to trig functions in cycle count on modern CPUs. So, replacing one with the other will not confer much benefit, as evidenced from the results. They're all typically implemented using LUTs, and it's hard to beat the performance of an optimized LUT, which is basically a multiplexer connected to some dedicated memory cells in hardware.
- kstrauser 6mo ago
  > The 4% improvement doesn't seem like it's worth the effort. People have gotten PhDs for smaller optimizations. I know. I've worked with them. > instructions like division and square root are roughly equal to trig functions in cycle count on modern CPUs. What's the x86-64 opcode for arcsin?
  - adrian_b 6mo ago
    Presumably the poster meant polynomial approximations of trigonometric functions not instructions for trigonometric functions, which are missing in most CPUs, though many GPUs have such instructions. x86-64 had instructions for the exponential and logarithmic functions in Xeon Phi, but those instructions have been removed in Skylake Server and the later Intel or AMD CPUs with AVX-512 support. However, instructions for trigonometric functions have no longer been added after Intel 80387, and those of 8087 and 80387 are deprecated.
  - glitchc 6mo ago
    > What's the x86-64 opcode for arcsin? Not required. ATAN and SQRTS(S|D) are sufficient, the half-angle approach in the article is the recommended way. > People have gotten PhDs for smaller optimizations. I know. I've worked with them. I understand the can, not sure about the should. Not trying to be snarky, we just seem to be producing PhDs with the slimmest of justifications. The bar needs to be higher.
    - kstrauser 6mo ago
      > I understand the can, not sure about the should. Not trying to be snarky, we just seem to be producing PhDs with the slimmest of justifications. The bar needs to be higher. I couldn't disagree more. Sure, making a 4% faster asin isn't going to change the world, but if it makes all callers a teensy bit faster, multiplied by the number of callers using it, then it adds up. Imagine the savings for a hyperscaler if they managed to made a more common instruction 4% faster.
jason_s 6mo ago
While I'm glad to see the OP got a good minimax solution at the end, it seems like the article missed clarifying one of the key points: error waveforms over a specified interval are critical, and if you don't see the characteristic minimax-like wiggle, you're wasting easy opportunity for improvement. Taylor series in general are a poor choice, and Pade approximants of Taylor series are equally poor. If you're going to use Pade approximants, they should be of the original function. I prefer Chebyshev approximation: https://www.embeddedrelated.com/showarticle/152.php https://www.embeddedrelated.com/showarticle/152.php which is often close enough to the more complicated Remez algorithm.
- ogogmad 6mo ago
  Chebyshev polynomials cos(n arcos(x)) provide one of the proofs that every continuous function f:[0,1]->R can be uniformly approximated by polynomial functions. Bernstein polynomials provide a shorter proof, but perhaps not the best numerical method: https://en.wikipedia.org/wiki/Bernstein_polynomial#See_also https://en.wikipedia.org/wiki/Bernstein_polynomial#See_also
  - AlotOfReading 6mo ago
    Those don't guarantee that that they can be well approximated by a polynomial of degree N though, like we have here. You can apply Jackson's inequality to calculate a maximum error bound, but the epsilon for degree 5 is pretty atrocious.
- rockmeamedee 6mo ago
  I had no idea, but this "wiggle" is required for an optimal approximation, it's called the "equioscillation property" [https://en.wikipedia.org/wiki/Equioscillation_theorem https://en.wikipedia.org/wiki/Equioscillation_theorem]. For a polynomial P (of degree n) to approximate a function F on the real numbers with minimal absolute error, the max error value of |P - F| needs to be hit multiple times, (n+2 times to be precise). You need to have the polynomial "wiggle" back and forth between the top of the error bound and the bottom. And even more surprisingly, this is a necessary _and sufficient_! condition for optimality. If you find a polynomial whose error alternates and it hits its max error bound n+2 times, you know that no other polynomial of degree n can do better, that is the best error bound you can get for degree n. Very cool!
exmadscientist 6mo ago
This line: > This amazing snippet of code was languishing in the docs of dead software, which in turn the original formula was scrawled away in a math textbook from the 60s. was kind of telling for me. I have some background in this sort of work (and long ago concluded that there was pretty much nothing you can do to improve on existing code, unless either you have some new specific hardware or domain constraint, or you're just looking for something quick-n-dirty for whatever reason, or are willing to invest research-paper levels of time and effort) and to think that someone would call Abramowitz and Stegun "a math textbook from the 60s" is kind of funny. It's got a similar level of importance to its field as Knuth's Art of Computer Programming or stuff like that. It's not an obscure text. Yeah, you might forget what all is in it if you don't use it often, but you'd go "oh, of course that would be in there, wouldn't it...."
- wolfi1 6mo ago
  Abramowitz/Stegun has been updated 2010 and resides now here: https://dlmf.nist.gov/ https://dlmf.nist.gov/
  - rerdavies 6mo ago
    Doesn't seem to be terribly up to date though. It seems to use almost exclusively taylor series, and seems to be completely uninterested in error analysis of any kind. Unless I'm missing something.
    - exmadscientist 6mo ago
      It's a general-purpose reference for mathematicians, not specifically for numerical analysis. Mathematicians are usually interested in the boring old power series centered at zero (Maclaurin series), so that's what gets prominence.
- def-pri-pub 6mo ago
  These are books that my uni courses never had me read. I'm a little shocked at times at how my degree program skimped on some of the more famous texts.
ok123456 6mo ago
Chebyshev approximation for asin is sum(2T_n(x) / (pi*n*n),n), the even terms are 0.
empiricus 6mo ago
Does anyone knows the resources for the algos used in the HW implementations of math functions? I mean the algos inside the CPUs and GPUs. How they make a tradeoff between transistor number, power consumption, cycles, which algos allow this.
- groos 6mo ago
  Here's one way to do it. https://en.wikipedia.org/wiki/CORDIC https://en.wikipedia.org/wiki/CORDIC
  - empiricus 6mo ago
    Thanks, but this seems to be optimized for the smallest number of gates, so it applies for simple microcontrollers and FPGA, and with limited precision. I was interested in actual state of the art used in modern CPUs and GPUs.
xt00 6mo ago
To be accurate, this is originally from Hastings 1955, Princeton "APPROXIMATIONS FOR DIGITAL COMPUTERS BY CECIL HASTINGS", page 159-163, there are actually multiple versions of the approximation with different constants used. So the original work was done with the goal of being performant for computers of the 1950's. Then the famous Abramowitz and Stegun guys put that in formula 4.4.45 with permission, then the nvidia CG library wrote some code that was based upon the formula, likely with some optimizations.
- rerdavies 6mo ago
  I ran this down, because I have a particular interest in vectorizable function approximations. Particular those that exploit bit-banging to handle range normalization. (Anyone have a good reference for that?) Regrettably, this is NOT from Hastings 1955. Hastings provides Taylor series and Chebyshev polynomial approximations. The OP's solution is a Pade approximation, which are not covered at all in Hastings.
  - xt00 6mo ago
    When you say "this is NOT from Hastings" I had to double check my post again -- I guess you are saying that the Pade approximation is not from Hastings, but the polynomial approximation that the OP referenced from nvidia from A&S and ultimately from Hastings, definitely is in Hastings on page 159 -- I think you were referring to the Pade approximation not being in Hastings, which appears to be true yes. In the article it is interesting that the OP tried taylor expansion and pade approximation, but not the fairly standard "welp lets just fit a Nth order polynomial to the arcsin" which is what Hastings did back in the day.
  - adampunk 6mo ago
    >Particular those that exploit bit-banging to handle range normalization https://userpages.cs.umbc.edu/phatak/645/supl/Ng-ArgReduction.pdf https://userpages.cs.umbc.edu/phatak/645/supl/Ng-ArgReductio... That's tiny but weird.
sixo 6mo ago
It appears that the real lesson here was to lean quite a bit more on theory than a programmer's usual roll-your-own heuristic would suggest. A fantastic amount of collective human thought has been dedicated to function approximations in the last century; Taylor methods are over 200 years old and unlikely to come close to state-of-the-art.
patchnull 6mo ago
[flagged]
- Sesse__ 6mo ago
  And similarly, entire generations of programmers were never taught Horner's scheme. You can see it in the article, where they write stuff like A * x * x * x * x * x * x + B * x * x * x * x + C * x * x + D (10 muls, 3 muladds) instead of the faster tmp = x * x; ((A * tmp + B) * tmp + C) * tmp + D (1 mul, 3 muladds)
  - zahlman 6mo ago
    Is this outside of what compilers can do nowadays? (Or do they refuse because it's floating-point?)
  - 3371 6mo ago
    Isn't that for... readability...?
  - arkmm 6mo ago
    Didn't know this technique had a name, but I would think a modern compiler could make this optimization on its own, no?
    - Sesse__ 6mo ago
      No, it's not equivalent for floating point, so a compiler won't do it unless you do -fassociative-math (or a superset, such as -ffast-math), in which case all correctness bets are off.
  - woadwarrior01 6mo ago
    The common subexpression elimination (CSE) pass in compilers takes care of that.
    - cmovq 6mo ago
      Compilers cannot do this optimization for floating point [1] unless you're compiling with -ffast-math. In general, don't rely on compilers to optimize floating point sub-expressions. [1]: https://godbolt.org/z/8bEjE9Wxx https://godbolt.org/z/8bEjE9Wxx
cmovq 6mo ago
> After all of the above work and that talk in mind, I decided to ask an LLM. Impressive that an LLM managed to produce the answer from a 7 year old stack overflow answer all on its own! [1] This would have been the first search result for “fast asin” before this article was published. [1]: https://stackoverflow.com/a/26030435 https://stackoverflow.com/a/26030435
- def-pri-pub 6mo ago
  I did see that, but isn't the vast majority of that page talking about acos() instead?
  - seanhunter 6mo ago
    That’s equivalent right? acos x = pi/2 - asin x So if you’ve got one that’s fast you have them both.
- deleted 6mo ago
  [deleted]
varispeed 6mo ago
If you are interested in such "tricks", you should check out the classic Hacker's Delight by Henry Warren
debo_ 6mo ago
https://bash-org-archive.com/?427792 https://bash-org-archive.com/?427792
deleted 6mo ago
[deleted]
andyjohnson0 6mo ago
Interesting article. A few years back I implemented a bunch of maths primitives, including trig functions, using Taylor sequences etc, to see how it was done. An interesting challenge, even at the elementary level I was working at. So this article got me wondering how much accuracy is needed before computing a series beats pre-computed lookup tables and interpolation. Anyone got any relevant experience to share? How much accuracy does ray tracing require?
peterabbitcook 6mo ago
I am curious, did you check how much your benchmarks moved (time and errors) if at all if you told the compiler to use —-use_fast_math or -ffast-math? There’s generally not a faster version of inverse trig functions to inline, but it might optimize some other stuff out. Unrelated to that, I’ve seen implementations (ie julia/base/special/trig) that use a “rational approximation” to asin, did you go down that road at any point?
veltas 6mo ago
Just a point that the constexpr/const use in that C++ code makes no difference to output, and is just noise really.
Skeime 6mo ago
Wouldn't it also be much better to evaluate the Taylor polynomials using Horner's method, instead? (Maybe C++ can do this automatically, but given that there might be rounding differences, it probably won't.)
WithinReason 6mo ago
> In any graphics application trigonometric functions are frequently used. Counterpoint from the man himself, "avoiding trigonometry": https://iquilezles.org/articles/noacos/ https://iquilezles.org/articles/noacos/
- djmips 6mo ago
  And further to that. https://fgiesen.wordpress.com/2010/10/21/finish-your-derivations-please/ https://fgiesen.wordpress.com/2010/10/21/finish-your-derivat...
coloneljelly 6mo ago
In DSP math, it is common to use Chebyshev polynomial approximation. You can get incredibly precise results within your required bounds.
djmips 6mo ago
To the blog poster: Robin Green is an excellent resource Faster Math Functions: https://basesandframes.wordpress.com/wp-content/uploads/2016/05/fast-math-functions_p1.pdf https://basesandframes.wordpress.com/wp-content/uploads/2016... https://basesandframes.wordpress.com/wp-content/uploads/2016/05/fast-math-functions_p2.pdf https://basesandframes.wordpress.com/wp-content/uploads/2016... Even faster math functions GDC 2020: https://www.gdcvault.com/play/1027337/Math-in-Game-Development-Summit https://www.gdcvault.com/play/1027337/Math-in-Game-Developme...
Ono-Sendai 6mo ago
Here's my fast acos, which I think can be converted to an asin: https://forwardscattering.org/post/66 https://forwardscattering.org/post/66
tyleo 6mo ago
I had to do an atan() on an slow embedded device once for an autonomous robot competition. Fastest impl I came up with was rounding and big switch statement.
kazinator 6mo ago
The glibc implementation already has tests for several ranges and hacks for them including Taylor series: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/e_asin.c https://github.com/lattera/glibc/blob/master/sysdeps/ieee754... The smallest range is |x| < 1.49011611938477e-8. In this case the routine just returns x, after calling some underflow-checking routine. So right there, if we neglect this detail in our own wrapper, we may be able to get a speedup, at the cost of sending very small values to the Tayor series. The next smallest range tested is |x| < 0.125, and after that |x| < 0.5 and 0.75. The authors are cetainly not missing any brilliant trick hiding in plain sight; they are doing a more assiduous job.
tzs 6mo ago
In that Padé approximant I think you can save a couple multiplications. As written it does this: n = 1 - 367/714 * x**2 d = 1 - 81/119 * x**2 + 183/4760 * x**4 return x * (n/d) That's got 7 multiplies (I'm counting divide as a multiply) and 3 additions. (I'm assuming the optimizer only computes x^2 once and computes x^2 by squaring x^2, and that all the constants are calculated at compile time). Replace n/d with 1/(d/n) and then replace d/n with q + r/n where q is the quotient of polynomial d divided by polynomial n and r is the remainder. This is the result: n = 1 - 367/714 * x**2 q = 1587627/1346890 - 549/7340 * x**2 r = -240737/1346890 return x / (q + r/n) That's got 5 multiplies and 3 additions.

263 points by def-pri-pub 6mo ago

deleted 6mo ago

[deleted]

erichocean 6mo ago

Ideal HN content, thanks!

orangepanda 6mo ago

> Nobody likes throwing away work they've done I like throwing away work I've done. Frees up my mental capacity for other work to throw away.

patchnull 6mo ago

[flagged]

stephencanon 6mo ago
These sorts of approximations (and more sophisticated methods) are fairly widely used in systems programming, as seen by the fact that Apple's asin is only a couple percent slower and sub-ulp accurate (https://members.loria.fr/PZimmermann/papers/accuracy.pdf https://members.loria.fr/PZimmermann/papers/accuracy.pdf). I would expect to get similar performance on non-Apple x86 using Intel's math library, which does not seem to have been measured, and significantly better performance while preserving accuracy using a vectorized library call. The approximation reported here is slightly faster but only accurate to about 2.7e11 ulp. That's totally appropriate for the graphics use in question, but no one would ever use it for a system library; less than half the bits are good. Also worth noting that it's possible to go faster without further loss of accuracy--the approximation uses a correctly rounded square root, which is much more accurate than the rest of the approximation deserves. An approximate square root will deliver the same overall accuracy and much better vectorized performance.
- Pannoniae 6mo ago
  Yeah, the only big problem with approx. sqrt is that it's not consistent across systems, for example Intel and AMD implement RSQRT differently... Fine for graphics, but if you need consistency, that messes things up.
  - stephencanon 6mo ago
    Newer rsqrt approximations (ARM NEON and SVE, and the AVX512F approximations on x86) make the behavior architectural so this is somewhat less of a problem (it still varies between _architectures_, however).
  - def-pri-pub 6mo ago
    Wait, what? Do you have a resource I could read up on about that? That is moderately concerning if your math isn't portable across chips.
    - stephencanon 6mo ago
      When Intel specced the rsqrt[ps]s and rcp[ps]s instructions ~30 years ago, they didn't fully specify their behavior. They just said their relative error is "smaller than 1.5 * 2⁻¹²," which someone thought was very clever because it gave them leeway to use tables or piecewise linear approximations or digit-by-digit computation or whatever was best suited to future processors. Since these are not IEEE 754 correctly-rounded operations, and there was (by definition) no software that currently used them, this was "fine". And mostly it has been OK, except for some cases like games or simulations that want to get bitwise identical results across HW, which (if they're lucky) just don't use these operations or (if they're unlucky) use them and have to handle mismatches somehow. Compilers never generate these operations implicitly unless you're compiling with some sort of fast-math flag, so you mostly only get to them by explicitly using an intrinsic, and in theory you know what you're signing up for if you do that. However, this did make them unusable for some scenarios where you would otherwise like to use them, so a bunch of graphics and scientific computing and math library developers said "please fully specify these operations next time" and now NEON/SVE and AVX512 have fully-specified reciprocal estimates,¹ which solves the problem unless you have to interoperate between x86 and ARM. ¹ e.g. Intel "specifies" theirs here: https://www.intel.com/content/www/us/en/developer/articles/code-sample/reference-implementations-for-ia-approximation-instructions-vrcp14-vrsqrt14-vrcp28-vrsqrt28-vexp2.html https://www.intel.com/content/www/us/en/developer/articles/c... ARM's is a little more readable: https://developer.arm.com/documentation/ddi0596/2021-03/Shared-Pseudocode/Shared-Functions?lang=en#impl-shared.RecipSqrtEstimate.2 https://developer.arm.com/documentation/ddi0596/2021-03/Shar...
def-pri-pub 6mo ago
I did scan some (major) open source games and graphics related project and found a few of them using `std::asin()`. I plan on submitting some patches.

adampunk 6mo ago

We love to leave faster functions languishing in library code. The basis for Q3A’s fast inverse square root had been sitting in fdlibm since 1986, on the net since 1993: https://www.netlib.org/fdlibm/e_sqrt.c https://www.netlib.org/fdlibm/e_sqrt.c

def-pri-pub 6mo ago
Funny enough that fdlimb implementation of asin() did come up in my research. I believe it might have been more performant in the past. But taking a quick scan of `e_asin.c`, I see it doing something similar to the Cg asin() implementation (but with more terms and more multiplications, which my guess is that it's slower). I think I see it also taking more branches (which could also lead to more of a slowdown).
- adampunk 6mo ago
  Yeah Ng’s work in fdlibm is cool and really clever in parts but a lot of branching. Some of the ways they reach correct rounding are…so cool.

drsopp 6mo ago

Did some quick calculations, and at this precision, it seems a table lookup might be able to fit in the L1 cache depending on the CPU model.

Pannoniae 6mo ago
Microbenchmarks. A LUT will win many of them but you pessimise the rest of the code. So unless a significant (read: 20+%) portion of your code goes into the LUT, there isn't that much point to bother. For almost any pure calculation without I/O, it's better to do the arithmetic than to do memory access.
- jcalvinowens 6mo ago
  Locality within the LUT matters too: if you know you're looking up identical or nearby-enough values to benefit from caching, an LUT can be more of a win. You only pay the cache cost for the portion you actually touch at runtime. I could imagine some graphics workloads tend compute asin() repeatedly with nearby input values. But I'd guess the locality isn't local enough to matter, only eight double precision floats fit in a cache line.
- hrmtst93837 6mo ago
  [flagged]
groundzeros2015 6mo ago
I don’t want to fill up L1 for sin.
jcalvinowens 6mo ago
Surely the loss in precision of a 32KB LUT for double precision asin() would be unacceptable?
- Jyaif 6mo ago
  By interpolating between values you can get excellent results with LUTs much smaller than 32KB. Will it be faster than the computation from op, that I don't know.

scottlamb 6mo ago

Isn't the faster approach SIMD [edit: or GPU]? A 1.05x to 1.90x speedup is great. A 16x speedup is better! They could be orthogonal improvements, but if I were prioritizing, I'd go for SIMD first. I searched for asin on Intel's intrinsics guide. They have a AVX-512 instrinsic `_mm512_asin_ps` but it says "sequence" rather than single-instruction. Presumably the actual sequence they use is in some header file somewhere, but I don't know off-hand where to look, so I don't know how it compares to a SIMDified version of `fast_asin_cg`. https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#text=asin&ig_expand=359 https://www.intel.com/content/www/us/en/docs/intrinsics-guid...

deleted 6mo ago
[deleted]
TimorousBestie 6mo ago
I don’t know much about raytracing but it’s probably tricky to orchestrate all those asin calls so that the input and output memory is aligned and contiguous. My uneducated intuition is that there’s little regularity as to which pixels will take which branches and will end up requiring which asin calls, but I might be wrong.
- scottlamb 6mo ago
  I'd expect it to come down to data-oriented design: SoA (structure of arrays) rather than AoS (array of structures). I skimmed the author's source code, and this is where I'd start: https://github.com/define-private-public/PSRayTracing/blob/8dea5113f6b00e1ef6bb5c0c117562e971280271/render_library/Objects/HittableList.cpp#L39 https://github.com/define-private-public/PSRayTracing/blob/8... Instead of an `_objects`, I might try for a `_spheres`, `_boxes`, etc. (Or just `_lists` still using the virtual dispatch but for each list, rather than each object.) The `asin` seems to be used just for spheres. Within my `Spheres::closest_hit` (note plural), I'd work to SIMDify it. (I'd try to SIMDify the others too of course but apparently not with `asin`.) I think it's doable: https://github.com/define-private-public/PSRayTracing/blob/8dea5113f6b00e1ef6bb5c0c117562e971280271/render_library/Objects/Sphere.cpp#L34 https://github.com/define-private-public/PSRayTracing/blob/8... I don't know much about ray tracers either (having only written a super-naive one back in college) but this is the general technique used to speed up games, I believe. Besides enabling SIMD, it's more cache-efficient and minimizes dispatch overhead. edit: there's also stuff that you can hoist in this impl. Restructuring as SoA isn't strictly necessary to do that, but it might make it more obvious and natural. As an example, this `ray_dir.length_squared()` is the same for the whole list. You'd notice that when iterating over the spheres. https://github.com/define-private-public/PSRayTracing/blob/8dea5113f6b00e1ef6bb5c0c117562e971280271/render_library/Objects/Sphere.cpp#L43 https://github.com/define-private-public/PSRayTracing/blob/8...

AlotOfReading 6mo ago

I'm pretty sure it's not faster, but it was fun to write: float asin(float x) { float x2 = 1.0f-fabs(x); u32 i = bitcast(x2); i = 0x5f3759df - (i>>1); float inv = bitcast(i); return copysign(pi/2-pi/2*(x2*inv),x); } Courtesy of evil floating point bithacks.

adampunk 6mo ago
// what the fuck
def-pri-pub 6mo ago
> floating point bithacks The forbidden magic
chuckadams 6mo ago
You brought Zalgo. I blame this decade on you.
moffkalast 6mo ago
> float asinine(float x) { FTFY :P
jacquesm 6mo ago
That could do with some subtitles.
- irishcoffee 6mo ago
  https://en.wikipedia.org/wiki/Fast_inverse_square_root https://en.wikipedia.org/wiki/Fast_inverse_square_root
teo_zero 6mo ago
The bad thing about this method is that it's slower than native CPU instructions. The good thing is that the result is very precise for at least 2 values of x, namely 1.0 and -1.0 JK
- 6mo ago

LegionMammal978 6mo ago

In general, I find that minimax approximation is an underappreciated tool, especially the quite simple Remez algorithm to generate an optimal polynomial approximation [0]. With some modifications, you can adapt it to optimize for either absolute or relative error within an interval, or even come up with rational-function approximations. (Though unfortunately, many presentations of the algorithm use overly-simple forms of sample point selection that can break down on nontrivial input curves, especially if they contain small oscillations.) [0] https://en.wikipedia.org/wiki/Remez_algorithm https://en.wikipedia.org/wiki/Remez_algorithm

jason_s 6mo ago
Not sure I would call Remez "simple"... it's all relative; I prefer Chebyshev approximation which is simpler than Remez.
- stephencanon 6mo ago
  Ideally either one is just a library call to generate the coefficients. Remez can get into trouble near the endpoints of the interval for asin and require a little bit of manual intervention, however.
- LegionMammal978 6mo ago
  Perhaps, but at least I find it very simple for the optimality properties it gives: there is no inherent need to say, "I know that better parameters likely exist, but the algorithm to find them would be hopelessly expensive," as is the case in many forms of minimax optimization.
herf 6mo ago
They teach a lot of Taylor/Maclaurin series in Math classes (and trig functions are sometimes called "CORDIC" which is an old method too) but these are not used much in actual FPUs and libraries. Maybe we should update the curricula so people know better ways.
- bee_rider 6mo ago

deleted 6mo ago

[deleted]

stephc_int13 6mo ago

My favorite tool to experiment with math approximation is lolremez. And you can easily ask your llm to do it for you.

glitchc 6mo ago

The 4% improvement doesn't seem like it's worth the effort. On a general note, instructions like division and square root are roughly equal to trig functions in cycle count on modern CPUs. So, replacing one with the other will not confer much benefit, as evidenced from the results. They're all typically implemented using LUTs, and it's hard to beat the performance of an optimized LUT, which is basically a multiplexer connected to some dedicated memory cells in hardware.

kstrauser 6mo ago
> The 4% improvement doesn't seem like it's worth the effort. People have gotten PhDs for smaller optimizations. I know. I've worked with them. > instructions like division and square root are roughly equal to trig functions in cycle count on modern CPUs. What's the x86-64 opcode for arcsin?
- adrian_b 6mo ago
  Presumably the poster meant polynomial approximations of trigonometric functions not instructions for trigonometric functions, which are missing in most CPUs, though many GPUs have such instructions. x86-64 had instructions for the exponential and logarithmic functions in Xeon Phi, but those instructions have been removed in Skylake Server and the later Intel or AMD CPUs with AVX-512 support. However, instructions for trigonometric functions have no longer been added after Intel 80387, and those of 8087 and 80387 are deprecated.
- glitchc 6mo ago
  > What's the x86-64 opcode for arcsin? Not required. ATAN and SQRTS(S|D) are sufficient, the half-angle approach in the article is the recommended way. > People have gotten PhDs for smaller optimizations. I know. I've worked with them. I understand the can, not sure about the should. Not trying to be snarky, we just seem to be producing PhDs with the slimmest of justifications. The bar needs to be higher.
  - kstrauser 6mo ago
    > I understand the can, not sure about the should. Not trying to be snarky, we just seem to be producing PhDs with the slimmest of justifications. The bar needs to be higher. I couldn't disagree more. Sure, making a 4% faster asin isn't going to change the world, but if it makes all callers a teensy bit faster, multiplied by the number of callers using it, then it adds up. Imagine the savings for a hyperscaler if they managed to made a more common instruction 4% faster.

jason_s 6mo ago

While I'm glad to see the OP got a good minimax solution at the end, it seems like the article missed clarifying one of the key points: error waveforms over a specified interval are critical, and if you don't see the characteristic minimax-like wiggle, you're wasting easy opportunity for improvement. Taylor series in general are a poor choice, and Pade approximants of Taylor series are equally poor. If you're going to use Pade approximants, they should be of the original function. I prefer Chebyshev approximation: https://www.embeddedrelated.com/showarticle/152.php https://www.embeddedrelated.com/showarticle/152.php which is often close enough to the more complicated Remez algorithm.

ogogmad 6mo ago
Chebyshev polynomials cos(n arcos(x)) provide one of the proofs that every continuous function f:[0,1]->R can be uniformly approximated by polynomial functions. Bernstein polynomials provide a shorter proof, but perhaps not the best numerical method: https://en.wikipedia.org/wiki/Bernstein_polynomial#See_also https://en.wikipedia.org/wiki/Bernstein_polynomial#See_also
- AlotOfReading 6mo ago
  Those don't guarantee that that they can be well approximated by a polynomial of degree N though, like we have here. You can apply Jackson's inequality to calculate a maximum error bound, but the epsilon for degree 5 is pretty atrocious.
rockmeamedee 6mo ago
I had no idea, but this "wiggle" is required for an optimal approximation, it's called the "equioscillation property" [https://en.wikipedia.org/wiki/Equioscillation_theorem https://en.wikipedia.org/wiki/Equioscillation_theorem]. For a polynomial P (of degree n) to approximate a function F on the real numbers with minimal absolute error, the max error value of |P - F| needs to be hit multiple times, (n+2 times to be precise). You need to have the polynomial "wiggle" back and forth between the top of the error bound and the bottom. And even more surprisingly, this is a necessary _and sufficient_! condition for optimality. If you find a polynomial whose error alternates and it hits its max error bound n+2 times, you know that no other polynomial of degree n can do better, that is the best error bound you can get for degree n. Very cool!

exmadscientist 6mo ago

This line: > This amazing snippet of code was languishing in the docs of dead software, which in turn the original formula was scrawled away in a math textbook from the 60s. was kind of telling for me. I have some background in this sort of work (and long ago concluded that there was pretty much nothing you can do to improve on existing code, unless either you have some new specific hardware or domain constraint, or you're just looking for something quick-n-dirty for whatever reason, or are willing to invest research-paper levels of time and effort) and to think that someone would call Abramowitz and Stegun "a math textbook from the 60s" is kind of funny. It's got a similar level of importance to its field as Knuth's Art of Computer Programming or stuff like that. It's not an obscure text. Yeah, you might forget what all is in it if you don't use it often, but you'd go "oh, of course that would be in there, wouldn't it...."

wolfi1 6mo ago
Abramowitz/Stegun has been updated 2010 and resides now here: https://dlmf.nist.gov/ https://dlmf.nist.gov/
- rerdavies 6mo ago
  Doesn't seem to be terribly up to date though. It seems to use almost exclusively taylor series, and seems to be completely uninterested in error analysis of any kind. Unless I'm missing something.
  - exmadscientist 6mo ago
    It's a general-purpose reference for mathematicians, not specifically for numerical analysis. Mathematicians are usually interested in the boring old power series centered at zero (Maclaurin series), so that's what gets prominence.
def-pri-pub 6mo ago
These are books that my uni courses never had me read. I'm a little shocked at times at how my degree program skimped on some of the more famous texts.

ok123456 6mo ago

Chebyshev approximation for asin is sum(2T_n(x) / (pi*n*n),n), the even terms are 0.

empiricus 6mo ago

Does anyone knows the resources for the algos used in the HW implementations of math functions? I mean the algos inside the CPUs and GPUs. How they make a tradeoff between transistor number, power consumption, cycles, which algos allow this.

groos 6mo ago
Here's one way to do it. https://en.wikipedia.org/wiki/CORDIC https://en.wikipedia.org/wiki/CORDIC
- empiricus 6mo ago
  Thanks, but this seems to be optimized for the smallest number of gates, so it applies for simple microcontrollers and FPGA, and with limited precision. I was interested in actual state of the art used in modern CPUs and GPUs.

xt00 6mo ago

To be accurate, this is originally from Hastings 1955, Princeton "APPROXIMATIONS FOR DIGITAL COMPUTERS BY CECIL HASTINGS", page 159-163, there are actually multiple versions of the approximation with different constants used. So the original work was done with the goal of being performant for computers of the 1950's. Then the famous Abramowitz and Stegun guys put that in formula 4.4.45 with permission, then the nvidia CG library wrote some code that was based upon the formula, likely with some optimizations.

rerdavies 6mo ago
I ran this down, because I have a particular interest in vectorizable function approximations. Particular those that exploit bit-banging to handle range normalization. (Anyone have a good reference for that?) Regrettably, this is NOT from Hastings 1955. Hastings provides Taylor series and Chebyshev polynomial approximations. The OP's solution is a Pade approximation, which are not covered at all in Hastings.
- xt00 6mo ago
  When you say "this is NOT from Hastings" I had to double check my post again -- I guess you are saying that the Pade approximation is not from Hastings, but the polynomial approximation that the OP referenced from nvidia from A&S and ultimately from Hastings, definitely is in Hastings on page 159 -- I think you were referring to the Pade approximation not being in Hastings, which appears to be true yes. In the article it is interesting that the OP tried taylor expansion and pade approximation, but not the fairly standard "welp lets just fit a Nth order polynomial to the arcsin" which is what Hastings did back in the day.
- adampunk 6mo ago
  >Particular those that exploit bit-banging to handle range normalization https://userpages.cs.umbc.edu/phatak/645/supl/Ng-ArgReduction.pdf https://userpages.cs.umbc.edu/phatak/645/supl/Ng-ArgReductio... That's tiny but weird.

sixo 6mo ago

It appears that the real lesson here was to lean quite a bit more on theory than a programmer's usual roll-your-own heuristic would suggest. A fantastic amount of collective human thought has been dedicated to function approximations in the last century; Taylor methods are over 200 years old and unlikely to come close to state-of-the-art.

patchnull 6mo ago

[flagged]

Sesse__ 6mo ago
And similarly, entire generations of programmers were never taught Horner's scheme. You can see it in the article, where they write stuff like A * x * x * x * x * x * x + B * x * x * x * x + C * x * x + D (10 muls, 3 muladds) instead of the faster tmp = x * x; ((A * tmp + B) * tmp + C) * tmp + D (1 mul, 3 muladds)
- zahlman 6mo ago
  Is this outside of what compilers can do nowadays? (Or do they refuse because it's floating-point?)
- 3371 6mo ago
  Isn't that for... readability...?
- arkmm 6mo ago
  Didn't know this technique had a name, but I would think a modern compiler could make this optimization on its own, no?
  - Sesse__ 6mo ago
    No, it's not equivalent for floating point, so a compiler won't do it unless you do -fassociative-math (or a superset, such as -ffast-math), in which case all correctness bets are off.
- woadwarrior01 6mo ago
  The common subexpression elimination (CSE) pass in compilers takes care of that.
  - cmovq 6mo ago
    Compilers cannot do this optimization for floating point [1] unless you're compiling with -ffast-math. In general, don't rely on compilers to optimize floating point sub-expressions. [1]: https://godbolt.org/z/8bEjE9Wxx https://godbolt.org/z/8bEjE9Wxx

cmovq 6mo ago

> After all of the above work and that talk in mind, I decided to ask an LLM. Impressive that an LLM managed to produce the answer from a 7 year old stack overflow answer all on its own! [1] This would have been the first search result for “fast asin” before this article was published. [1]: https://stackoverflow.com/a/26030435 https://stackoverflow.com/a/26030435

def-pri-pub 6mo ago
I did see that, but isn't the vast majority of that page talking about acos() instead?
- seanhunter 6mo ago
  That’s equivalent right? acos x = pi/2 - asin x So if you’ve got one that’s fast you have them both.
deleted 6mo ago
[deleted]

varispeed 6mo ago

If you are interested in such "tricks", you should check out the classic Hacker's Delight by Henry Warren

debo_ 6mo ago

https://bash-org-archive.com/?427792 https://bash-org-archive.com/?427792

deleted 6mo ago

[deleted]

andyjohnson0 6mo ago

Interesting article. A few years back I implemented a bunch of maths primitives, including trig functions, using Taylor sequences etc, to see how it was done. An interesting challenge, even at the elementary level I was working at. So this article got me wondering how much accuracy is needed before computing a series beats pre-computed lookup tables and interpolation. Anyone got any relevant experience to share? How much accuracy does ray tracing require?

peterabbitcook 6mo ago

I am curious, did you check how much your benchmarks moved (time and errors) if at all if you told the compiler to use —-use_fast_math or -ffast-math? There’s generally not a faster version of inverse trig functions to inline, but it might optimize some other stuff out. Unrelated to that, I’ve seen implementations (ie julia/base/special/trig) that use a “rational approximation” to asin, did you go down that road at any point?

veltas 6mo ago

Just a point that the constexpr/const use in that C++ code makes no difference to output, and is just noise really.

Skeime 6mo ago

Wouldn't it also be much better to evaluate the Taylor polynomials using Horner's method, instead? (Maybe C++ can do this automatically, but given that there might be rounding differences, it probably won't.)

WithinReason 6mo ago

> In any graphics application trigonometric functions are frequently used. Counterpoint from the man himself, "avoiding trigonometry": https://iquilezles.org/articles/noacos/ https://iquilezles.org/articles/noacos/

djmips 6mo ago
And further to that. https://fgiesen.wordpress.com/2010/10/21/finish-your-derivations-please/ https://fgiesen.wordpress.com/2010/10/21/finish-your-derivat...

coloneljelly 6mo ago

In DSP math, it is common to use Chebyshev polynomial approximation. You can get incredibly precise results within your required bounds.

djmips 6mo ago

To the blog poster: Robin Green is an excellent resource Faster Math Functions: https://basesandframes.wordpress.com/wp-content/uploads/2016/05/fast-math-functions_p1.pdf https://basesandframes.wordpress.com/wp-content/uploads/2016... https://basesandframes.wordpress.com/wp-content/uploads/2016/05/fast-math-functions_p2.pdf https://basesandframes.wordpress.com/wp-content/uploads/2016... Even faster math functions GDC 2020: https://www.gdcvault.com/play/1027337/Math-in-Game-Development-Summit https://www.gdcvault.com/play/1027337/Math-in-Game-Developme...

Ono-Sendai 6mo ago

Here's my fast acos, which I think can be converted to an asin: https://forwardscattering.org/post/66 https://forwardscattering.org/post/66

tyleo 6mo ago

I had to do an atan() on an slow embedded device once for an autonomous robot competition. Fastest impl I came up with was rounding and big switch statement.

kazinator 6mo ago

The glibc implementation already has tests for several ranges and hacks for them including Taylor series: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/e_asin.c https://github.com/lattera/glibc/blob/master/sysdeps/ieee754... The smallest range is |x| < 1.49011611938477e-8. In this case the routine just returns x, after calling some underflow-checking routine. So right there, if we neglect this detail in our own wrapper, we may be able to get a speedup, at the cost of sending very small values to the Tayor series. The next smallest range tested is |x| < 0.125, and after that |x| < 0.5 and 0.75. The authors are cetainly not missing any brilliant trick hiding in plain sight; they are doing a more assiduous job.

tzs 6mo ago

In that Padé approximant I think you can save a couple multiplications. As written it does this: n = 1 - 367/714 * x**2 d = 1 - 81/119 * x**2 + 183/4760 * x**4 return x * (n/d) That's got 7 multiplies (I'm counting divide as a multiply) and 3 additions. (I'm assuming the optimizer only computes x^2 once and computes x^2 by squaring x^2, and that all the constants are calculated at compile time). Replace n/d with 1/(d/n) and then replace d/n with q + r/n where q is the quotient of polynomial d divided by polynomial n and r is the remainder. This is the result: n = 1 - 367/714 * x**2 q = 1587627/1346890 - 549/7340 * x**2 r = -240737/1346890 return x / (q + r/n) That's got 5 multiplies and 3 additions.