4 ms·

This is a well done post, but please note that computers do not in fact use Taylor series for evaluating polynomials. They do use polynomials! (Or, sometimes, r

This is a well done post, but please note that computers do not in fact use Taylor series for evaluating polynomials. They do use polynomials! (Or, sometimes, rational functions, but I think the recent trend has been toward polynomials only.) They also use transformations like the log(1+x) - log(1-x) one described in this article. (That one is used in fdlibm, a very influential math library, though more modern libraries seem to not use it.) But the polynomials themselves are usually generated through some other technique, especially the Remez algorithm (https://en.wikipedia.org/wiki/Remez_algorithm https://en.wikipedia.org/wiki/Remez_algorithm) and modern improvements to it like Sollya's (https://sollya.org/ https://sollya.org/) LLL-based floating-point polynomial fitting. It's still an active area of research; the RLibM project (https://people.cs.rutgers.edu/~sn349/rlibm/ https://people.cs.rutgers.edu/~sn349/rlibm/) from Rutgers has introduced a totally new way to fit low-precision polynomials (say, up to 32 bits) using massive linear programming problems (roughly, 4 constraints on 10 variables). Source: am researcher in this field.

pclmulqdq 2y ago
Hi Pavel. For people who don't know, this is the person whose lab produced Herbie (https://herbie.uwplse.org/ https://herbie.uwplse.org/). Tacking on to this, I have seen rational approximants in "fast approximate" math libraries, while the ones targeting the highest numerical accuracy are often still using polynomials. Setting up a solver for 0.5 ULP when you have a rational function is definitely a lot harder. In my own work, it also empirically seems that the division creates some problems for precision of the last bit, so your polynomials are longer than you might expect. One interesting approach for bit-accurate rational approximants is to have a rational approximant get most of the way there while fitting the error of your rational approximant to a polynomial, but I don't think there's a good way to figure out where the computationally-efficient split of the problem is (ie how big to make the rational side and the polynomial side).
- pavpanchekha 2y ago
  Yeah, you're right on all of that. I think there have been recent breakthroughs on quadratic (I think at ARITH 23? Yep, just looked, it's the best paper: https://arith2023.arithsymposium.org/program.html https://arith2023.arithsymposium.org/program.html) so hopefully this will become more accessible with time. Though OTOH RLibM-style techniques can get you very short polynomials, to the point that it's hard to imagine beating them by much given the high cost of division operations (like 3x a multiplication).
tech_ken 2y ago
> computers do not in fact use Taylor series for evaluating polynomials. I mean I get what you're saying, but the smartass in me observes that since any polynomial is it's own Taylor expansion...
- KeplerBoy 2y ago
  I think it's a typo and supposed to read "logarithms".
  - tech_ken 2y ago
    Ahhhh that makes more sense
Sesse__ 2y ago
They used to, though! The AMD K5 paper is a classic, and it describes how they used Taylor series because they could compute them on-the-fly (they had limited ROM space). https://www.researchgate.net/publication/3612479_The_K5_transcendental_functions https://www.researchgate.net/publication/3612479_The_K5_tran...
- pclmulqdq 2y ago
  "Used to" is sort of an overstatement, IMO. The use of Taylor series here is remarkable because they need to generate coefficients on the fly. Chebyshev and Remez have been around for quite a long time, since the early 1930's, and people doing numerics in a serious fashion have generally used them for polynomial approximation since the birth of computers (except in unique circumstances like the paper you cited). The new engineering that has come recently with Sollya and similar solvers is explicit design around floating point operations. Chebyshev and Remez use the real numbers, but Sollya uses floats in its minimax function, and the coefficients you get are actually somewhat different due to rounding. Fast LP solvers have also enabled this approach more generally.
  - Sesse__ 2y ago
    It should be said that Sollya doesn't _really_ use floats. It restricts its coefficients to rationals that are (mostly) representable by floats, but the minimax is still run in full precision. Which means you can often beat it by brute force or similar.
    - pclmulqdq 2y ago
      Yes, that is a good clarification - it uses wider precision internally, but it takes floating point restrictions and operations into account. If you otherwise used Remez, you would have to just quantize the coefficients blindly and then tweak (probably manually) if something was off. Shamelessly plugging, you can sort of see the old-school process here (using integers and mixed-precision fixed point is harder with Sollya): https://specbranch.com/posts/faster-div8/ https://specbranch.com/posts/faster-div8/ The results of a follow-up quantization process that was programmatic using an LP solver to narrow the operands went into a paper I wrote for ARITH last year: https://arith2023.arithsymposium.org/papers/Newton-Raphson%20Integer%20Division%20for%20Area-Constrained%20Microcontrollers.pdf https://arith2023.arithsymposium.org/papers/Newton-Raphson%2... The coefficients and constants in both cases are substantially higher than what you get from Remez, to allow for truncation. When you have quantized operands, they are quite a bit higher. The same goes for approximations generated by Sollya - the optimal coefficients are relatively far from what Remez would tell you to do because of the error you get from rounding.
nullc 2y ago
Do you think we'll ever get to a point where there is tooling to generate implementations for user specified functions? E.g. say I want log2(x!) (Picked that example because the size of the intermediate value blocks computing it naively, and because it's a function I've had to implement before).
- thayne 2y ago
  Technically, that isn't a continuous function, so I don't think the Remez algorithm and similar would work directly. But if you replace the factorial with the gamma function, you get a continuous function that is equivalent for positive integers. And you might be able to get something from that. But there might be something more efficient that can take advantage of the fact that the range of your function is positive integers. If you just need something that doesn't overflow, you can take advantage of the properties of logarithms to compute it as the sum of log2(n) for all n from 1 to x.
  - thayne 2y ago
    And if you want an integer result, log2 can be quickly computed using the number of leading zeros.
- Cold_Miserable 2y ago
  Its not possible to automate. Log is one example. You can't just use log(x) but log(x+1). There are other problems like the "zero problem" and its still possible to devise better approximations with other elementary operations such as |x| or ABS(x).

pclmulqdq 2y ago

Hi Pavel. For people who don't know, this is the person whose lab produced Herbie (https://herbie.uwplse.org/ https://herbie.uwplse.org/). Tacking on to this, I have seen rational approximants in "fast approximate" math libraries, while the ones targeting the highest numerical accuracy are often still using polynomials. Setting up a solver for 0.5 ULP when you have a rational function is definitely a lot harder. In my own work, it also empirically seems that the division creates some problems for precision of the last bit, so your polynomials are longer than you might expect. One interesting approach for bit-accurate rational approximants is to have a rational approximant get most of the way there while fitting the error of your rational approximant to a polynomial, but I don't think there's a good way to figure out where the computationally-efficient split of the problem is (ie how big to make the rational side and the polynomial side).

pavpanchekha 2y ago
Yeah, you're right on all of that. I think there have been recent breakthroughs on quadratic (I think at ARITH 23? Yep, just looked, it's the best paper: https://arith2023.arithsymposium.org/program.html https://arith2023.arithsymposium.org/program.html) so hopefully this will become more accessible with time. Though OTOH RLibM-style techniques can get you very short polynomials, to the point that it's hard to imagine beating them by much given the high cost of division operations (like 3x a multiplication).

tech_ken 2y ago

> computers do not in fact use Taylor series for evaluating polynomials. I mean I get what you're saying, but the smartass in me observes that since any polynomial is it's own Taylor expansion...

KeplerBoy 2y ago
I think it's a typo and supposed to read "logarithms".
- tech_ken 2y ago
  Ahhhh that makes more sense

Sesse__ 2y ago

They used to, though! The AMD K5 paper is a classic, and it describes how they used Taylor series because they could compute them on-the-fly (they had limited ROM space). https://www.researchgate.net/publication/3612479_The_K5_transcendental_functions https://www.researchgate.net/publication/3612479_The_K5_tran...

pclmulqdq 2y ago
"Used to" is sort of an overstatement, IMO. The use of Taylor series here is remarkable because they need to generate coefficients on the fly. Chebyshev and Remez have been around for quite a long time, since the early 1930's, and people doing numerics in a serious fashion have generally used them for polynomial approximation since the birth of computers (except in unique circumstances like the paper you cited). The new engineering that has come recently with Sollya and similar solvers is explicit design around floating point operations. Chebyshev and Remez use the real numbers, but Sollya uses floats in its minimax function, and the coefficients you get are actually somewhat different due to rounding. Fast LP solvers have also enabled this approach more generally.
- Sesse__ 2y ago
  It should be said that Sollya doesn't _really_ use floats. It restricts its coefficients to rationals that are (mostly) representable by floats, but the minimax is still run in full precision. Which means you can often beat it by brute force or similar.
  - pclmulqdq 2y ago
    Yes, that is a good clarification - it uses wider precision internally, but it takes floating point restrictions and operations into account. If you otherwise used Remez, you would have to just quantize the coefficients blindly and then tweak (probably manually) if something was off. Shamelessly plugging, you can sort of see the old-school process here (using integers and mixed-precision fixed point is harder with Sollya): https://specbranch.com/posts/faster-div8/ https://specbranch.com/posts/faster-div8/ The results of a follow-up quantization process that was programmatic using an LP solver to narrow the operands went into a paper I wrote for ARITH last year: https://arith2023.arithsymposium.org/papers/Newton-Raphson%20Integer%20Division%20for%20Area-Constrained%20Microcontrollers.pdf https://arith2023.arithsymposium.org/papers/Newton-Raphson%2... The coefficients and constants in both cases are substantially higher than what you get from Remez, to allow for truncation. When you have quantized operands, they are quite a bit higher. The same goes for approximations generated by Sollya - the optimal coefficients are relatively far from what Remez would tell you to do because of the error you get from rounding.

nullc 2y ago

Do you think we'll ever get to a point where there is tooling to generate implementations for user specified functions? E.g. say I want log2(x!) (Picked that example because the size of the intermediate value blocks computing it naively, and because it's a function I've had to implement before).

thayne 2y ago
Technically, that isn't a continuous function, so I don't think the Remez algorithm and similar would work directly. But if you replace the factorial with the gamma function, you get a continuous function that is equivalent for positive integers. And you might be able to get something from that. But there might be something more efficient that can take advantage of the fact that the range of your function is positive integers. If you just need something that doesn't overflow, you can take advantage of the properties of logarithms to compute it as the sum of log2(n) for all n from 1 to x.
- thayne 2y ago
  And if you want an integer result, log2 can be quickly computed using the number of leading zeros.
Cold_Miserable 2y ago
Its not possible to automate. Log is one example. You can't just use log(x) but log(x+1). There are other problems like the "zero problem" and its still possible to devise better approximations with other elementary operations such as |x| or ABS(x).