3

I am using Python 2 and the fairly simple method given in Wikipedia's article "Cubic function". This could also be a problem with the cube root function I have to define in order to create the function mentioned in the title.

# Cube root and cubic equation solver
#
# Copyright (c) 2013 user2330618
#
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, you can obtain one at http://www.mozilla.org/MPL/2.0/.

from __future__ import division
import cmath
from cmath import log, sqrt

def cbrt(x):
    """Computes the cube root of a number."""
    if x.imag != 0:
        return cmath.exp(log(x) / 3)
    else:
        if x < 0:
            d = (-x) ** (1 / 3)
            return -d
        elif x >= 0:
            return x ** (1 / 3)

def cubic(a, b, c, d):
    """Returns the real roots to cubic equations in expanded form."""
    # Define the discriminants
    D = (18 * a * b * c * d) - (4 * (b ** 3) * d) + ((b ** 2) * (c ** 2)) - \
    (4 * a * (c ** 3)) - (27 * (a ** 2) * d ** 2)
    D0 = (b ** 2) - (3 * a * c)
    i = 1j  # Because I prefer i over j
    # Test for some special cases
    if D == 0 and D0 == 0:
        return -(b / (3 * a))
    elif D == 0 and D0 != 0:
        return [((b * c) - (9 * a * d)) / (-2 * D0), ((b ** 3) - (4 * a * b * c)
        + (9 * (a ** 2) * d)) / (-a * D0)]
        else:
            D1 = (2 * (b ** 3)) - (9 * a * b * c) + (27 * (a ** 2) * d)
            # More special cases
            if D != 0 and D0 == 0 and D1 < 0:
                C = cbrt((D1 - sqrt((D1 ** 2) - (4 * (D0 ** 3)))) / 2)
            else:
                C = cbrt((D1 + sqrt((D1 ** 2) - (4 * (D0 ** 3)))) / 2)
                u_2 = (-1 + (i * sqrt(3))) / 2
                u_3 = (-1 - (i * sqrt(3))) / 2
                x_1 = (-(b + C + (D0 / C))) / (3 * a)
                x_2 = (-(b + (u_2 * C) + (D0 / (u_2 * C)))) / (3 * a)
                x_3 = (-(b + (u_3 * C) + (D0 / (u_3 * C)))) / (3 * a)
                if D > 0:
                    return [x_1, x_2, x_3]
                else:
                    return x_1

I've found that this function is capable of solving some simple cubic equations:

print cubic(1, 3, 3, 1)
-1.0

And a while ago I had gotten it to a point where it could solve equations with two roots. I've just done a rewrite and now it's gone haywire. For example, these coefficients are the expanded form of (2x - 4)(x + 4)(x + 2) and it should return [4.0, -4.0, -2.0] or something similar:

print cubic(2, 8, -8, -32)
[(-4+1.4802973661668753e-16j), (2+2.9605947323337506e-16j), (-2.0000000000000004-1.1842378929335002e-15j)]

Is this more a mathematical or a programming mistake I'm making?

Update: Thank you, everyone, for your answers, but there are more problems with this function than I've iterated so far. For example, I often get an error relating to the cube root function:

print cubic(1, 2, 3, 4)  # Correct solution: about -1.65
...
    if x > 0:
TypeError: no ordering relation is defined for complex numbers
print cubic(1, -3, -3, -1)  # Correct solution: about 3.8473
    if x > 0:
TypeError: no ordering relation is defined for complex numbers
martineau
  • 119,623
  • 25
  • 170
  • 301
  • 11
    The values you show are only infinitesimally different from what they should be (they are real numbers with tiny tiny imaginary parts). It is possibly due to a rounding error or some floating point imprecision in the calculations. – BrenBarn Apr 29 '13 at 03:59
  • 1
    What's with the special cases in the cube root function? I think I see a bug caused by it: it won't work if `x.imag` is negative. – icktoofay Apr 29 '13 at 04:06
  • 3
    If you are the author of the code then it's no problem to publish it somewhere else with a totally different license (more permissive or more restrictive, even though the latter makes little sense). However, you can't re-post other people's code published under a more restrictive license here. – Florian Brucker Apr 29 '13 at 04:06
  • 1
    «This site's terms of service do not allow you to use any content posted here for commercial purposes,» That's [simply not true](http://creativecommons.org/licenses/by-sa/3.0/). It's also completely off-topic. If you want to get help here, then do so without trying to end-run around the terms. Questions here should be self-contained, containing all the code that's necessary to describe the problem such that it can be answered. Please edit your question to include your code; you do not lose the rights to it, and you may use any code posted in an answer for any purpose, with attribution. – jscs Apr 29 '13 at 05:26
  • "Any other downloading, copying, or storing any Content for other than personal, noncommercial use is expressly prohibited without prior written permission from Stack Exchange, or from the copyright holder identified in such Content's copyright notice." I am aware that I will not lose my rights to the code, but I'll post the code here and give such permission. –  Apr 29 '13 at 14:36
  • That paragraph covers stuff that's copyright of Stack Exchange Inc, which is called "Content". You need to look at the terms for "Subscriber Content", the questions and answers posted by users. I am not a lawyer; this is not legal advice, etc., etc. – jscs Apr 29 '13 at 17:13
  • Nevertheless, code snippets on this side can only be used under a license that Creative Commons has explicitly discouraged for software. –  Apr 30 '13 at 04:04

3 Answers3

7

Wolfram Alpha confirms that the roots to your last cubic are indeed

(-4, -2, 2)

and not as you say

... it should return [4.0, -4.0, -2.0]

Not withstanding that (I presume) typo, your program gives

[(-4+1.4802973661668753e-16j), (2+2.9605947323337506e-16j), (-2.0000000000000004-1.1842378929335002e-15j)]

Which to accuracy of 10**(-15) are the exact same roots as the correct solution. The tiny imaginary part is probably due, as others have said, to rounding.

Note that you'll have to use exact arithmetic to always correctly cancel if you are using a solution like Cardano's. This one of the reasons why programs like MAPLE or Mathematica exist, there is often a disconnect from the formula to the implementation.

To get only the real portion of a number in pure python you call .real. Example:

a = 3.0+4.0j
print a.real
>> 3.0
Hooked
  • 84,485
  • 43
  • 192
  • 261
  • That was not a typo. Python returns floats, and I've seen no need thus far to return integers. In any case, how can I filter out that tiny imaginary part? –  Apr 29 '13 at 14:45
  • 3
    No the typo was that you wrote the WRONG ROOTS of (4,-4,2) not (-4,-2,2). I'll edit the answer to show how to pull out only the real portion. – Hooked Apr 29 '13 at 14:48
  • Kindly look again at what you posted. I wrote [4.0, -4.0, -2.0], the same roots but with 4 moved to the beginning. In any case, I've set the function to return only the real part of the roots, and now I'm getting [-4.0, 2.0, -2.0000000000000004]. –  Apr 29 '13 at 23:43
  • 1
    @user2330618 [4,-4,-2] is not the same as [-4,-2,2] there is multiplicity on the two not the four. It looks like the fix solved your original problem. As for your update, I would ask a _new question_ and link to this one. This site is not a discussion board, but a Q&A (check out the FAQ!). We will be happy to help you from there. – Hooked Apr 30 '13 at 00:59
  • 1
    @user2330618 Side note: it's worth understanding _why_ complex numbers are not well ordered, it's a _really_ important concept. For example "1" and "i" both have unit magnitude, but which ones comes first if I want to "order" them? – Hooked Apr 30 '13 at 01:01
  • I'll ask a new question. Is this relevant to my problem? –  Apr 30 '13 at 03:54
  • @user2330618 Can't tell without the details, please link to the new question when you ask it. If the question is solved please mark it as such. Try to avoid an extended discussion in comments. – Hooked Apr 30 '13 at 04:09
  • [New question here.](http://stackoverflow.com/questions/16291902/why-am-i-receiving-a-no-ordering-relation-defined-for-complex-numbers-error) –  Apr 30 '13 at 14:29
3

Hooked's answer is the way to go if you want to do this numerically. You can also do it symbolically using sympy:

>>> from sympy import roots
>>> roots('2*x**3 + 8*x**2 - 8*x - 32')
{2: 1, -4: 1, -2: 1}

This gives you the roots and their multiplicity.

Community
  • 1
  • 1
Florian Brucker
  • 9,621
  • 3
  • 48
  • 81
  • Are these really just numerical roots or the exact symbolic expressions if the factorization is not so simple? +1 for the python library suggestion! – Hooked Apr 30 '13 at 01:04
  • @Hooked: These are exact symbolic expressions. Take a look at the documentation I've linked to, it contains examples where the roots are (symbolically represented) irrationals. – Florian Brucker Apr 30 '13 at 10:34
-5

You are using integer values - which are not automatically converted to floats by Python. The more generic solution will be to write coefficients in the function as float numbers - 18.0 instead of 18, etc. That will do the trick An illustration - from the code:

>>> 2**(1/3)
1
>>> 2**(1/3.)
1.2599210498948732
>>>
volcano
  • 3,578
  • 21
  • 28
  • 5
    I don't think that would cause the problem. Most operations (adding, subtracting, multiplying) automatically promote the other operand to a float if one is already a float. Division normally wouldn't, but `from __future__ import division` is present, so it will. – icktoofay Apr 29 '13 at 04:04
  • @ictoofay Exactly - actually, 18 is the very first coefficient in the code. If a,b c and d are all integers, the first operation that will cause conversion to float is sqrt – volcano Apr 29 '13 at 04:08
  • 1
    Try adding `from __future__ import division` (as exists in the code) to the top of your test and you'll find that `2 ** (1/3)` yields the same result as `2 ** (1/3.)`. – icktoofay Apr 29 '13 at 04:13
  • Sorry, I didn't see the code because it was posted behind a link. I just naturally assumed that since you obviously didn't try your answer to see if it was correct, you must not have seen the code. – Gabe Apr 29 '13 at 04:13
  • Funny people here at stackoverflow - my post was downvoted by some shmucks at least twice. What a lovely community of brothers :-( – volcano Apr 29 '13 at 04:23
  • 2
    @volcano: though I didn't downvote -- why waste the rep?, is my view -- your answer can't be right, as icktoofay has explained twice. This fits one of the criteria for [when to downvote](http://stackoverflow.com/privileges/vote-down), namely when an answer is clearly incorrect. – DSM Apr 29 '13 at 04:27
  • 2
    @volcano I think you misunderstand the downvote. It doesn't mean that _you_ or your knowledge is incomplete/incorrect. It means that the downvoters think that **this answer** is incomplete or incorrect. To take it personally would be a folly, especially when there is information on _why_ this may be incorrect in the comments. Think of this as constructive criticism. – Hooked Apr 29 '13 at 14:54
  • 1
    I didn't downvote *you*; I downvoted *this answer*. And the reason I did so is that you can verify yourself that it is wrong simply by downloading the sample code, making all the integer literals into floating point, and seeing that those changes don't make a difference. – Gabe Apr 29 '13 at 17:01