Resolve numeric instability with the quartic equation solver

[Lathreas]

The quartic equation solver (which solves equations of the form of a x^4 + b x^3 + c x^2 + d x + e = 0) requires a division by a before any calculation can be made. If a is very small, this division by a small number leads to a numeric instability.

This mainly affects the raytracer, which requires accurate solutions to quartic equations. As a result this issue should be fixed before the raytracer can be considered complete.

This might be resolved by using 'perturbation methods', or even something like the Herbie tool to optimize floating point operations. Alternatively, we could use an internal normalized integer representation (such as C#'s Decimal) to reduce numeric instability overall, if appropriate.

[AdamNorberg]

I wonder if there are opportunities to embed precision-preserving logic in the Real type. I'm starting to think I should first define IReal and then implement DoubleReal and FloatReal as underlying implementations, with GetReal(double) as the canonical factory (that we switch out between Double or Float depending on which we want to use - or perhaps other specialized implementations). One serious downside to this is that a variable of type IReal is a reference, so all IReal math would take indirection costs, boxing and unboxing, etc., even when the underlying implementations are struct. But, I was thinking we might want to make them classes anyway because it's not clear to me that we don't need polymorphic behavior to solve a bunch of actual problems.

I may be badly overthinking this. I specialize in overthinking things.

[AdamNorberg]

Apparently I messed up with a ` somewhere. Oops. Anyway, I'm overthinking this. I'm going to just implement Real as a struct and we can refactor later.

[Lathreas]

(I edited the ` away. :P)

It's certainly an interesting idea, and an automatically optimizing Real certainly does sound interesting. I do feel like it would be too much work to implement that immediately, though, as it's probably easier to refactor the code that has numeric instability rather than make a super smart class that handles all that. And, since we will be dealing with a lot of values, keeping them structs will probably be very important to keep the memory and performance overhead from becoming too much. So yeah, let's indeed make it a simple struct for now.

[AdamNorberg]

WIP: https://github.com/AdamNorberg/anatomy3d/blob/create_real/engine/Real.cs

Remaining: .Pow, .IsNegative, .SigNum (treating negative zero as zero for this purpose), binary arithmetic operators, comparison operators, .IsNear, .IsProportionallyNear, .IsAbsolutelyNear.

[Lathreas]

I've solved the singularities when any of the parameters a, b, c, d, or e become zero by resorting to solving a lower-order equation. See this issue as well.

However, even then some spikes occurred, at locations completely unrelated to the singularity.

I've tracked down the issue to the following point:

In engine/calculus/QuarticFunction.cs, we calculate a few intermediate values to solve the quartic equation. However, for some values of a, b, c, d, or e, p1 and q become very large:

double p1 = c*c*c - 4.5*b*c*d + 13.5*a*d*d + 13.5*b*b*e - 36.0*a*c*e;

double q = c*c - 3.0*b*d + 12.0*a*e;

Complex p2 = p1 + Complex.Sqrt(-q*q*q + p1*p1);

Complex pow = Complex.Pow(p2, (1.0/3.0));

This is likely not due to a singularity, but instead simply because of the many multiplications that occur. This also makes sense, because there was nothing special ('singularity-like') about the positions where the spikes were. I've optimized the code a bit, and converted the intermediate calculations to double. (The intermediates temporarily used single precision to highlight any singularities more easily.) Since the numbers are large but unlikely to explode to infinity, this means I'm quite confident now that double should suffice for all real-world calculations with this function.

To make the claim rigorous that we are not dealing with a singularity, I'll spend a bit longer on finding out the exact details of the specific angle at which it occurs, but for now I consider the numeric instability issues resolved.

[Lathreas]

We should probably attempt to find any additional weak points using extensive unit tests with both large and small values.

[AdamNorberg]

Can you add unit tests that fail when using float for intermediate calculations but pass when using double?

[Lathreas]

Since these issues have not popped up anymore after the fixes, and with the final numeric accuracy PR being reviewed as well, I'll mark this issue as resolved!