davideberly / GeometricTools

A collection of source code for computing in the fields of mathematics, geometry, graphics, image analysis and physics.
Boost Software License 1.0
1.14k stars 214 forks source link

Cubic root not found #48

Closed tdewolff closed 1 year ago

tdewolff commented 1 year ago

Solving for the cubic polynomial 16t^3 - 24t^2 + 24t - 8 = 0 should give 0.5 as a real root, but looking at https://github.com/davideberly/GeometricTools/blob/6238789eaf9653e9eec4432ce06821552760c65e/GTL/Mathematics/RootFinders/RootsQuadratic.h#L98 it seems that this case is not covered. For this particular polynomial, the depressed cubic polynomial is x^3 + 0.75x + 0 = 0 (i.e. c1 = 0.75 and c0 = 0), which gets passed to the code for a depressed quadratic where c0 is the c1 from before. Clearly, c0 = 0.75 > 0 which goes to the line above, but that is false!

davideberly commented 1 year ago

I added your example to my unit-test environment for RootsCubic. The solver computes a single real-valued root of 0.5.

void UnitTestRootsCubic::InvestigateRegression()
{
    float p0 = -8.0f, p1 = 24.0f, p2 = -24.0f, p3 = 16.0f;
    std::map<float, size_t> rootMultiplicity{};
    RootsCubic::Solve(p0, p1, p2, p3, rootMultiplicity);
    float maxValue = GetMaxError(p0, p1, p2, p3, rootMultiplicity);
    // rootMultiplicity has a single element (0.5, 1)
    // maxValue = polynomial(0.5), which is computed to be 0
}

It appears the code is working correctly. In the cubic SolveDepressed, c0 is 0 and c1 is 0.75. The execution enters the block "if (c0 == rat0)" and then into "RootsQuadratic::SolveDepressed(c1, rootMultiplicity)". This function has local c0 = 0.75 (c1 from the cubic was passed). The depressed quadratic is x^2+c0 = 0. Because c0 > 0, the depressed quadratic has no real-valued roots. Comments in the code indicate this: "// A complex-conjugate pair of roots.". On return to the cubic SolveDepressed, the "root" is 0.0 (for the depressed cubic). On return to the cubic Solve function, the depressed cubic root is translated back to the original polynomial variable by subtracting q2third = -0.5. This produces the original cubic root of 0.5.

tdewolff commented 1 year ago

I'm so sorry to have wasted your time, I was sure to have erased the issue! The issue in my code was not adding the zero root for the cubic case after having no roots of the quadratic, which I found after re-reading your code carefully. In any case, thank you for this excellent repository and the accompanying papers.

davideberly commented 1 year ago

No worries. I always assume my code is incorrect when people file bugs and will investigate, so you did not waste my time.