What's wrong with this Python function to solv

问题:

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

回答1:

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

回答2:

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.

回答3:

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
>>>