2. Computers actually use this recursion to calculate sines and cosines! So, it's a little inefficient to code your Chebyshev polynomials using `math.sin`.
3. Using generating functions, you can get a closed form for T_n(x) that only takes O(log n) time to calculate. (Note: assuming you count multiplications as constant. However, you actually need O(log n) bits to accurately represent x, so it's more accurately O((log n)^2 log log n).)
If you are ready to spend some precomputation time to compute a good approximation, you can use the Remez algorithm [1]. It is implemented in the Sollya library for machine precision [2,3]. It has notably been used to implement the Core Math library [4] to provide correct rounding for the math functions in the libc library.
This is so strange: a few days ago I commented on an HN post (https://news.ycombinator.com/item?id=40582712) about when "programmer" became an acknowledged job title, and, in my comment, mentioned how I used a Chebyshev approximation followed by two iterations of Newton's method to compute sqrt, and then today this article shows up with exactly that use case!
I wrote that code (to compute sqrt 2) in 1974 or 1975 in IBM 360 Assembler. I used a conditional macro constant that increased the number of iterations of Newton's from 2 to 3 just in case the client wanted double precision.
chebyshev approximations are fucking awesome, but this article gives too short shrift to table lookup; it does go a bit beyond nearest-neighbor interpolation to linear interpolation, and correctly points out that this gives you error that is quadratic in the distance from the x-coordinate of the nearest table entry (and therefore worst-case error quadratic in your point spacing), and that this gives you half the error of the fifth-order chebyshev approximation. it says that this is a 'rare case', but in fact you will always get a lower error from table lookup if you use enough points. it's just that with only linear interpolation, the number of points rapidly becomes impractical
as i understand it, other commonly-used strategies include spline interpolation (using second-, third-, or even fourth-order interpolation, requiring respectively three, four, and five multiplications, which can be done concurrently) and, in suitable cases like this example, newton iteration from an initial table-lookup guess
unlike the piecewise-taylor approach outlined early in the article, spline interpolation only requires storing a tiny amount more data than simple nearest-neighbor table lookup (potentially three more points for fourth-order interpolation, so a 256-entry table becomes 259 entries)
on a different topic, i think it's easy to find embedded dsp applications where the easiest solution uses fourier transforms, which usually do require high-precision floating point. machine vision, radio communication, musical applications, etc.
incidentally, if you find yourself in a situation where you actually need the taylor expansion of √(1+x) or √(½+x) or something, and you don't want to do a bunch of pencil-and-paper algebra (or don't trust yourself), pari/gp has your back:
How does that feature work without a solution to the tablemaker's dilemma? Does the compiler just give up and give you arbitrary precision every time you use transcendentals or does it give up and use lower precision if the estimate exceeds some arbitrary bound?
Looking at the spec [0], it only demands 1 ULP of precision for the elementary operations, and 2, 4, or 8 ULPs for special functions. So there's no magic 1/2-ULP guarantee.
aaui that Ada line would declare a fixed-point type, so decimal/binary shouldn't really come up except for overflow handling, here. (Haven't tried it though.)
1. You can use a Fourier transform to get the coefficients in O(n log n) time.
2. So, multiplying two approximations only takes O(n log n) time.
3. Also, adding, integrating, or taking the derivative only take O(n) time.
This is why chebfun/chebpy can run so fast while magically finding roots/derivatives/etc. A couple other interesting facts:
1. Remember the double-angle formula? There's a more general recursion for the Chebyshev polynomials:
\[ T_n(x) = 2x T_{n-1}(x) - T_{n-2}. \]
So, e.g.
\[ T_2(cos(theta)) = cos(2*theta) = 2cos(theta)^2 - 1 = 2cos(theta)T_1(cos(theta)) - T_0(cos(theta)) \]
2. Computers actually use this recursion to calculate sines and cosines! So, it's a little inefficient to code your Chebyshev polynomials using `math.sin`.
3. Using generating functions, you can get a closed form for T_n(x) that only takes O(log n) time to calculate. (Note: assuming you count multiplications as constant. However, you actually need O(log n) bits to accurately represent x, so it's more accurately O((log n)^2 log log n).)