I've been coding a ray tracing toolkit that does a lot of quadratics, to find the intersections of rays with spherical surfaces. If the roots are pathologically closely spaced, computing the discriminant (b**2-4ac) becomes a problem, but if the roots are very different in magnitude (which is the other, and more common, source of problems with quadratics), you just flip the equation around and solve for 1/x (taking sgn(0) = 1)
-b - sgn(b)*sqrt(b**2-4*a*c) x1 = --------------------------- 2*a
2*c x2 = ---------------------------- -b -sgn(b)*sqrt(b**2-4*a*c)The second quadratic equation is one of those obvious things that lots of people don't think of because they learned the first formula in their youth. (I think I found it in Acton about 20 years ago, an event closely followed by the sound of a loud forehead slap.)
The normal equation system for the basis set {1, x, x**2, ....} is so ill-conditioned for high orders that it rapidly becomes numerically indeterminate even if the data are exact. (It's the famous Hilbert matrix.)
One way of seeing this is that the Chebyshev polynomial of order N is bounded between y = +-1, but has a leading coefficient of 2**N.
If you rearrange terms and divide by 2**N, you get
x**N = (some polynomial of order N-2) + E/2**N
where -1 packages, and a neat saturating S32.32 math package for the 68K, as useful as
I usually use either singular value decomposition or just the ratio of the transforms, followed by windowing in the time and/or frequency domains to control the noise gain. Alternatively you can just smooth off the target function to keep the noise gain reasonable. If you just go for it, the noise gain is what it is.
I mostly agree. But there are a bunch of other pretty things that are actually useful, especially in the relationships between infinite series and contour integrals, and in asymptotic methods such as steepest descents, WKB theory, and higher-order stationary phase.
One of these times I have to try out an FPGA.
Cheers
Phil Hobbs