12 ms·

Speeding up atan2f

245 points by rostayob 5y ago

drej 5y ago
Nice. Reminds me of an optimisation trick from a while ago: I remember being bottlenecked by one of these trigonometric functions years ago when working with a probabilistic data structure... then I figured the input domain was pretty small (a couple dozen values), so I precomputed those and used an array lookup instead. A huge win in terms of perf, obviously only applicable in these extreme cases.
- tantalor 5y ago
  https://en.wikipedia.org/wiki/Memoization https://en.wikipedia.org/wiki/Memoization
  - bluedino 5y ago
    Technically it's a lookup table if you pre-compute them. Memoization would just be caching them as you do them.
    - tantalor 5y ago
      Not necessarily, you could "cache" them in a compilation step and then use the table at runtime.
  - bruce343434 5y ago
    Tangential at best, but why was the 'r' dropped from that term? Or why not call it caching? Why the weird "memo-ization"? It makes me think of a mass extinction event where everything is turned into a memo.
    - franciscop 5y ago
      It's explained right in the linked Wikipedia page: > The term "memoization" was coined by Donald Michie in 1968[3] and is derived from the Latin word "memorandum" ("to be remembered"), usually truncated as "memo" in American English, and thus carries the meaning of "turning [the results of] a function into something to be remembered". While "memoization" might be confused with "memorization" (because they are etymological cognates), "memoization" has a specialized meaning in computing.
    - wongarsu 5y ago
      The term memoization likely precedes the word caching (as related to computing, obviously weapon caches are far older). Memoization was coined in 1968. CPU caches only came about in the 80s as registers became significantly faster than main memory. As wikipedia outlines, the r was dropped because of the memo. It's derived from the latin word memorandum that does contain the r, just like memory, but apparently it was more meant as an analogy to written memos.
- ThePadawan 5y ago
  One of the things I recently learned that sounded the most "that can't possibly work well enough" is an optimization for sin(x): If abs(x) < 0.1, "sin(x)" is approximated really well by "x". That's it. For small x, just return x. (Obviously, there is some error involved, but for the speedup gained, it's a very good compromise)
  - sharikone 5y ago
    I think that you will find that for subnormal numbers any math library will use the identity function for sin(x) and 1 for cos(x)
    - ThePadawan 5y ago
      Right, but the largest subnormal number in single-precision floats is ~ 10^-38. That the sin(x) approximation still works well for 10^-1 (with an error of ~0.01%) is the really cool thing!
  - chriswarbo 5y ago
    This is a very common assumption in Physics, e.g. https://en.wikipedia.org/wiki/Pendulum_(mathematics)#Small-angle_approximation https://en.wikipedia.org/wiki/Pendulum_(mathematics)#Small-a... Whether it's appropriate in a numerical calculation obviously depends on the possible inputs and the acceptable error bars :)
  - quietbritishjim 5y ago
    That is precisely the technique discussed in the article: it's the first term of the Taylor expansion. Except that the article used more terms of the expansion, and also used very slightly "wrong" coefficients to improve the overall accuracy within the small region.
  - Bostonian 5y ago
    Why wouldn't you at least include the x^3 term in the Taylor series for abs(x) < 0.1?
    - 5y ago
- mooman219 5y ago
  This would be a decent lookup for the atan2 function: https://gist.github.com/mooman219/19b18ff07bb9d609a103ef0cd059422c https://gist.github.com/mooman219/19b18ff07bb9d609a103ef0cd0...
drfuchs 5y ago
Wouldn’t CORDIC have done the trick faster? There’s no mention that they even considered it, even though it’s been around for half a century or so.
- adrian_b 5y ago
  CORDIC was great for devices that were too simple to include hardware multipliers. CORDIC, in its basic form, produces just 1 bit of result per clock cycle. Only for very low precisions (i.e. with few bits for each result) you would have the chance to overlap enough parallel atan2 computations to achieve a throughput comparable to what you can do with polynomial approximations on modern CPUs, with pipelined multipliers. CORDIC remains useful only for ASIC/FPGA implementations, in the cases when the area and the power consumption are much more important than the speed.
- mzs 5y ago
  CORDIC unlikely to be faster than this article's: "the approximation produces a result every 2 clock cycles"
Const-me 5y ago
I wonder how does it compare with Microsoft’s implementation, there: https://github.com/microsoft/DirectXMath/blob/jan2021/Inc/DirectXMathVector.inl#L4915-L5104 https://github.com/microsoft/DirectXMath/blob/jan2021/Inc/Di... Based on the code your version is probably much faster. It would be interesting to compare precision still, MS uses 17-degree polynomial there.
- pclmulqdq 5y ago
  The "small error" in this article isn't so small when people want exact results. Unfortunately, a high-degree polynomial is necessary if you want 24 bit precision across the input range for functions like this. That said, I wonder if a rational approximation might not be so bad on modern machines...
  - Const-me 5y ago
    > I wonder if a rational approximation might not be so bad on modern machines Uops.info says the throughput of 32-byte FP32 vector division is 5 cycles on Intel, 3.5 cycles on AMD. The throughput of 32-byte vector multiplication is 0.5 cycles. The same is for FMA, which is a single instruction computing a*b+c for 3 float vectors. If the approximation formula only has 1 division where both numerator and denominator are polynomials, and the degree of both polynomials is much lower than 17, the rational version might be not too bad compared to the 17-degree polynomial version.
    - pclmulqdq 5y ago
      A Pade approximation is generally like that. Especially when you are dealing with singularities, it should be a much smaller polynomial than 17 terms in the numerator and denominator. I think I'm going to have to start a blog, and this will be on the list.
      - janwas 5y ago
        Indeed, degree-4 rational polynomials are generally enough for JPEG XL. Code for evaluating them using Horner's scheme (just FMAs): https://github.com/libjxl/libjxl/blob/bdde644b94c125a15e532b2572b96306371a7d4e/lib/jxl/rational_polynomial-inl.h https://github.com/libjxl/libjxl/blob/bdde644b94c125a15e532b... I initially used Chebfun to approximate them, then a teammate implemented the same Caratheodory-Fejer method in C. We subsequently convert from Chebyshev basis to monomials. Blog sounds good :)
gumby 5y ago
Because of the disregard for the literature common in CS I loved this part: > This is achieved through ...and some cool documents from the 50s. A bit of anecdote: back when I was a research scientist (corporate lab) 30+ years ago I would in fact go downstairs to the library and read — I was still a kid with a lot to learn (and still am). When I came across (by chance or by someone’s suggestion) something useful to my work I’d photocopy the article and and try it out. I’d put a comment in the code with the reference. My colleagues in my group would give me (undeserved) credit for my supposed brilliance even though I said in the code where the idea came from and would determinedly point to the very paper on my desk. This attitude seemed bizarre as the group itself was producing conference papers and even books. (Obviously this was not a universal phenomenon as there were other people in the lab overall, and friends on the net, suggesting papers to read. But I’ve seen it a lot, from back then up to today)
- teddyh 5y ago
  I’ve seen an odd use of language, mostly from U.S. Americans, where “intelligent” = “smart” = “knowledgable”. And conversely, assuming that knowing many things is what makes you smart. I suspect this is what happened there; they thought you “brilliant” because you knew many things – in their minds, there might not be a difference between the two concepts.
- whatshisface 5y ago
  Not brilliant? Use of the written word to pass ideas down between generations is one of the most brilliant ideas in human history, and since nobody around you was doing it, the only explanation is that you independently invented the ethos of the scholar yourself.
  - gumby 5y ago
    Heh heh, Ockham’s razor provides an explanation!
- dleslie 5y ago
  I find that disregard for education to be tasteless, but perhaps warranted due to the modern prevalence of "bootcamps" and similar fast-but-shallow programmes. Personally, I love when I get a chance to apply something from a published paper; I will leave a citation as a comment in the code and a brief explanation of how it works and why it was chosen. I have no regrets for having achieved a computing degree, so many years ago.
prionassembly 5y ago
I wonder whether Padé approximants are well known by this kind of researcher. E.g. http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents/EfficientApproximationArctgFunction.pdf http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents...
- nimish 5y ago
  Numerical analysts are definitely familiar. Library authors may not be unless they took a course in function approximation. Personally I prefer chebyshev approximants like chebfun but pade are really good at the cost of a division.
  - stephencanon 5y ago
    A math library implementor will generally be familiar with at least minimax and Chebyshev approximations, and generally Padé and Carathéodory-Fejér approximations as well. (Source: I implemented math libraries for a decade. We used all of these frequently.)
zokier 5y ago
I would have wished to see the error analysis section expanded a bit, or maybe seeing some sort of tests to validate the max error. In particular if the mathematical approximation function arctan* has max error of 1/10000 degrees then I'd naively expect that the float-based implementation to have worse error. Furthermore it's not obvious if additional error could be introduced by the division float atan_input = (swap ? x : y) / (swap ? y : x);
sorenjan 5y ago
How do you handle arrays of values where the array lengths are not a multiple of 8 in this kind of vectorized code? Do you zero pad them before handling them to the vectorized function, or do you run a second loop element by element on the remaining elements after the main one? What happens if you try to do `_mm256_load_ps(&ys[i])` with < 8 elements remaining?
- TonyTrapp 5y ago
  Typically you have have the vectorized loop followed by a non-vectorized loop handling the remaining items, or even just explicit if-then-else statements for each possible number of remaining items.
- zbjornson 5y ago
  You could pad the input as you said, which avoids a "tail loop," but otherwise you usually do a partial load (load <8 elements into a vector) and store. Some instruction set extensions provide "masked load/store" instructions for this, but there are ways to do it without those too. To your last question specifically, if you _mm256_load_ps(&ys[i]) and you're at the edge of a page, you'll get a segfault. Otherwise you'll get undefined values (which can be okay, you could ignore them instead of dealing with a partial load).
- vlovich123 5y ago
  > Do you zero pad them before handling them to the vectorized function, or do you run a second loop element by element on the remaining elements after the main one? Typically the latter. That's why I find SVEs so interesting. Should improve code density.
- twic 5y ago
  RISC-V has a solution to this which seems quite elegant - the vector length is set by the program, up to some maximum, and for the final block, it is simply whatever is left: https://gms.tf/riscv-vector.html https://gms.tf/riscv-vector.html
deleted 5y ago
[deleted]
h0mie 5y ago
Love posts like this!
aconz2 5y ago
Nice writeup and interesting results. I hadn't seen the use of perf_event_open(2) before directly in code which looks cool. The baseline is at a huge disadvantage here because the call to atan2 in the loop never gets inlined and the loop doesn't seem to get unrolled (which is surprising actually). Manually unrolling by 8 gives me an 8x speedup. Maybe I'm missing something with the `-static` link but unless they're using musl I didn't think -lm could get statically linked.
- xxpor 5y ago
  Could it be that calls to glibc never get inlined, since inlining is essentially static linking by another name? Or since they're in separate compilation units, without LTO you'd never perform the analysis to figure out if it's worth it. Would LTO inspect across a dynamic loading boundary? Just speculating, I really have no idea. Everything I work with is statically linked, so I've never really had to think about it.
- aconz2 5y ago
  You've made a silly mistake and just did 1/8th the work, so of course it was 8x speedup. I am still wondering how the numbers would look if you could inline atan2f from libc. Too bad that isn't easier
jvz01 5y ago
I have developed very fast, accurate, and vectorizable atan() and atan2() implementations, leveraging AVX/SSE capabilities. You can find them here [warning: self-signed SSL-Cert]. https://fox-toolkit.org/wordpress/?p=219 https://fox-toolkit.org/wordpress/?p=219
- ghusbands 5y ago
  Your result is significantly slower than the versions presented in the article, though yours has more terms and so may be more accurate.
  - jvz01 5y ago
    Yes, aim was to be acurate down to 1 lsb while significantly faster. Feel free to drop terms from the polynomial if you can live with less accurate results! The coefficients were generated by a package called Sollya, I've used it a few times to develop accurate chebyshev approximations for functions. Abramowitz & Stegun is another good reference.
    - touisteur 5y ago
      Please, Would you mind one of these days updating your blog post with the instructions you gave to sollya? I'm trying something stupid with log1p and can't get sollya to help, mostly because I'm not putting enough time to read all the docs...
- jvz01 5y ago
  Little side-note: algorithm as given is scalar; however, its branch-free, and defined entirely in the header file. So, compilers will typically be able to vectorize it, and thus achieve speed up directly based on the vector size. I see potential [but architecture-dependent] optimization using Estrin scheme for evaluating the polynomial.
jacobolus 5y ago
If someone wants a fast version of x ↦ tan(πx/2), let me recommend the approximation: tanpi_2 = function tanpi_2(x) { var y = (1 - x*x); return x * (((-0.000221184 * y + 0.0024971104) * y - 0.02301937096) * y + 0.3182994604 + 1.2732402998 / y); } (valid for -1 <= x <= 1) https://observablehq.com/@jrus/fasttan https://observablehq.com/@jrus/fasttan with error: https://www.desmos.com/calculator/hmncdd6fuj https://www.desmos.com/calculator/hmncdd6fuj But even better is to avoid trigonometry and angle measures as much as possible. Almost everything can be done better (faster, with fewer numerical problems) with vector methods; if you want a 1-float representation of an angle, use the stereographic projection: stereo = (x, y) => y/(x + Math.hypot(x, y)); stereo_to_xy = (s) => { var q = 1/(1 + s*s); return !q ? [-1, 0] : [(1 - s*s)/q, 2*s/q]; }
- leowbattle 5y ago
  > Almost everything can be done better ... use the stereographic projection I was just learning about stereographic projection earlier. Isn't it odd how when you know something you notice it appearing in places. Can you give an example of an operation that could be performed better using stereographic projection rather than angles?
  - jacobolus 5y ago
    Generally you can just stick to storing 2-coordinate vectors and using vector operations. The places where you might want to convert to a 1-number representation are when you have a lot of numbers you want to store or transmit somewhere. Using the stereographic projection (half-angle tangent) instead of angle measure works better with the floating point number format and uses only rational arithmetic instead of transcendental functions.
hdersch 5y ago
A comparison with one of the many SIMD-mathlibraries would have been fairer than with plain libm. Long time ago I wrote such a dual-platform library for the PS3 (cell-processor) and x86 architecture (outdated, but still available here [1]). Depending on how standard libm implements atan2f, a speedup of 3x to 15x is achieved, without sacrifying accuracy. 1. https://webuser.hs-furtwangen.de/~dersch/libsimdmath.pdf https://webuser.hs-furtwangen.de/~dersch/libsimdmath.pdf
aj7 5y ago
Around 1988, I added phase shift to my optical thin film design program written in Excel 4.0 for the Mac. At the time, this was utterly unique: each spreadsheet row represented a layer and the matrices describing each layer could be calculated right in that row by squashing them down horizontally. The S- and P-polarization matrices could be recorded this way, and the running matrix products similarly maintained. Finally, using a simple one input table, the reflectance of a typically 25-31 layer laser mirror could be calculated. And in less than a second on a 20 MHz 68020 (?) Mac II for about 50 wavelengths. The best part were the graphics which were instantaneous, beautiful, publishable, and pasteable into customer quotations. Semi-technical people could be trained to use the whole thing. Now about the phase shift. In 1988, atan2 didn’t exist. Anywhere. Not in FORTRAN, Basic, Excel, or a C library. I’m sure phase shift calculators implemented it, each working alone. For us, it was critical. You see we actually cared not about the phase shift, but its second derivative, the group delay dispersion. This was the beginning of the femtosecond laser era, and people needed to check whether these broadband laser pulses would be inadvertently stretched by reflection off or transmission through the mirror coating. So atan2, the QUADRANT-PRESERVING arc tangent, is required for a continuous,differential phase function. An Excel function macro did this, with IF statements correcting the quadrant. And the irony of all this? I CALLED it atan2.
- floxy 5y ago
  Table 8-3, page 8-7 (pdf page 97): http://www.bitsavers.org/www.computer.museum.uq.edu.au/pdf/DEC-10-AFDO-D%20decsystem10%20FORTRAN%20IV%20Programmer's%20Reference%20Manual.pdf http://www.bitsavers.org/www.computer.museum.uq.edu.au/pdf/D...
  - raphlinus 5y ago
    1965: https://www.ecma-international.org/publications-and-standards/standards/ecma-9/ https://www.ecma-international.org/publications-and-standard... pdf page 36
    - floxy 5y ago
      Nice find. Anyone have a way to search for old versions "math.h"? The 1989 version of ANSI C had atan2: https://web.archive.org/web/20161223125339/http://flash-gordon.me.uk/ansi.c.txt https://web.archive.org/web/20161223125339/http://flash-gord... ...but I wonder when it first hit the scene.
nice2meetu 5y ago
I did something similar for tanh once, though I found I could get to 1 ulp. Part of the motivation was that I could get 10x faster than libc. However, I then tried on my FreeBSD and could only get 4x faster. After a lot of head scratching and puzzling it turned out there was a bug in the version of libc on my linux box that slowed things down. It kind of took the wind out of the achievement, but it was still a great learning experience.
spaetzleesser 5y ago
I am envious of people who can deal with such problems. The problem is defined clearly and can be measured easily. This is so much more fun than figuring out why some SAAS service is misbehaving.
unemphysbro 5y ago
coolest blog post I've seen here in a while. :)
cogman10 5y ago
I'm actually a bit surprised that the x86 SIMD instructions don't support trig functions.
azhenley 5y ago
This is pretty similar to my quest to make my own cos() when my friend didn't have access to libc. It was fun! Though I don't have the math or low-level knowledge that this author does. https://web.eecs.utk.edu/~azh/blog/cosine.html https://web.eecs.utk.edu/~azh/blog/cosine.html
pklausler 5y ago
(Undoubtedly) stupid question: would it be any faster to project (x, y) to the unit circle (x', y'), then compute acos(x') or asin(y'), and then correct the result based on the signs of x & y? When converting Cartesian coordinates to polar, the value of r=HYPOT(x, y) is needed anyway, so the projection to the unit circle would be a single division by r.
stephencanon 5y ago
> if we’re working with batches of points and willing to live with tiny errors, we can produce an atan2 approximation which is 50 times faster than the standard version provided by libc. Which libc, though? I assume glibc, but it's frustrating when people talk about libc as though there were a single implementation. Each vendor supplies their own implementation, libc is just a common interface defined by the C library. There is no "standard version" provided by libc. In particular, glibc's math functions are not especially fast--Intel's and Apple's math libraries are 4-5x faster for some functions[1], and often more accurate as well, for example (and both vendors provide vectorized implementations). Even within glibc versions, there have been enormous improvements over the last decade or so, and for some functions there are big performance differences depending on whether or not -fno-math-errno is specified. (I would also note that atan2 has a lot of edge cases, and more than half the work in a standards-compliant libc is in getting those edge cases with zeros and infinities right, which this implementation punts on. There's nothing wrong with that, but that's a bigger tradeoff for most users than the small loss of accuracy and important to note.) So what are we actually comparing against here? Comparing against a clown-shoes baseline makes for eye-popping numbers, but it's not very meaningful. None of this should really take away from the work presented, by the way--the techniques described here are very useful for people interested in this stuff. [1] I don't know the current state of atan2f in glibc specifically; it's possible that it's been improved since last I looked at its performance. But the blog post cites "105.98 cycles / element", which would be glacially slow on any semi-recent hardware, which makes me think something is up here.
- rostayob 5y ago
  (I'm the author) You're right, I should have specified -- it is glibc 2.32-48 . This the source specifying how glibc is built: https://github.com/NixOS/nixpkgs/blob/97c5d0cbe76901da0135b05cdbdfc5b068a7942c/pkgs/development/libraries/glibc/default.nix https://github.com/NixOS/nixpkgs/blob/97c5d0cbe76901da0135b0... . I've amended the article so that it says `glibc` rather than `libc`, and added a sidenote specifying the version. I link to it statically as indicated in the gist, although I believe that shouldn't matter. Also see https://gist.github.com/bitonic/d0f5a0a44e37d4f0be03d34d47acb6cf#gistcomment-3864005 https://gist.github.com/bitonic/d0f5a0a44e37d4f0be03d34d47ac... . Note that the hardware is not particularly recent (Q3 2017), but we tend to rent servers which are not exactly on the bleeding edge, so that was my platform.
shoo 5y ago
Related -- there's a 2011 post from Paul Minero with fast approximations for logarithm, exponential, power, inverse root. http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html http://www.machinedlearnings.com/2011/06/fast-approximate-lo... Minero's faster approximate log2, < 1.4% relative error for x in [1/100, 10]. Here's the simple non-sse version: static inline float fasterlog2 (float x) { union { float f; uint32_t i; } vx = { x }; float y = vx.i; y *= 1.1920928955078125e-7f; return y - 126.94269504f; } This fastapprox library also includes fast approximations of some other functions that show up in statistical / probabilistic calculations -- gamma, digamma, lambert w function. It is BSD licensed, originally lived in google code, copies of the library live on in github, e.g. https://github.com/etheory/fastapprox https://github.com/etheory/fastapprox It's also interesting to read through libm. E.g. compare Sun's ~1993 atan2 & atan: https://github.com/JuliaMath/openlibm/blob/master/src/e_atan2.c https://github.com/JuliaMath/openlibm/blob/master/src/e_atan... https://github.com/JuliaMath/openlibm/blob/master/src/s_atan.c https://github.com/JuliaMath/openlibm/blob/master/src/s_atan...

245 points by rostayob 5y ago

drej 5y ago

Nice. Reminds me of an optimisation trick from a while ago: I remember being bottlenecked by one of these trigonometric functions years ago when working with a probabilistic data structure... then I figured the input domain was pretty small (a couple dozen values), so I precomputed those and used an array lookup instead. A huge win in terms of perf, obviously only applicable in these extreme cases.

tantalor 5y ago
https://en.wikipedia.org/wiki/Memoization https://en.wikipedia.org/wiki/Memoization
- bluedino 5y ago
  Technically it's a lookup table if you pre-compute them. Memoization would just be caching them as you do them.
  - tantalor 5y ago
    Not necessarily, you could "cache" them in a compilation step and then use the table at runtime.
- bruce343434 5y ago
  Tangential at best, but why was the 'r' dropped from that term? Or why not call it caching? Why the weird "memo-ization"? It makes me think of a mass extinction event where everything is turned into a memo.
  - franciscop 5y ago
    It's explained right in the linked Wikipedia page: > The term "memoization" was coined by Donald Michie in 1968[3] and is derived from the Latin word "memorandum" ("to be remembered"), usually truncated as "memo" in American English, and thus carries the meaning of "turning [the results of] a function into something to be remembered". While "memoization" might be confused with "memorization" (because they are etymological cognates), "memoization" has a specialized meaning in computing.
  - wongarsu 5y ago
    The term memoization likely precedes the word caching (as related to computing, obviously weapon caches are far older). Memoization was coined in 1968. CPU caches only came about in the 80s as registers became significantly faster than main memory. As wikipedia outlines, the r was dropped because of the memo. It's derived from the latin word memorandum that does contain the r, just like memory, but apparently it was more meant as an analogy to written memos.
ThePadawan 5y ago
One of the things I recently learned that sounded the most "that can't possibly work well enough" is an optimization for sin(x): If abs(x) < 0.1, "sin(x)" is approximated really well by "x". That's it. For small x, just return x. (Obviously, there is some error involved, but for the speedup gained, it's a very good compromise)
- sharikone 5y ago
  I think that you will find that for subnormal numbers any math library will use the identity function for sin(x) and 1 for cos(x)
  - ThePadawan 5y ago
    Right, but the largest subnormal number in single-precision floats is ~ 10^-38. That the sin(x) approximation still works well for 10^-1 (with an error of ~0.01%) is the really cool thing!
- chriswarbo 5y ago
  This is a very common assumption in Physics, e.g. https://en.wikipedia.org/wiki/Pendulum_(mathematics)#Small-angle_approximation https://en.wikipedia.org/wiki/Pendulum_(mathematics)#Small-a... Whether it's appropriate in a numerical calculation obviously depends on the possible inputs and the acceptable error bars :)
- quietbritishjim 5y ago
  That is precisely the technique discussed in the article: it's the first term of the Taylor expansion. Except that the article used more terms of the expansion, and also used very slightly "wrong" coefficients to improve the overall accuracy within the small region.
- Bostonian 5y ago
  Why wouldn't you at least include the x^3 term in the Taylor series for abs(x) < 0.1?
  - 5y ago
mooman219 5y ago
This would be a decent lookup for the atan2 function: https://gist.github.com/mooman219/19b18ff07bb9d609a103ef0cd059422c https://gist.github.com/mooman219/19b18ff07bb9d609a103ef0cd0...

drfuchs 5y ago

Wouldn’t CORDIC have done the trick faster? There’s no mention that they even considered it, even though it’s been around for half a century or so.

adrian_b 5y ago
CORDIC was great for devices that were too simple to include hardware multipliers. CORDIC, in its basic form, produces just 1 bit of result per clock cycle. Only for very low precisions (i.e. with few bits for each result) you would have the chance to overlap enough parallel atan2 computations to achieve a throughput comparable to what you can do with polynomial approximations on modern CPUs, with pipelined multipliers. CORDIC remains useful only for ASIC/FPGA implementations, in the cases when the area and the power consumption are much more important than the speed.
mzs 5y ago
CORDIC unlikely to be faster than this article's: "the approximation produces a result every 2 clock cycles"

Const-me 5y ago

I wonder how does it compare with Microsoft’s implementation, there: https://github.com/microsoft/DirectXMath/blob/jan2021/Inc/DirectXMathVector.inl#L4915-L5104 https://github.com/microsoft/DirectXMath/blob/jan2021/Inc/Di... Based on the code your version is probably much faster. It would be interesting to compare precision still, MS uses 17-degree polynomial there.

pclmulqdq 5y ago
The "small error" in this article isn't so small when people want exact results. Unfortunately, a high-degree polynomial is necessary if you want 24 bit precision across the input range for functions like this. That said, I wonder if a rational approximation might not be so bad on modern machines...
- Const-me 5y ago
  > I wonder if a rational approximation might not be so bad on modern machines Uops.info says the throughput of 32-byte FP32 vector division is 5 cycles on Intel, 3.5 cycles on AMD. The throughput of 32-byte vector multiplication is 0.5 cycles. The same is for FMA, which is a single instruction computing a*b+c for 3 float vectors. If the approximation formula only has 1 division where both numerator and denominator are polynomials, and the degree of both polynomials is much lower than 17, the rational version might be not too bad compared to the 17-degree polynomial version.
  - pclmulqdq 5y ago
    A Pade approximation is generally like that. Especially when you are dealing with singularities, it should be a much smaller polynomial than 17 terms in the numerator and denominator. I think I'm going to have to start a blog, and this will be on the list.
    - janwas 5y ago
      Indeed, degree-4 rational polynomials are generally enough for JPEG XL. Code for evaluating them using Horner's scheme (just FMAs): https://github.com/libjxl/libjxl/blob/bdde644b94c125a15e532b2572b96306371a7d4e/lib/jxl/rational_polynomial-inl.h https://github.com/libjxl/libjxl/blob/bdde644b94c125a15e532b... I initially used Chebfun to approximate them, then a teammate implemented the same Caratheodory-Fejer method in C. We subsequently convert from Chebyshev basis to monomials. Blog sounds good :)

gumby 5y ago

Because of the disregard for the literature common in CS I loved this part: > This is achieved through ...and some cool documents from the 50s. A bit of anecdote: back when I was a research scientist (corporate lab) 30+ years ago I would in fact go downstairs to the library and read — I was still a kid with a lot to learn (and still am). When I came across (by chance or by someone’s suggestion) something useful to my work I’d photocopy the article and and try it out. I’d put a comment in the code with the reference. My colleagues in my group would give me (undeserved) credit for my supposed brilliance even though I said in the code where the idea came from and would determinedly point to the very paper on my desk. This attitude seemed bizarre as the group itself was producing conference papers and even books. (Obviously this was not a universal phenomenon as there were other people in the lab overall, and friends on the net, suggesting papers to read. But I’ve seen it a lot, from back then up to today)

teddyh 5y ago
I’ve seen an odd use of language, mostly from U.S. Americans, where “intelligent” = “smart” = “knowledgable”. And conversely, assuming that knowing many things is what makes you smart. I suspect this is what happened there; they thought you “brilliant” because you knew many things – in their minds, there might not be a difference between the two concepts.
whatshisface 5y ago
Not brilliant? Use of the written word to pass ideas down between generations is one of the most brilliant ideas in human history, and since nobody around you was doing it, the only explanation is that you independently invented the ethos of the scholar yourself.
- gumby 5y ago
  Heh heh, Ockham’s razor provides an explanation!
dleslie 5y ago
I find that disregard for education to be tasteless, but perhaps warranted due to the modern prevalence of "bootcamps" and similar fast-but-shallow programmes. Personally, I love when I get a chance to apply something from a published paper; I will leave a citation as a comment in the code and a brief explanation of how it works and why it was chosen. I have no regrets for having achieved a computing degree, so many years ago.

prionassembly 5y ago

I wonder whether Padé approximants are well known by this kind of researcher. E.g. http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents/EfficientApproximationArctgFunction.pdf http://www-labs.iro.umontreal.ca/~mignotte/IFT2425/Documents...

nimish 5y ago
Numerical analysts are definitely familiar. Library authors may not be unless they took a course in function approximation. Personally I prefer chebyshev approximants like chebfun but pade are really good at the cost of a division.
- stephencanon 5y ago
  A math library implementor will generally be familiar with at least minimax and Chebyshev approximations, and generally Padé and Carathéodory-Fejér approximations as well. (Source: I implemented math libraries for a decade. We used all of these frequently.)

zokier 5y ago

I would have wished to see the error analysis section expanded a bit, or maybe seeing some sort of tests to validate the max error. In particular if the mathematical approximation function arctan* has max error of 1/10000 degrees then I'd naively expect that the float-based implementation to have worse error. Furthermore it's not obvious if additional error could be introduced by the division float atan_input = (swap ? x : y) / (swap ? y : x);

sorenjan 5y ago

How do you handle arrays of values where the array lengths are not a multiple of 8 in this kind of vectorized code? Do you zero pad them before handling them to the vectorized function, or do you run a second loop element by element on the remaining elements after the main one? What happens if you try to do `_mm256_load_ps(&ys[i])` with < 8 elements remaining?

TonyTrapp 5y ago
Typically you have have the vectorized loop followed by a non-vectorized loop handling the remaining items, or even just explicit if-then-else statements for each possible number of remaining items.
zbjornson 5y ago
You could pad the input as you said, which avoids a "tail loop," but otherwise you usually do a partial load (load <8 elements into a vector) and store. Some instruction set extensions provide "masked load/store" instructions for this, but there are ways to do it without those too. To your last question specifically, if you _mm256_load_ps(&ys[i]) and you're at the edge of a page, you'll get a segfault. Otherwise you'll get undefined values (which can be okay, you could ignore them instead of dealing with a partial load).
vlovich123 5y ago
> Do you zero pad them before handling them to the vectorized function, or do you run a second loop element by element on the remaining elements after the main one? Typically the latter. That's why I find SVEs so interesting. Should improve code density.
twic 5y ago
RISC-V has a solution to this which seems quite elegant - the vector length is set by the program, up to some maximum, and for the final block, it is simply whatever is left: https://gms.tf/riscv-vector.html https://gms.tf/riscv-vector.html

deleted 5y ago

[deleted]

h0mie 5y ago

Love posts like this!

aconz2 5y ago

Nice writeup and interesting results. I hadn't seen the use of perf_event_open(2) before directly in code which looks cool. The baseline is at a huge disadvantage here because the call to atan2 in the loop never gets inlined and the loop doesn't seem to get unrolled (which is surprising actually). Manually unrolling by 8 gives me an 8x speedup. Maybe I'm missing something with the `-static` link but unless they're using musl I didn't think -lm could get statically linked.

xxpor 5y ago
Could it be that calls to glibc never get inlined, since inlining is essentially static linking by another name? Or since they're in separate compilation units, without LTO you'd never perform the analysis to figure out if it's worth it. Would LTO inspect across a dynamic loading boundary? Just speculating, I really have no idea. Everything I work with is statically linked, so I've never really had to think about it.
aconz2 5y ago
You've made a silly mistake and just did 1/8th the work, so of course it was 8x speedup. I am still wondering how the numbers would look if you could inline atan2f from libc. Too bad that isn't easier

jvz01 5y ago

I have developed very fast, accurate, and vectorizable atan() and atan2() implementations, leveraging AVX/SSE capabilities. You can find them here [warning: self-signed SSL-Cert]. https://fox-toolkit.org/wordpress/?p=219 https://fox-toolkit.org/wordpress/?p=219

ghusbands 5y ago
Your result is significantly slower than the versions presented in the article, though yours has more terms and so may be more accurate.
- jvz01 5y ago
  Yes, aim was to be acurate down to 1 lsb while significantly faster. Feel free to drop terms from the polynomial if you can live with less accurate results! The coefficients were generated by a package called Sollya, I've used it a few times to develop accurate chebyshev approximations for functions. Abramowitz & Stegun is another good reference.
  - touisteur 5y ago
    Please, Would you mind one of these days updating your blog post with the instructions you gave to sollya? I'm trying something stupid with log1p and can't get sollya to help, mostly because I'm not putting enough time to read all the docs...
jvz01 5y ago
Little side-note: algorithm as given is scalar; however, its branch-free, and defined entirely in the header file. So, compilers will typically be able to vectorize it, and thus achieve speed up directly based on the vector size. I see potential [but architecture-dependent] optimization using Estrin scheme for evaluating the polynomial.

jacobolus 5y ago

If someone wants a fast version of x ↦ tan(πx/2), let me recommend the approximation: tanpi_2 = function tanpi_2(x) { var y = (1 - x*x); return x * (((-0.000221184 * y + 0.0024971104) * y - 0.02301937096) * y + 0.3182994604 + 1.2732402998 / y); } (valid for -1 <= x <= 1) https://observablehq.com/@jrus/fasttan https://observablehq.com/@jrus/fasttan with error: https://www.desmos.com/calculator/hmncdd6fuj https://www.desmos.com/calculator/hmncdd6fuj But even better is to avoid trigonometry and angle measures as much as possible. Almost everything can be done better (faster, with fewer numerical problems) with vector methods; if you want a 1-float representation of an angle, use the stereographic projection: stereo = (x, y) => y/(x + Math.hypot(x, y)); stereo_to_xy = (s) => { var q = 1/(1 + s*s); return !q ? [-1, 0] : [(1 - s*s)/q, 2*s/q]; }

leowbattle 5y ago
> Almost everything can be done better ... use the stereographic projection I was just learning about stereographic projection earlier. Isn't it odd how when you know something you notice it appearing in places. Can you give an example of an operation that could be performed better using stereographic projection rather than angles?
- jacobolus 5y ago
  Generally you can just stick to storing 2-coordinate vectors and using vector operations. The places where you might want to convert to a 1-number representation are when you have a lot of numbers you want to store or transmit somewhere. Using the stereographic projection (half-angle tangent) instead of angle measure works better with the floating point number format and uses only rational arithmetic instead of transcendental functions.

hdersch 5y ago

A comparison with one of the many SIMD-mathlibraries would have been fairer than with plain libm. Long time ago I wrote such a dual-platform library for the PS3 (cell-processor) and x86 architecture (outdated, but still available here [1]). Depending on how standard libm implements atan2f, a speedup of 3x to 15x is achieved, without sacrifying accuracy. 1. https://webuser.hs-furtwangen.de/~dersch/libsimdmath.pdf https://webuser.hs-furtwangen.de/~dersch/libsimdmath.pdf

aj7 5y ago

Around 1988, I added phase shift to my optical thin film design program written in Excel 4.0 for the Mac. At the time, this was utterly unique: each spreadsheet row represented a layer and the matrices describing each layer could be calculated right in that row by squashing them down horizontally. The S- and P-polarization matrices could be recorded this way, and the running matrix products similarly maintained. Finally, using a simple one input table, the reflectance of a typically 25-31 layer laser mirror could be calculated. And in less than a second on a 20 MHz 68020 (?) Mac II for about 50 wavelengths. The best part were the graphics which were instantaneous, beautiful, publishable, and pasteable into customer quotations. Semi-technical people could be trained to use the whole thing. Now about the phase shift. In 1988, atan2 didn’t exist. Anywhere. Not in FORTRAN, Basic, Excel, or a C library. I’m sure phase shift calculators implemented it, each working alone. For us, it was critical. You see we actually cared not about the phase shift, but its second derivative, the group delay dispersion. This was the beginning of the femtosecond laser era, and people needed to check whether these broadband laser pulses would be inadvertently stretched by reflection off or transmission through the mirror coating. So atan2, the QUADRANT-PRESERVING arc tangent, is required for a continuous,differential phase function. An Excel function macro did this, with IF statements correcting the quadrant. And the irony of all this? I CALLED it atan2.

floxy 5y ago
Table 8-3, page 8-7 (pdf page 97): http://www.bitsavers.org/www.computer.museum.uq.edu.au/pdf/DEC-10-AFDO-D%20decsystem10%20FORTRAN%20IV%20Programmer's%20Reference%20Manual.pdf http://www.bitsavers.org/www.computer.museum.uq.edu.au/pdf/D...
- raphlinus 5y ago
  1965: https://www.ecma-international.org/publications-and-standards/standards/ecma-9/ https://www.ecma-international.org/publications-and-standard... pdf page 36
  - floxy 5y ago
    Nice find. Anyone have a way to search for old versions "math.h"? The 1989 version of ANSI C had atan2: https://web.archive.org/web/20161223125339/http://flash-gordon.me.uk/ansi.c.txt https://web.archive.org/web/20161223125339/http://flash-gord... ...but I wonder when it first hit the scene.

nice2meetu 5y ago

I did something similar for tanh once, though I found I could get to 1 ulp. Part of the motivation was that I could get 10x faster than libc. However, I then tried on my FreeBSD and could only get 4x faster. After a lot of head scratching and puzzling it turned out there was a bug in the version of libc on my linux box that slowed things down. It kind of took the wind out of the achievement, but it was still a great learning experience.

spaetzleesser 5y ago

I am envious of people who can deal with such problems. The problem is defined clearly and can be measured easily. This is so much more fun than figuring out why some SAAS service is misbehaving.

unemphysbro 5y ago

coolest blog post I've seen here in a while. :)

cogman10 5y ago

I'm actually a bit surprised that the x86 SIMD instructions don't support trig functions.

azhenley 5y ago

This is pretty similar to my quest to make my own cos() when my friend didn't have access to libc. It was fun! Though I don't have the math or low-level knowledge that this author does. https://web.eecs.utk.edu/~azh/blog/cosine.html https://web.eecs.utk.edu/~azh/blog/cosine.html

pklausler 5y ago

(Undoubtedly) stupid question: would it be any faster to project (x, y) to the unit circle (x', y'), then compute acos(x') or asin(y'), and then correct the result based on the signs of x & y? When converting Cartesian coordinates to polar, the value of r=HYPOT(x, y) is needed anyway, so the projection to the unit circle would be a single division by r.

stephencanon 5y ago

> if we’re working with batches of points and willing to live with tiny errors, we can produce an atan2 approximation which is 50 times faster than the standard version provided by libc. Which libc, though? I assume glibc, but it's frustrating when people talk about libc as though there were a single implementation. Each vendor supplies their own implementation, libc is just a common interface defined by the C library. There is no "standard version" provided by libc. In particular, glibc's math functions are not especially fast--Intel's and Apple's math libraries are 4-5x faster for some functions[1], and often more accurate as well, for example (and both vendors provide vectorized implementations). Even within glibc versions, there have been enormous improvements over the last decade or so, and for some functions there are big performance differences depending on whether or not -fno-math-errno is specified. (I would also note that atan2 has a lot of edge cases, and more than half the work in a standards-compliant libc is in getting those edge cases with zeros and infinities right, which this implementation punts on. There's nothing wrong with that, but that's a bigger tradeoff for most users than the small loss of accuracy and important to note.) So what are we actually comparing against here? Comparing against a clown-shoes baseline makes for eye-popping numbers, but it's not very meaningful. None of this should really take away from the work presented, by the way--the techniques described here are very useful for people interested in this stuff. [1] I don't know the current state of atan2f in glibc specifically; it's possible that it's been improved since last I looked at its performance. But the blog post cites "105.98 cycles / element", which would be glacially slow on any semi-recent hardware, which makes me think something is up here.

rostayob 5y ago
(I'm the author) You're right, I should have specified -- it is glibc 2.32-48 . This the source specifying how glibc is built: https://github.com/NixOS/nixpkgs/blob/97c5d0cbe76901da0135b05cdbdfc5b068a7942c/pkgs/development/libraries/glibc/default.nix https://github.com/NixOS/nixpkgs/blob/97c5d0cbe76901da0135b0... . I've amended the article so that it says `glibc` rather than `libc`, and added a sidenote specifying the version. I link to it statically as indicated in the gist, although I believe that shouldn't matter. Also see https://gist.github.com/bitonic/d0f5a0a44e37d4f0be03d34d47acb6cf#gistcomment-3864005 https://gist.github.com/bitonic/d0f5a0a44e37d4f0be03d34d47ac... . Note that the hardware is not particularly recent (Q3 2017), but we tend to rent servers which are not exactly on the bleeding edge, so that was my platform.

shoo 5y ago

Related -- there's a 2011 post from Paul Minero with fast approximations for logarithm, exponential, power, inverse root. http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html http://www.machinedlearnings.com/2011/06/fast-approximate-lo... Minero's faster approximate log2, < 1.4% relative error for x in [1/100, 10]. Here's the simple non-sse version: static inline float fasterlog2 (float x) { union { float f; uint32_t i; } vx = { x }; float y = vx.i; y *= 1.1920928955078125e-7f; return y - 126.94269504f; } This fastapprox library also includes fast approximations of some other functions that show up in statistical / probabilistic calculations -- gamma, digamma, lambert w function. It is BSD licensed, originally lived in google code, copies of the library live on in github, e.g. https://github.com/etheory/fastapprox https://github.com/etheory/fastapprox It's also interesting to read through libm. E.g. compare Sun's ~1993 atan2 & atan: https://github.com/JuliaMath/openlibm/blob/master/src/e_atan2.c https://github.com/JuliaMath/openlibm/blob/master/src/e_atan... https://github.com/JuliaMath/openlibm/blob/master/src/s_atan.c https://github.com/JuliaMath/openlibm/blob/master/src/s_atan...