In the previous post in this series I looked at translating a routine to solve quadratic equations from VBA to Python.  Using the same process I have also created Python code to solve Cubic equations.  The only added difficulty for the Cubic routine was that the VBA function included a call to a VBA function to sort the results in ascending order.  It would have been quite easy to translate the VBA routine, but since Python already includes a sort routine it was easier to use this, so the VBA function call:

BubbleSort CubicA, 1, 3

became in Python:

CubicA[0:3,0:1] = sort(CubicA[0:3,0:1],0)

Note that this is the numpy sort routine, operating on the first 3 rows of a numpy 2D array.

Having generated the code that was returning numerical results in the desired format, it was just necessary to check that the results were correct, however a quick check soon showed that they were not.

For the Quadratic routine a quick check through the code soon found the source of the error.  I had converted the VBA:

If (Abs(B) >= Abs(C)) Then
            E = 1# - (A / B) * (C / B)
            D = Sqr(Abs(E)) * Abs(B)
        Else
            E = A
            If (C < 0#) Then E = -A
            E = B * (B / Abs(C)) - E
            D = Sqr(Abs(E)) * Sqr(Abs(C))
        End If

To Python

     if (abs(B) >= abs(C)):
            E = 1. - (A / B) * (C / B)
            D = ((abs(E)) * abs(B))**0.5
        else:
            E = A
            if (C < 0.):
                 E = -A
             E = B * (B / abs(C)) - E
             D = ((abs(E)) * (abs(C))**0.5)**0.5

The Python math module has a sqrt function, so it would have been simpler to use that, but using the equivalent **0.5, it was just necessary to correct the placement of brackets and exponents:

     if (abs(B) >= abs(C)):
            E = 1. - (A / B) * (C / B)
            D = ((abs(E))**0.5) * abs(B)
        else:
            E = A
            if (C < 0.):
                E = -A
            E = B * (B / abs(C)) - E
            D = ((abs(E))**0.5) * (abs(C))**0.5

The problem with the cubic code was not so obvious, but some experimentation revealed that the problem was occurring when there were one real and two complex roots, and was caused because:

A_ = -sgnR * (abs(R_) + (R2 - Q3) ** 0.5) ** (1 / 3)

was always evaluating to 1 or -1, i.e. to -sgnR.  The problem was due to the fact that in Python version 2 division of two integers always returns an integer, so 1/3 = 0.  To get a decimal result in Python 2 at least one of the two values must be non-integer, so 1./3 = 0.3333333333333333, and yields the correct result.

In this case changing the one instance of 1/3 to 1./3 was sufficient to fix the problem, but a better solution can be achieved by adding:

from __future__ import division

to the top of the code module.  This line causes Python Version 2 to use the Version 3 rules for division, which yields a decimal result when two integers are divide, rather than an integer.

Results from the two functions are shown in the screenshot below:

py_Quad and py_CubicC results