4 ms·

I wonder how does it compare with Microsoft’s implementation, there: https://github.com/microsoft/DirectXMath/blob/jan2021/Inc/DirectXMathVector.inl#L4915-L5104

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 :)
- azalemeth 5y ago
  Indeed. I love Padé approximants -- they're very underused. Another really fun trick for approximations of expensive functions are turning them into Chebsyshev functions, like Nick Trefethen's chebfun package. You can trivially differentiate, integrate and compute chebfun approximations of really horrible special functions (e.g. DE solutions) with fantastic accuracy. Criminally underused.
  - jacobolus 5y ago
    If you like Padé approximants, have you been following Trefethen & al.’s recent work on the "AAA" algorithm etc.? https://arxiv.org/pdf/1612.00337.pdf https://arxiv.org/pdf/1612.00337.pdf https://arxiv.org/pdf/1705.10132.pdf https://arxiv.org/pdf/1705.10132.pdf See also the lecture https://www.youtube.com/watch?v=S1upJPMIFfg https://www.youtube.com/watch?v=S1upJPMIFfg
    - azalemeth 5y ago
      No, I haven't -- thank you very much for sharing. I've only ever used them, e.g. for spectral decomposition in magnetic resonance in a medical setting (1) with the fast Padé transformation, which has a wonderfully deep relationship with the DEs describing magnetic resonance (2). Thanks also for the lecture - I'll enjoy. Nick Trefethen is such a good speaker. He taught me in a graduate course in numerical methods -- honestly, I'd watch those lectures on a loop again. (1) https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5179227/ https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5179227/ (2) https://iopscience.iop.org/article/10.1088/0031-9155/51/5/001/meta https://iopscience.iop.org/article/10.1088/0031-9155/51/5/00...
- petschge 5y ago
  24 bits precision would be possible with the Hastings approximation with 15 terms instead of 6 terms. A naive estimate would be that this requires 2.5 times as much work and could run at 6GB/s and take a little less than 5 cycles per element.
  - Someone 5y ago
    I think that remains to be seen. Exact calculations will give you 24 bits, but in computers, floating point calculations typically aren’t exact. If, as this article does, you compute all your intermediate values as 32-bit IEEE floats, it would hugely surprise me if you still got 24 bits of precision after 1 division (to compute y/x or x/y), 15 multiplications and 14 additions (to evaluate the polynomial) and, sometimes, a subtraction from π/2. (That’s a weakness in this article. They point out how fast their code is, but do not show any test that it actually computes approximations to atan2, let alone about how accurate their results are) Getting all the bits right in floating point is expensive.
    - adgjlsfhk1 5y ago
      The good news is that the high order terms of the polynomial don't notably affect the error. with fma, horner's method converges to .5 ULP accuracy. The subtraction with pi/2 can also be done without adding inaccuracy by using compensated summation.

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 :)
azalemeth 5y ago
Indeed. I love Padé approximants -- they're very underused. Another really fun trick for approximations of expensive functions are turning them into Chebsyshev functions, like Nick Trefethen's chebfun package. You can trivially differentiate, integrate and compute chebfun approximations of really horrible special functions (e.g. DE solutions) with fantastic accuracy. Criminally underused.
- jacobolus 5y ago
  If you like Padé approximants, have you been following Trefethen & al.’s recent work on the "AAA" algorithm etc.? https://arxiv.org/pdf/1612.00337.pdf https://arxiv.org/pdf/1612.00337.pdf https://arxiv.org/pdf/1705.10132.pdf https://arxiv.org/pdf/1705.10132.pdf See also the lecture https://www.youtube.com/watch?v=S1upJPMIFfg https://www.youtube.com/watch?v=S1upJPMIFfg
  - azalemeth 5y ago
    No, I haven't -- thank you very much for sharing. I've only ever used them, e.g. for spectral decomposition in magnetic resonance in a medical setting (1) with the fast Padé transformation, which has a wonderfully deep relationship with the DEs describing magnetic resonance (2). Thanks also for the lecture - I'll enjoy. Nick Trefethen is such a good speaker. He taught me in a graduate course in numerical methods -- honestly, I'd watch those lectures on a loop again. (1) https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5179227/ https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5179227/ (2) https://iopscience.iop.org/article/10.1088/0031-9155/51/5/001/meta https://iopscience.iop.org/article/10.1088/0031-9155/51/5/00...
petschge 5y ago
24 bits precision would be possible with the Hastings approximation with 15 terms instead of 6 terms. A naive estimate would be that this requires 2.5 times as much work and could run at 6GB/s and take a little less than 5 cycles per element.
- Someone 5y ago
  I think that remains to be seen. Exact calculations will give you 24 bits, but in computers, floating point calculations typically aren’t exact. If, as this article does, you compute all your intermediate values as 32-bit IEEE floats, it would hugely surprise me if you still got 24 bits of precision after 1 division (to compute y/x or x/y), 15 multiplications and 14 additions (to evaluate the polynomial) and, sometimes, a subtraction from π/2. (That’s a weakness in this article. They point out how fast their code is, but do not show any test that it actually computes approximations to atan2, let alone about how accurate their results are) Getting all the bits right in floating point is expensive.
  - adgjlsfhk1 5y ago
    The good news is that the high order terms of the polynomial don't notably affect the error. with fma, horner's method converges to .5 ULP accuracy. The subtraction with pi/2 can also be done without adding inaccuracy by using compensated summation.