Cubic formula gives the wrong result (triple checked)

The name of the pictureThe name of the pictureThe name of the pictureClash Royale CLAN TAG#URR8PPP











up vote
25
down vote

favorite
5












I'd like to solve $ax^3 + bx^2 + cx + d = 0$ using the cubic formula.



I coded three versions of this formula, described in three sources:
MathWorld, EqWorld,

and in the book, "The Unattainable Attempt to Avoid the Casus Irreducibilis for Cubic Equations".



While I get identical results across all versions, these results are incorrect.



For example, for $a=1$, $b=2$, $c=3$, $d=4$,



I find incorrect roots:

$x_1 = -0.1747 - 0.8521i$,

$x_2 = 0.4270 + 1.1995i$,

$x_3 = -2.2523 - 0.3474i$.



The correct roots are:

$x_1 = -1.6506$,

$x_2 = -0.1747 + 1.5469i$,

$x_3 = -0.1747 - 1.5469i$



In case you're interested, the actual code is below.
Thank you for your help!



%% Wolfram version

Q = (3*c - b^2) / 9;
R = (9*b*c - 27*d - 2*b^3) / 54;

D = Q^3 + R^2;
S = (R + sqrt(D))^(1/3);
T = (R - sqrt(D))^(1/3);

x1 = - b/3 + (S + T);
x2 = - b/3 - (S + T) / 2 + sqrt(-3) * (S - T) / 2;
x3 = - b/3 - (S + T) / 2 - sqrt(-3) * (S - T) / 2;

%% Book version

omega1 = - 1/2 + sqrt(-3)/2;
omega2 = - 1/2 - sqrt(-3)/2;

p = (3*a*c - b^2) / (3*a^2);
q = (2*b^3 - 9*a*b*c + 27*(a^3)*d) / (27*a^3);

r = sqrt(q^2/4 + p^3/27);
s = (-q/2 + r)^(1/3);
t = (-q/2 - r)^(1/3);

x1 = s + t - b/(3*a);
x2 = omega1*s + omega2*t - b/(3*a);
x3 = omega2*s + omega1*t - b/(3*a);

%% Eqworld version

p = - 1/3 * (b/a)^2 + (c/a);
q = 2/27 * (b/a)^3 - (b*c)/(3*a^2) + d/a;

D = (p/3)^3 + (q/2)^2;
A = (-q/2 + sqrt(D))^(1/3);
B = (-q/2 - sqrt(D))^(1/3);

y1 = A + B;
y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);

x1 = y1 - b / (3*a);
x2 = y2 - b / (3*a);
x3 = y3 - b / (3*a);









share|cite|improve this question























  • I haven't looked into your solution carefully, but a common error in cases like this is to inadvertently use integer division instead of fractional division, which is an easy mistake in many programming systems. Is it possible that you've done this? Can you check omega1 and make sure the -1/2 didn't get dropped?
    – MJD
    May 17 '17 at 14:32











  • Thank you for your help.This part seems Ok (I'm using Matlab).
    – Akim
    May 17 '17 at 14:44






  • 4




    For purposes of numerical stability, you might be interested in this formulation.
    – J. M. is not a mathematician
    May 17 '17 at 21:02










  • I think this is a great alternative; perhaps one downside is the need to compute trigonometric functions.
    – Akim
    May 18 '17 at 2:15






  • 1




    Akim, it's a tradeoff (which happens often in programming applications): in casus irreducibilis, you have the choice of either avoiding complex arithmetic and using trigonometry, or using only radicals and arithmetic, but won't be able to avoid complex evaluations.
    – J. M. is not a mathematician
    May 18 '17 at 7:13














up vote
25
down vote

favorite
5












I'd like to solve $ax^3 + bx^2 + cx + d = 0$ using the cubic formula.



I coded three versions of this formula, described in three sources:
MathWorld, EqWorld,

and in the book, "The Unattainable Attempt to Avoid the Casus Irreducibilis for Cubic Equations".



While I get identical results across all versions, these results are incorrect.



For example, for $a=1$, $b=2$, $c=3$, $d=4$,



I find incorrect roots:

$x_1 = -0.1747 - 0.8521i$,

$x_2 = 0.4270 + 1.1995i$,

$x_3 = -2.2523 - 0.3474i$.



The correct roots are:

$x_1 = -1.6506$,

$x_2 = -0.1747 + 1.5469i$,

$x_3 = -0.1747 - 1.5469i$



In case you're interested, the actual code is below.
Thank you for your help!



%% Wolfram version

Q = (3*c - b^2) / 9;
R = (9*b*c - 27*d - 2*b^3) / 54;

D = Q^3 + R^2;
S = (R + sqrt(D))^(1/3);
T = (R - sqrt(D))^(1/3);

x1 = - b/3 + (S + T);
x2 = - b/3 - (S + T) / 2 + sqrt(-3) * (S - T) / 2;
x3 = - b/3 - (S + T) / 2 - sqrt(-3) * (S - T) / 2;

%% Book version

omega1 = - 1/2 + sqrt(-3)/2;
omega2 = - 1/2 - sqrt(-3)/2;

p = (3*a*c - b^2) / (3*a^2);
q = (2*b^3 - 9*a*b*c + 27*(a^3)*d) / (27*a^3);

r = sqrt(q^2/4 + p^3/27);
s = (-q/2 + r)^(1/3);
t = (-q/2 - r)^(1/3);

x1 = s + t - b/(3*a);
x2 = omega1*s + omega2*t - b/(3*a);
x3 = omega2*s + omega1*t - b/(3*a);

%% Eqworld version

p = - 1/3 * (b/a)^2 + (c/a);
q = 2/27 * (b/a)^3 - (b*c)/(3*a^2) + d/a;

D = (p/3)^3 + (q/2)^2;
A = (-q/2 + sqrt(D))^(1/3);
B = (-q/2 - sqrt(D))^(1/3);

y1 = A + B;
y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);

x1 = y1 - b / (3*a);
x2 = y2 - b / (3*a);
x3 = y3 - b / (3*a);









share|cite|improve this question























  • I haven't looked into your solution carefully, but a common error in cases like this is to inadvertently use integer division instead of fractional division, which is an easy mistake in many programming systems. Is it possible that you've done this? Can you check omega1 and make sure the -1/2 didn't get dropped?
    – MJD
    May 17 '17 at 14:32











  • Thank you for your help.This part seems Ok (I'm using Matlab).
    – Akim
    May 17 '17 at 14:44






  • 4




    For purposes of numerical stability, you might be interested in this formulation.
    – J. M. is not a mathematician
    May 17 '17 at 21:02










  • I think this is a great alternative; perhaps one downside is the need to compute trigonometric functions.
    – Akim
    May 18 '17 at 2:15






  • 1




    Akim, it's a tradeoff (which happens often in programming applications): in casus irreducibilis, you have the choice of either avoiding complex arithmetic and using trigonometry, or using only radicals and arithmetic, but won't be able to avoid complex evaluations.
    – J. M. is not a mathematician
    May 18 '17 at 7:13












up vote
25
down vote

favorite
5









up vote
25
down vote

favorite
5






5





I'd like to solve $ax^3 + bx^2 + cx + d = 0$ using the cubic formula.



I coded three versions of this formula, described in three sources:
MathWorld, EqWorld,

and in the book, "The Unattainable Attempt to Avoid the Casus Irreducibilis for Cubic Equations".



While I get identical results across all versions, these results are incorrect.



For example, for $a=1$, $b=2$, $c=3$, $d=4$,



I find incorrect roots:

$x_1 = -0.1747 - 0.8521i$,

$x_2 = 0.4270 + 1.1995i$,

$x_3 = -2.2523 - 0.3474i$.



The correct roots are:

$x_1 = -1.6506$,

$x_2 = -0.1747 + 1.5469i$,

$x_3 = -0.1747 - 1.5469i$



In case you're interested, the actual code is below.
Thank you for your help!



%% Wolfram version

Q = (3*c - b^2) / 9;
R = (9*b*c - 27*d - 2*b^3) / 54;

D = Q^3 + R^2;
S = (R + sqrt(D))^(1/3);
T = (R - sqrt(D))^(1/3);

x1 = - b/3 + (S + T);
x2 = - b/3 - (S + T) / 2 + sqrt(-3) * (S - T) / 2;
x3 = - b/3 - (S + T) / 2 - sqrt(-3) * (S - T) / 2;

%% Book version

omega1 = - 1/2 + sqrt(-3)/2;
omega2 = - 1/2 - sqrt(-3)/2;

p = (3*a*c - b^2) / (3*a^2);
q = (2*b^3 - 9*a*b*c + 27*(a^3)*d) / (27*a^3);

r = sqrt(q^2/4 + p^3/27);
s = (-q/2 + r)^(1/3);
t = (-q/2 - r)^(1/3);

x1 = s + t - b/(3*a);
x2 = omega1*s + omega2*t - b/(3*a);
x3 = omega2*s + omega1*t - b/(3*a);

%% Eqworld version

p = - 1/3 * (b/a)^2 + (c/a);
q = 2/27 * (b/a)^3 - (b*c)/(3*a^2) + d/a;

D = (p/3)^3 + (q/2)^2;
A = (-q/2 + sqrt(D))^(1/3);
B = (-q/2 - sqrt(D))^(1/3);

y1 = A + B;
y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);

x1 = y1 - b / (3*a);
x2 = y2 - b / (3*a);
x3 = y3 - b / (3*a);









share|cite|improve this question















I'd like to solve $ax^3 + bx^2 + cx + d = 0$ using the cubic formula.



I coded three versions of this formula, described in three sources:
MathWorld, EqWorld,

and in the book, "The Unattainable Attempt to Avoid the Casus Irreducibilis for Cubic Equations".



While I get identical results across all versions, these results are incorrect.



For example, for $a=1$, $b=2$, $c=3$, $d=4$,



I find incorrect roots:

$x_1 = -0.1747 - 0.8521i$,

$x_2 = 0.4270 + 1.1995i$,

$x_3 = -2.2523 - 0.3474i$.



The correct roots are:

$x_1 = -1.6506$,

$x_2 = -0.1747 + 1.5469i$,

$x_3 = -0.1747 - 1.5469i$



In case you're interested, the actual code is below.
Thank you for your help!



%% Wolfram version

Q = (3*c - b^2) / 9;
R = (9*b*c - 27*d - 2*b^3) / 54;

D = Q^3 + R^2;
S = (R + sqrt(D))^(1/3);
T = (R - sqrt(D))^(1/3);

x1 = - b/3 + (S + T);
x2 = - b/3 - (S + T) / 2 + sqrt(-3) * (S - T) / 2;
x3 = - b/3 - (S + T) / 2 - sqrt(-3) * (S - T) / 2;

%% Book version

omega1 = - 1/2 + sqrt(-3)/2;
omega2 = - 1/2 - sqrt(-3)/2;

p = (3*a*c - b^2) / (3*a^2);
q = (2*b^3 - 9*a*b*c + 27*(a^3)*d) / (27*a^3);

r = sqrt(q^2/4 + p^3/27);
s = (-q/2 + r)^(1/3);
t = (-q/2 - r)^(1/3);

x1 = s + t - b/(3*a);
x2 = omega1*s + omega2*t - b/(3*a);
x3 = omega2*s + omega1*t - b/(3*a);

%% Eqworld version

p = - 1/3 * (b/a)^2 + (c/a);
q = 2/27 * (b/a)^3 - (b*c)/(3*a^2) + d/a;

D = (p/3)^3 + (q/2)^2;
A = (-q/2 + sqrt(D))^(1/3);
B = (-q/2 - sqrt(D))^(1/3);

y1 = A + B;
y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);

x1 = y1 - b / (3*a);
x2 = y2 - b / (3*a);
x3 = y3 - b / (3*a);






algebra-precalculus polynomials roots cubic-equations






share|cite|improve this question















share|cite|improve this question













share|cite|improve this question




share|cite|improve this question








edited May 18 '17 at 20:50

























asked May 17 '17 at 13:55









Akim

38110




38110











  • I haven't looked into your solution carefully, but a common error in cases like this is to inadvertently use integer division instead of fractional division, which is an easy mistake in many programming systems. Is it possible that you've done this? Can you check omega1 and make sure the -1/2 didn't get dropped?
    – MJD
    May 17 '17 at 14:32











  • Thank you for your help.This part seems Ok (I'm using Matlab).
    – Akim
    May 17 '17 at 14:44






  • 4




    For purposes of numerical stability, you might be interested in this formulation.
    – J. M. is not a mathematician
    May 17 '17 at 21:02










  • I think this is a great alternative; perhaps one downside is the need to compute trigonometric functions.
    – Akim
    May 18 '17 at 2:15






  • 1




    Akim, it's a tradeoff (which happens often in programming applications): in casus irreducibilis, you have the choice of either avoiding complex arithmetic and using trigonometry, or using only radicals and arithmetic, but won't be able to avoid complex evaluations.
    – J. M. is not a mathematician
    May 18 '17 at 7:13
















  • I haven't looked into your solution carefully, but a common error in cases like this is to inadvertently use integer division instead of fractional division, which is an easy mistake in many programming systems. Is it possible that you've done this? Can you check omega1 and make sure the -1/2 didn't get dropped?
    – MJD
    May 17 '17 at 14:32











  • Thank you for your help.This part seems Ok (I'm using Matlab).
    – Akim
    May 17 '17 at 14:44






  • 4




    For purposes of numerical stability, you might be interested in this formulation.
    – J. M. is not a mathematician
    May 17 '17 at 21:02










  • I think this is a great alternative; perhaps one downside is the need to compute trigonometric functions.
    – Akim
    May 18 '17 at 2:15






  • 1




    Akim, it's a tradeoff (which happens often in programming applications): in casus irreducibilis, you have the choice of either avoiding complex arithmetic and using trigonometry, or using only radicals and arithmetic, but won't be able to avoid complex evaluations.
    – J. M. is not a mathematician
    May 18 '17 at 7:13















I haven't looked into your solution carefully, but a common error in cases like this is to inadvertently use integer division instead of fractional division, which is an easy mistake in many programming systems. Is it possible that you've done this? Can you check omega1 and make sure the -1/2 didn't get dropped?
– MJD
May 17 '17 at 14:32





I haven't looked into your solution carefully, but a common error in cases like this is to inadvertently use integer division instead of fractional division, which is an easy mistake in many programming systems. Is it possible that you've done this? Can you check omega1 and make sure the -1/2 didn't get dropped?
– MJD
May 17 '17 at 14:32













Thank you for your help.This part seems Ok (I'm using Matlab).
– Akim
May 17 '17 at 14:44




Thank you for your help.This part seems Ok (I'm using Matlab).
– Akim
May 17 '17 at 14:44




4




4




For purposes of numerical stability, you might be interested in this formulation.
– J. M. is not a mathematician
May 17 '17 at 21:02




For purposes of numerical stability, you might be interested in this formulation.
– J. M. is not a mathematician
May 17 '17 at 21:02












I think this is a great alternative; perhaps one downside is the need to compute trigonometric functions.
– Akim
May 18 '17 at 2:15




I think this is a great alternative; perhaps one downside is the need to compute trigonometric functions.
– Akim
May 18 '17 at 2:15




1




1




Akim, it's a tradeoff (which happens often in programming applications): in casus irreducibilis, you have the choice of either avoiding complex arithmetic and using trigonometry, or using only radicals and arithmetic, but won't be able to avoid complex evaluations.
– J. M. is not a mathematician
May 18 '17 at 7:13




Akim, it's a tradeoff (which happens often in programming applications): in casus irreducibilis, you have the choice of either avoiding complex arithmetic and using trigonometry, or using only radicals and arithmetic, but won't be able to avoid complex evaluations.
– J. M. is not a mathematician
May 18 '17 at 7:13










2 Answers
2






active

oldest

votes

















up vote
43
down vote



accepted










The problem occurs when you are taking a cube root: in the Wolfram version, when you compute $$T = sqrt[3]R - sqrtD = sqrt[3]-frac3527-sqrtfrac5027$$ and in corresponding places in the other versions of the cubic formula.



The expression inside the cube root is negative (approximately $-2.66$), and the cubic formula works if you take the real cube root (approximately $-1.39$). But writing $(R - sqrt D)^1/3$ in many computer algebra systems instead computes the principal cube root: the complex root with largest real part. This will be the real cube root for a positive number, but here, it gives us approximately $0.69 - 1.2 i$, and that is the source of your error.



I can't identify the language you're using by looking at it, so I don't know what the best way to avoid this problem is. In Mathematica, there's a built-in CubeRoot command, which will always take the real root of a real number.



You can always make it go with some conditionals: define $T$ to be



  • $(R - sqrtD)^1/3$, if $R - sqrt D ge 0$, and

  • $-(sqrt D - R)^1/3$, if $R - sqrt D < 0$.

(And do a similar thing everywhere else you take a cube root.)




The above will work for real coefficients; for complex coefficients (or even when $D$ is negative), we want to be more careful, because then a real root might not exist.



The idea is that we ran into trouble here, not because we weren't taking the real root necessarily, but because we took two cube roots that "didn't match up with each other". We can rewrite the cubic formula in a few different ways so that we only take one cube root, and express the other cube root we'd have to take in terms of the first. Then the issue doesn't occur.



This gives us a better way to solve the problem in code: after computing $S = (R + sqrt D)^1/3$, set $T = -Q/S$. The justification is this: since $S^3 = R + sqrt D$ and $T^3 = R - sqrt D$, we have $S^3 T^3 = R^2 - D = -Q^3$. Naively, we can cancel the cube roots and assume $ST = -Q$, and this works out.



See also this question and its answer.






share|cite|improve this answer


















  • 1




    Thank you so much for this! I'm using Matlab but ultimately will convert the code to C. Let me check and see how I can fix it.
    – Akim
    May 17 '17 at 14:47






  • 2




    In Matlab, it seems like the nthroot command takes real roots, and C has the cbrt command.
    – Misha Lavrov
    May 17 '17 at 14:50










  • Thanks Misha! I will be able to test this properly in a couple of hours.
    – Akim
    May 17 '17 at 14:52






  • 1




    I agree, it seems most straightforward to proceed with $ST = - Q$. It's interesting that this ambiguity is not mentioned in any of the original sources.
    – Akim
    May 17 '17 at 16:54






  • 1




    @Akim I cannot reproduce the failure of the formula in your example. When I duplicate the Wolfram version of the formula in Mathematica, but replace the definition of $T$ by $T = -Q/S$, I get the correct roots of $x^3 + 20x^2 + 30x + 10 = 0$. Can you double-check that you haven't made a sign error somewhere? After all, your $x_1$ and $x_2$ seem correct, and it's only your $x_3$ that's giving you trouble.
    – Misha Lavrov
    May 18 '17 at 20:31

















up vote
-1
down vote













This does not look right.



y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);


Did you mean cbrt(-3)?






share|cite|improve this answer




















    Your Answer




    StackExchange.ifUsing("editor", function ()
    return StackExchange.using("mathjaxEditing", function ()
    StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
    StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
    );
    );
    , "mathjax-editing");

    StackExchange.ready(function()
    var channelOptions =
    tags: "".split(" "),
    id: "69"
    ;
    initTagRenderer("".split(" "), "".split(" "), channelOptions);

    StackExchange.using("externalEditor", function()
    // Have to fire editor after snippets, if snippets enabled
    if (StackExchange.settings.snippets.snippetsEnabled)
    StackExchange.using("snippets", function()
    createEditor();
    );

    else
    createEditor();

    );

    function createEditor()
    StackExchange.prepareEditor(
    heartbeatType: 'answer',
    convertImagesToLinks: true,
    noModals: false,
    showLowRepImageUploadWarning: true,
    reputationToPostImages: 10,
    bindNavPrevention: true,
    postfix: "",
    noCode: true, onDemand: true,
    discardSelector: ".discard-answer"
    ,immediatelyShowMarkdownHelp:true
    );



    );













     

    draft saved


    draft discarded


















    StackExchange.ready(
    function ()
    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2284996%2fcubic-formula-gives-the-wrong-result-triple-checked%23new-answer', 'question_page');

    );

    Post as a guest






























    2 Answers
    2






    active

    oldest

    votes








    2 Answers
    2






    active

    oldest

    votes









    active

    oldest

    votes






    active

    oldest

    votes








    up vote
    43
    down vote



    accepted










    The problem occurs when you are taking a cube root: in the Wolfram version, when you compute $$T = sqrt[3]R - sqrtD = sqrt[3]-frac3527-sqrtfrac5027$$ and in corresponding places in the other versions of the cubic formula.



    The expression inside the cube root is negative (approximately $-2.66$), and the cubic formula works if you take the real cube root (approximately $-1.39$). But writing $(R - sqrt D)^1/3$ in many computer algebra systems instead computes the principal cube root: the complex root with largest real part. This will be the real cube root for a positive number, but here, it gives us approximately $0.69 - 1.2 i$, and that is the source of your error.



    I can't identify the language you're using by looking at it, so I don't know what the best way to avoid this problem is. In Mathematica, there's a built-in CubeRoot command, which will always take the real root of a real number.



    You can always make it go with some conditionals: define $T$ to be



    • $(R - sqrtD)^1/3$, if $R - sqrt D ge 0$, and

    • $-(sqrt D - R)^1/3$, if $R - sqrt D < 0$.

    (And do a similar thing everywhere else you take a cube root.)




    The above will work for real coefficients; for complex coefficients (or even when $D$ is negative), we want to be more careful, because then a real root might not exist.



    The idea is that we ran into trouble here, not because we weren't taking the real root necessarily, but because we took two cube roots that "didn't match up with each other". We can rewrite the cubic formula in a few different ways so that we only take one cube root, and express the other cube root we'd have to take in terms of the first. Then the issue doesn't occur.



    This gives us a better way to solve the problem in code: after computing $S = (R + sqrt D)^1/3$, set $T = -Q/S$. The justification is this: since $S^3 = R + sqrt D$ and $T^3 = R - sqrt D$, we have $S^3 T^3 = R^2 - D = -Q^3$. Naively, we can cancel the cube roots and assume $ST = -Q$, and this works out.



    See also this question and its answer.






    share|cite|improve this answer


















    • 1




      Thank you so much for this! I'm using Matlab but ultimately will convert the code to C. Let me check and see how I can fix it.
      – Akim
      May 17 '17 at 14:47






    • 2




      In Matlab, it seems like the nthroot command takes real roots, and C has the cbrt command.
      – Misha Lavrov
      May 17 '17 at 14:50










    • Thanks Misha! I will be able to test this properly in a couple of hours.
      – Akim
      May 17 '17 at 14:52






    • 1




      I agree, it seems most straightforward to proceed with $ST = - Q$. It's interesting that this ambiguity is not mentioned in any of the original sources.
      – Akim
      May 17 '17 at 16:54






    • 1




      @Akim I cannot reproduce the failure of the formula in your example. When I duplicate the Wolfram version of the formula in Mathematica, but replace the definition of $T$ by $T = -Q/S$, I get the correct roots of $x^3 + 20x^2 + 30x + 10 = 0$. Can you double-check that you haven't made a sign error somewhere? After all, your $x_1$ and $x_2$ seem correct, and it's only your $x_3$ that's giving you trouble.
      – Misha Lavrov
      May 18 '17 at 20:31














    up vote
    43
    down vote



    accepted










    The problem occurs when you are taking a cube root: in the Wolfram version, when you compute $$T = sqrt[3]R - sqrtD = sqrt[3]-frac3527-sqrtfrac5027$$ and in corresponding places in the other versions of the cubic formula.



    The expression inside the cube root is negative (approximately $-2.66$), and the cubic formula works if you take the real cube root (approximately $-1.39$). But writing $(R - sqrt D)^1/3$ in many computer algebra systems instead computes the principal cube root: the complex root with largest real part. This will be the real cube root for a positive number, but here, it gives us approximately $0.69 - 1.2 i$, and that is the source of your error.



    I can't identify the language you're using by looking at it, so I don't know what the best way to avoid this problem is. In Mathematica, there's a built-in CubeRoot command, which will always take the real root of a real number.



    You can always make it go with some conditionals: define $T$ to be



    • $(R - sqrtD)^1/3$, if $R - sqrt D ge 0$, and

    • $-(sqrt D - R)^1/3$, if $R - sqrt D < 0$.

    (And do a similar thing everywhere else you take a cube root.)




    The above will work for real coefficients; for complex coefficients (or even when $D$ is negative), we want to be more careful, because then a real root might not exist.



    The idea is that we ran into trouble here, not because we weren't taking the real root necessarily, but because we took two cube roots that "didn't match up with each other". We can rewrite the cubic formula in a few different ways so that we only take one cube root, and express the other cube root we'd have to take in terms of the first. Then the issue doesn't occur.



    This gives us a better way to solve the problem in code: after computing $S = (R + sqrt D)^1/3$, set $T = -Q/S$. The justification is this: since $S^3 = R + sqrt D$ and $T^3 = R - sqrt D$, we have $S^3 T^3 = R^2 - D = -Q^3$. Naively, we can cancel the cube roots and assume $ST = -Q$, and this works out.



    See also this question and its answer.






    share|cite|improve this answer


















    • 1




      Thank you so much for this! I'm using Matlab but ultimately will convert the code to C. Let me check and see how I can fix it.
      – Akim
      May 17 '17 at 14:47






    • 2




      In Matlab, it seems like the nthroot command takes real roots, and C has the cbrt command.
      – Misha Lavrov
      May 17 '17 at 14:50










    • Thanks Misha! I will be able to test this properly in a couple of hours.
      – Akim
      May 17 '17 at 14:52






    • 1




      I agree, it seems most straightforward to proceed with $ST = - Q$. It's interesting that this ambiguity is not mentioned in any of the original sources.
      – Akim
      May 17 '17 at 16:54






    • 1




      @Akim I cannot reproduce the failure of the formula in your example. When I duplicate the Wolfram version of the formula in Mathematica, but replace the definition of $T$ by $T = -Q/S$, I get the correct roots of $x^3 + 20x^2 + 30x + 10 = 0$. Can you double-check that you haven't made a sign error somewhere? After all, your $x_1$ and $x_2$ seem correct, and it's only your $x_3$ that's giving you trouble.
      – Misha Lavrov
      May 18 '17 at 20:31












    up vote
    43
    down vote



    accepted







    up vote
    43
    down vote



    accepted






    The problem occurs when you are taking a cube root: in the Wolfram version, when you compute $$T = sqrt[3]R - sqrtD = sqrt[3]-frac3527-sqrtfrac5027$$ and in corresponding places in the other versions of the cubic formula.



    The expression inside the cube root is negative (approximately $-2.66$), and the cubic formula works if you take the real cube root (approximately $-1.39$). But writing $(R - sqrt D)^1/3$ in many computer algebra systems instead computes the principal cube root: the complex root with largest real part. This will be the real cube root for a positive number, but here, it gives us approximately $0.69 - 1.2 i$, and that is the source of your error.



    I can't identify the language you're using by looking at it, so I don't know what the best way to avoid this problem is. In Mathematica, there's a built-in CubeRoot command, which will always take the real root of a real number.



    You can always make it go with some conditionals: define $T$ to be



    • $(R - sqrtD)^1/3$, if $R - sqrt D ge 0$, and

    • $-(sqrt D - R)^1/3$, if $R - sqrt D < 0$.

    (And do a similar thing everywhere else you take a cube root.)




    The above will work for real coefficients; for complex coefficients (or even when $D$ is negative), we want to be more careful, because then a real root might not exist.



    The idea is that we ran into trouble here, not because we weren't taking the real root necessarily, but because we took two cube roots that "didn't match up with each other". We can rewrite the cubic formula in a few different ways so that we only take one cube root, and express the other cube root we'd have to take in terms of the first. Then the issue doesn't occur.



    This gives us a better way to solve the problem in code: after computing $S = (R + sqrt D)^1/3$, set $T = -Q/S$. The justification is this: since $S^3 = R + sqrt D$ and $T^3 = R - sqrt D$, we have $S^3 T^3 = R^2 - D = -Q^3$. Naively, we can cancel the cube roots and assume $ST = -Q$, and this works out.



    See also this question and its answer.






    share|cite|improve this answer














    The problem occurs when you are taking a cube root: in the Wolfram version, when you compute $$T = sqrt[3]R - sqrtD = sqrt[3]-frac3527-sqrtfrac5027$$ and in corresponding places in the other versions of the cubic formula.



    The expression inside the cube root is negative (approximately $-2.66$), and the cubic formula works if you take the real cube root (approximately $-1.39$). But writing $(R - sqrt D)^1/3$ in many computer algebra systems instead computes the principal cube root: the complex root with largest real part. This will be the real cube root for a positive number, but here, it gives us approximately $0.69 - 1.2 i$, and that is the source of your error.



    I can't identify the language you're using by looking at it, so I don't know what the best way to avoid this problem is. In Mathematica, there's a built-in CubeRoot command, which will always take the real root of a real number.



    You can always make it go with some conditionals: define $T$ to be



    • $(R - sqrtD)^1/3$, if $R - sqrt D ge 0$, and

    • $-(sqrt D - R)^1/3$, if $R - sqrt D < 0$.

    (And do a similar thing everywhere else you take a cube root.)




    The above will work for real coefficients; for complex coefficients (or even when $D$ is negative), we want to be more careful, because then a real root might not exist.



    The idea is that we ran into trouble here, not because we weren't taking the real root necessarily, but because we took two cube roots that "didn't match up with each other". We can rewrite the cubic formula in a few different ways so that we only take one cube root, and express the other cube root we'd have to take in terms of the first. Then the issue doesn't occur.



    This gives us a better way to solve the problem in code: after computing $S = (R + sqrt D)^1/3$, set $T = -Q/S$. The justification is this: since $S^3 = R + sqrt D$ and $T^3 = R - sqrt D$, we have $S^3 T^3 = R^2 - D = -Q^3$. Naively, we can cancel the cube roots and assume $ST = -Q$, and this works out.



    See also this question and its answer.







    share|cite|improve this answer














    share|cite|improve this answer



    share|cite|improve this answer








    edited May 17 '17 at 15:23

























    answered May 17 '17 at 14:39









    Misha Lavrov

    38.2k55093




    38.2k55093







    • 1




      Thank you so much for this! I'm using Matlab but ultimately will convert the code to C. Let me check and see how I can fix it.
      – Akim
      May 17 '17 at 14:47






    • 2




      In Matlab, it seems like the nthroot command takes real roots, and C has the cbrt command.
      – Misha Lavrov
      May 17 '17 at 14:50










    • Thanks Misha! I will be able to test this properly in a couple of hours.
      – Akim
      May 17 '17 at 14:52






    • 1




      I agree, it seems most straightforward to proceed with $ST = - Q$. It's interesting that this ambiguity is not mentioned in any of the original sources.
      – Akim
      May 17 '17 at 16:54






    • 1




      @Akim I cannot reproduce the failure of the formula in your example. When I duplicate the Wolfram version of the formula in Mathematica, but replace the definition of $T$ by $T = -Q/S$, I get the correct roots of $x^3 + 20x^2 + 30x + 10 = 0$. Can you double-check that you haven't made a sign error somewhere? After all, your $x_1$ and $x_2$ seem correct, and it's only your $x_3$ that's giving you trouble.
      – Misha Lavrov
      May 18 '17 at 20:31












    • 1




      Thank you so much for this! I'm using Matlab but ultimately will convert the code to C. Let me check and see how I can fix it.
      – Akim
      May 17 '17 at 14:47






    • 2




      In Matlab, it seems like the nthroot command takes real roots, and C has the cbrt command.
      – Misha Lavrov
      May 17 '17 at 14:50










    • Thanks Misha! I will be able to test this properly in a couple of hours.
      – Akim
      May 17 '17 at 14:52






    • 1




      I agree, it seems most straightforward to proceed with $ST = - Q$. It's interesting that this ambiguity is not mentioned in any of the original sources.
      – Akim
      May 17 '17 at 16:54






    • 1




      @Akim I cannot reproduce the failure of the formula in your example. When I duplicate the Wolfram version of the formula in Mathematica, but replace the definition of $T$ by $T = -Q/S$, I get the correct roots of $x^3 + 20x^2 + 30x + 10 = 0$. Can you double-check that you haven't made a sign error somewhere? After all, your $x_1$ and $x_2$ seem correct, and it's only your $x_3$ that's giving you trouble.
      – Misha Lavrov
      May 18 '17 at 20:31







    1




    1




    Thank you so much for this! I'm using Matlab but ultimately will convert the code to C. Let me check and see how I can fix it.
    – Akim
    May 17 '17 at 14:47




    Thank you so much for this! I'm using Matlab but ultimately will convert the code to C. Let me check and see how I can fix it.
    – Akim
    May 17 '17 at 14:47




    2




    2




    In Matlab, it seems like the nthroot command takes real roots, and C has the cbrt command.
    – Misha Lavrov
    May 17 '17 at 14:50




    In Matlab, it seems like the nthroot command takes real roots, and C has the cbrt command.
    – Misha Lavrov
    May 17 '17 at 14:50












    Thanks Misha! I will be able to test this properly in a couple of hours.
    – Akim
    May 17 '17 at 14:52




    Thanks Misha! I will be able to test this properly in a couple of hours.
    – Akim
    May 17 '17 at 14:52




    1




    1




    I agree, it seems most straightforward to proceed with $ST = - Q$. It's interesting that this ambiguity is not mentioned in any of the original sources.
    – Akim
    May 17 '17 at 16:54




    I agree, it seems most straightforward to proceed with $ST = - Q$. It's interesting that this ambiguity is not mentioned in any of the original sources.
    – Akim
    May 17 '17 at 16:54




    1




    1




    @Akim I cannot reproduce the failure of the formula in your example. When I duplicate the Wolfram version of the formula in Mathematica, but replace the definition of $T$ by $T = -Q/S$, I get the correct roots of $x^3 + 20x^2 + 30x + 10 = 0$. Can you double-check that you haven't made a sign error somewhere? After all, your $x_1$ and $x_2$ seem correct, and it's only your $x_3$ that's giving you trouble.
    – Misha Lavrov
    May 18 '17 at 20:31




    @Akim I cannot reproduce the failure of the formula in your example. When I duplicate the Wolfram version of the formula in Mathematica, but replace the definition of $T$ by $T = -Q/S$, I get the correct roots of $x^3 + 20x^2 + 30x + 10 = 0$. Can you double-check that you haven't made a sign error somewhere? After all, your $x_1$ and $x_2$ seem correct, and it's only your $x_3$ that's giving you trouble.
    – Misha Lavrov
    May 18 '17 at 20:31










    up vote
    -1
    down vote













    This does not look right.



    y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
    y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);


    Did you mean cbrt(-3)?






    share|cite|improve this answer
























      up vote
      -1
      down vote













      This does not look right.



      y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
      y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);


      Did you mean cbrt(-3)?






      share|cite|improve this answer






















        up vote
        -1
        down vote










        up vote
        -1
        down vote









        This does not look right.



        y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
        y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);


        Did you mean cbrt(-3)?






        share|cite|improve this answer












        This does not look right.



        y2 = - 1/2 * (A + B) + sqrt(-3)/2 * (A - B);
        y3 = - 1/2 * (A + B) - sqrt(-3)/2 * (A - B);


        Did you mean cbrt(-3)?







        share|cite|improve this answer












        share|cite|improve this answer



        share|cite|improve this answer










        answered Sep 8 at 10:28









        Narf the Mouse

        280148




        280148



























             

            draft saved


            draft discarded















































             


            draft saved


            draft discarded














            StackExchange.ready(
            function ()
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2284996%2fcubic-formula-gives-the-wrong-result-triple-checked%23new-answer', 'question_page');

            );

            Post as a guest













































































            這個網誌中的熱門文章

            How to combine Bézier curves to a surface?

            Carbon dioxide

            Why am i infinitely getting the same tweet with the Twitter Search API?