Cardano's solution of the cubic equation

Hi all
I am going to implement exercise below in fortran:

i wrote code below :

program main
    implicit none
    integer :: k
    real :: a_1, a_2, a_3, Q, R, discriminant, theta, s, x_1
    complex :: z_2, z_3, omega
    real , dimension(3) :: x_k
    real, parameter :: pi = 3.1415
    print *, "Please enter coefficients for cubic equation:"
    print *, "Cubic equation -> x^3 + a_1 * x^2 + a_2 * x + a_3 = 0"
    read *, a_1, a_2, a_3
    Q = (a_1 ** 2 - 3 * a_2) / 9
    R = ((2 * a_1 ** 3) - 9 * a_1 * a_2 + 27 * a_3) / 54
    discriminant = Q ** 3 - R ** 2
    if (discriminant >= 0) then
        theta = acos(R / sqrt(Q ** 3))
        do k = 1, 3
            x_k(k) = (-2 * sqrt(Q)) * cos((theta + 2 * pi * (k - 1)) / 3) - (a_1 / 3)
        end do
        print *, "Real roots are:", x_k
    else
        s = -1 * sign(R, R) * (abs(R) + sqrt(R ** 2 - Q ** 3))
        omega = (-1, 1.7320) / 2
        x_1 = s + Q / s - a_1 / 3
        z_2 = (omega * s) + (omega ** 2 * (Q / s)) - a_1 / 3
        z_3 = (omega ** 2 * s) + (omega * Q / s) - a_1 / 3
        print *, "Real root is:", x_1, "Complex roots are:", z_2, z_3
    end if
 end program main

For a_1 = 1, a_2 = 2, a_3 = 3 ,
it should print
x_1 = -1.27568
z_2 = 0.13784 + i * 1.52731
z_3 = 0.13784 - i * 1.52731
But it prints
x_1 = -3.12811947
z_2 = ( 1.06406796 , -2.74305439 )
z_3 = ( 1.06392860 , 2.74305439 )

1 Like

Hi @ELNS, I do not know where the problem is, but one error might be at sign(R,R), this makes no sense, since sign(R,R)= R. The intrinsic function sign(A,B) returns A with the sign of B. eg sign(3,-5) = -3. Perhaps what you need is sign(1,R) ?

3 Likes

With sign(1, R), It also prints wrong answer.

1 Like

I created an issue for this:

2 Likes

I could not find the error, but I had a go at adapting your program to any order polynomial using the polynomial companion matrix (for the theory click here or here) and the standard eigenvalue solver from LAPACK:

program main
    implicit none
    integer, parameter :: sp = kind(1.0)
    integer, parameter :: n = 3
    real(sp) :: a(n), M(n,n),wr(n),wi(n)
    real(sp) :: work(3*n), vl(n),vr(n)
    integer :: k, info
    character(len=20) :: num

    if (command_argument_count() == n) then
        do k = 1, n
            call get_command_argument(k,num)
            read(num,*) a(k)
        end do
        print *, "Coefficients: ",a
    else
        write(*,"(A,I0,A)") "Please enter coefficients for the ", n, "-th order polynomial"
        write(*,"(A)") "(the highest order coefficient is assumed to be 1)"
        read(*,*) a(:)
    end if
    print *, "Polynomial coefficients (a_1,a_2,...):", a

    ! Form the companion matrix
    M = 0
    do k = 1, n-1
        M(k,k+1) = 1
    end do
    do k = n,1,-1
        M(n,n-k+1) = -a(k)
    end do

    print *, "Polynomial companion matrix:"
    do k = 1, n
        print *, M(k,:)
    end do

    ! Call LAPACK routine xGEEV
    call sgeev('N','N',n,M,n,wr,wi,vl,1,vr,1,work,3*n,info)
    if (info /= 0) then
        print *, "[sgeev] info = ", info
        if (info < 0) then
            print *, "The ",abs(info),"-th argument had an illegal value."
        else
            print *, "The QR algorithm failed to compute all the eigenvalues."
        end if
        error stop 1
    end if

    print *, "Real and Imaginary parts of the roots:"
    do k = 1, n
        print *, wr(k), wi(k)
    end do

 end program main

Example of running the program:

$ gfortran -Wall main.f90 -o main -llapack
$ ./main 1 2 3
 Coefficients:    1.00000000       2.00000000       3.00000000    
 Polynomial coefficients (a_1,a_2,...):   1.00000000       2.00000000       3.00000000    
 Polynomial companion matrix:
   0.00000000       1.00000000       0.00000000    
   0.00000000       0.00000000       1.00000000    
  -3.00000000      -2.00000000      -1.00000000    
 Real and Imaginary parts of the roots:
  -1.27568197       0.00000000    
  0.137840927       1.52731216    
  0.137840927      -1.52731216    

For higher order polynomials, fourth, fifth, sixth, … the value of the parameter n should be increased and the program should be recompiled.

4 Likes

Hi,
The problem is not (completely) the fortran code, it is the exercise !!
Take the following parameters: a1=a2=0 and a3=8.
The three solutions are the three cubic roots of -8 so proportional to 2.

With those, parameters:
Q=0, R=a3/2=4 and Q3-R2<0
Then s=-(4+4)=-8

With the exercise expression, the first root, x1=s+Q/s-a1/3=s=-8

Obviously, x1 is not the root of the polynomial.

=> I’m guessing, instead of s you have to use s**(1/3)
and correct sign(R,R) in sign(1,R) (as mentioned)

And now, with your parameters:
Real root is: -1.27568233 Complex roots are: (0.137859225,-1.52726758) (0.137781680,1.52726758)

4 Likes

Since this is still an open issue in your Github, I tried the alternative Cardano’s formula that you posted and it worked, here is one possible implementation:

module m_cardano
implicit none
!Cardano's solutions:
! x1 = s + t - b/(3a) ,  real root
! x2 = -(s+t)/2 - b/(3a) + i(sqrt(3)/2)(s-t) , first complex root
! x3 = -(s+t)/2 - b/(3a) - i(sqrt(3)/2)(s-t) , second complex root
!where:
! s = ( r + sqrt(q**3 + r**2) )**(1/3)
! t = ( r - sqrt(q**3 + r**2) )**(1/3)
! q = (3ac - b**2)/(9a**2)
! r = (9abc - 27da**2 - 2b**3)/(54a**3)
private 
public t_cubic_solution, solve


type t_cubic_solution
    real :: x1
    complex :: x2
    complex :: x3
end type t_cubic_solution


contains
    pure type(t_cubic_solution) function solve(a, b, c, d) result(res)
        real, intent(in)    :: a, b, c, d
        real                :: s, t, q, r, rpart, ipart, temp

        q       = (3.0*a*c - b**2)/(9.0*a**2) 
        r       = (9.0*a*b*c - 27.0*d*a**2 - 2.0*b**3)/(54.0*a**3)

        temp    = r + sqrt(q**3 + r**2)
        s       = sign(1.0, temp) * abs(temp)**(1.0/3.0)

        temp    =  r - sqrt(q**3 + r**2)
        t       = sign(1.0, temp) * abs(temp)**(1.0/3.0)


        res%x1  = s + t - b/(3.0*a)

        rpart   = -(s+t)/2.0 - b/(3.0*a)
        ipart   = (sqrt(3.0)/2.0)*(s-t) 

        res%x2  = cmplx( rpart, ipart )
        res%x3  = cmplx( rpart, -ipart )
    end function solve
end module m_cardano
program cubic_polynomial
use m_cardano
implicit none
!cubic polynomial: P: ax**3 + bx**2 + cx + d = 0
! given polynomial:
! x**3 + x**2 + 2x + 3 = 0
real, parameter :: a = 1.0, b = 1.0, c = 2.0, d = 3.0
character(len=*), parameter :: fmt = '(f8.5, 2(1x,f8.5,sp,f8.5,"i"))'
type(t_cubic_solution) :: p

p = solve(a, b, c, d)

write(*,fmt)p%x1, p%x2, p%x3
end program cubic_polynomial

Note, that inside the function solve you still need to separate the temp variable from its sign to avoid raising a negative number to a real number, which happens in this example when calculating the parameter t, otherwise you would get a NaN value.

2 Likes

Thank you for helping me :pray: :pray:
Can i use this in my repository?

1 Like

I got this error:
m_cardano.f90:14:34:

 public t_cubic_solution, solve
                              1

Error: Symbol ‘solve’ at (1) has no IMPLICIT type

1 Like

Hi @ELNS, it works for me, check here https://godbolt.org/z/guE64H.
The compiler does indeed throw a Note about ABI break from gcc4 but the result is correct.

Of course, you can use it as you like, but keep in mind that this is not a general solution, it only works when the discriminant is greater than zero (one real and two complex roots). It needs further refactoring to be complete. When I find some time I will work on this further.

3 Likes

So i use this in my repository as it is , feel free to submit a PR when you find some time in order to complete it.
Thank you very much.

1 Like

programar en (https://www.mycompiler.io/)

program cubicequation

    implicit none
    
    integer (8):: k
    real (8) :: a_1, a_2, a_3, Q, R, discriminante, theta, s, x_1
    complex :: z_2, z_3, omega
    real(8) , dimension(3) :: x_k
    real(8), parameter :: pi = 3.1415
    write(*,*) "Insira os coeficientes para a equação cúbica:"
    write(*,*) "Equação cúbica -> x^3 + a_1 * x^2 + a_2 * x + a_3 = 0"
    read(*,*) a_1,a_2,a_3
    Q = (a_1 ** 2.0 - 3.0 * a_2) / 9.0
    R = ((2.0 * a_1 ** 3.0) - 9.0 * a_1 * a_2 + 27.0 * a_3) / 54.0
    discriminante = Q * 3.0 - R * 2.0
    if (discriminante >= 0) then
        theta = acos(R / sqrt(Q ** 3.0))
        do k = 1, 3
            x_k(k) = (-2.0 * sqrt(Q)) * cos((theta + 2.0 * pi * (k - 1.0)) / 3.0) - (a_1 / 3.0)
        end do
        write(*,*) "As raízes reais são:", x_k
    else
        s = -1.0 * sign(R, R) * (abs(R) + sqrt(R * 2.0 - Q * 3.0))
        omega = (-1, 1.7320) / 2.0
        x_1 = s + Q / s - a_1 / 3.0
        z_2 = (omega * s) + (omega ** 2.0 * (Q / s)) - a_1 / 3.0
        z_3 = (omega ** 2 * s) + (omega * Q / s) - a_1 / 3.0
        write(*,*) "A verdadeira raiz é:", x_1, "As raízes complexas são:", z_2, z_3
    end if
 end program cubicequation
1 Like

@Ana91, welcome to the forum. The code you posted looks like a straightforward implementation of Cardano’s formulae. Just a few remarks, if you allow me ;).

  • You are mixing double-precision and single-precision numbers variables. Usually you will want to avoid that.
  • The value for pi is inaccurate: only 4 decimals and a single-precision literal at that which is assigned to a double-precision real.
  • The use of kind=8 is non-portable. Not all compilers actually this scheme, though it is common ;).
  • Integer powers are best written in that way: Q ** 3 instead of Q ** 3.0
  • Numbers likes 3.0 are single precision, use 3.0d0 or 3.0_8 (or better: 3.0_dp, where dp is an integer parameter that indicates double precision in a platform-independent way. Lots of possibilities there)

Not trying to be pedantic, just some friendly comments :slight_smile:

1 Like

@Ana91 , there are several things wrong with the program code that you posted, and @Arjen has pointed out the Fortran related errors. It is troubling that you merely posted non-working code, with no comments or questions!

I wish to point out some even more important errors, namely, many of the formulae are wrong or were coded incorrectly. With wrong formulae, no programming language is likely to yield a working program, except by sheer accident. Where did you find the incorrect formulae? Do you have a reference? For instance, the formula for the discriminant is wrong. Next, it is impossible to apply Cardano’s rule without being able to evaluate the cube roots of real numbers, and I do not see any “CBRT” or “**(1.0/3)” in your code.

Here is a link to a solution that I have verified.

1 Like
!  Cardino.f90 
!
!  FUNCTIONS:
!  Cardino - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Cardino
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Cardino
    use Base

    implicit none
    
    COMPLEX(KIND=DP) A
    
    COMPLEX(KIND=DP)  B
    COMPLEX(KIND=DP)  C
    COMPLEX(KIND=DP)  D
    
    ! Set to unlikely values.
    
    COMPLEX(KIND=DP) R 
    COMPLEX(KIND=DP) Q 
    COMPLEX(KIND=DP) S 
    COMPLEX(KIND=DP) SD 
    COMPLEX(KIND=DP) T 
    COMPLEX(KIND=DP) i 
    COMPLEX(KIND=DP) x1 
    COMPLEX(KIND=DP) x2
    COMPLEX(KIND=DP) x3 
    
    
    A = (1.0d0,0.0d0)
    B = (1.0d0,0.0d0)
    C = (1.0d0,0.0d0)
    D = (1.0d0,0.0d0)
    i = (0.0d0,1.0d0)
    !COMPLEX(KIND=DP) :: QQ = (0.0d0,0.0d0)
    
   
    
    Q = ((THREE*a*c) - (b*b))/(NINE*a*a)
    R = ((NINE*a*b*c) - (TWENTYSEVEN*a*a*d) - (TWO*b*b*b))/(TWENTYSEVEN*TWO*a*a*a)
    
    SD = sqrt((Q*Q*Q) + (R*R))
    S = (SD + R)**THIRDROOT
    
    T = (R - SD)**THIRDROOT
    
    x1 = S+T - (b/(THREE*a))
    x2 = (-((S+T)/TWO)) -(b/(THREE*a)) + ((i*(sqrt(THREE)/TWO))*(S-T))
    
    x3 = (-((S+T)/TWO)) -(b/(THREE*a)) - ((i*(sqrt(THREE)/TWO))*(S-T))
    
    ! Variables

    ! Body of Cardino
    print *, 'Cardinao Solver'
    
    write(*,100)a%RE,a%IM, b%RE, b%IM,c%RE, c%IM, d%RE,d%IM ,Q%RE,Q%IM, R%RE, R%IM, S%RE, S%IM,SD%RE, SD%IM ,T%RE, T%IM,x1%RE, x1%IM,x2%RE, x2%IM,x3%RE, x3%IM
100 Format( "  a  Real  :: ",F10.4,"  a  Imaginary  :: ",F10.4,/ & 
            "  b  Real  :: ",F10.4,"  b  Imaginary  :: ",F10.4,/& 
            "  c  Real  :: ",F10.4,"  c  Imaginary  :: ",F10.4,/& 
            "  d  Real  :: ",F10.4,"  d  Imaginary  :: ",F10.4,/& 
            "  Q  Real  :: ",F10.4,"  Q  Imaginary  :: ",F10.4,/& 
            "  R  Real  :: ",F10.4,"  R  Imaginary  :: ",F10.4,/& 
            "  S  Real  :: ",F10.4,"  S  Imaginary  :: ",F10.4,/& 
            "  SD Real  :: ",F10.4,"  SD Imaginary  :: ",F10.4,/& 
            "  T  Real  :: ",F10.4,"  T  Imaginary  :: ",F10.4,/& 
            "  x1 Real  :: ",F10.4,"  x1 Imaginary  :: ",F10.4,/& 
            "  x2 Real  :: ",F10.4,"  x2 Imaginary  :: ",F10.4,/& 
            "  x3 Real  :: ",F10.4,"  x3 Imaginary  :: ",F10.4,/) 

    end program Cardino

I do not include base it is trivial. If I have made a mistake then some one will tell me, the beauty of QA.

It appears to me that many of these students have been afflicted by a typographical error in the textbook page of formulae that they are programming. There is no complete reference or attribution, so we do not know which book this page came from – it is provided as a scanned image.

I believe that the equation

s = -{\rm sgn}(R)(|R|+\sqrt{R^2-Q^3})

should have been

s = -{\rm sgn}(R)(|R|+\sqrt{R^2-Q^3})^\frac{1}{3}

That is, the (1/3) exponent at the end is missing. @gardhor and @stavros spotted this error, but their remarks may not have been noticed by the authors of subsequent posts.

P.S.: More searching allowed me to locate the typo in the following text book, see numbered page 53.

Numerical Methods of Mathematics Implemented in Fortran, S.K.Bose, Springer, 2019.

I sent the author and type-setter, aged 84, an email pointing out the typo, and he promptly replied:

“Thanks for your note. It is greatly appreciated that students are savouring the nuances of numerical math where errors are constant companions!”

You made all the variables complex, and as a result the program will not work for some values of the coefficients of the cubic. For example, take the case chosen, in which a = 1, b = 1, c = 2, d = 3, all real. For this case, your program outputs three complex roots, which is impossible.

When you make Z = cmplx(-27,0), Z**(1.0/3) may not come out as -3, but as a complex number, namely, the principal value of the root, 3(1+i\sqrt 3)/2. When working with complex variables, we have to make sure that the correct pick is made from the multiple values that a single expression can take. See, for example, a related post in Math Stackexchange.

Solving cubic equations is a problem that was solved hundreds of years ago, but we have seen in this thread that implementing the methods in a program may turn out to be non-trivial. Here are a couple of references that help us see how solving for polynomial roots can be an important topic and a goal of current research in applications such as control theory, lens design, etc.

A paper, General Complex Polynomial Root Solver and Its Further Optimization for Binary Microlenses, by J. Skowron and A. Gould, Archive-org, 2012.

Their complex polynomial solver, cmplx_roots_sg, can be applied to the polynomial example of this thread using the following driver:

program solvecubic
implicit none
integer, parameter :: rk=kind(0d0)
complex(rk) :: roots(3),poly(4)
integer :: degree = 3
logical :: polish_roots_after = .false., use_roots_as_starting_points = .false.
poly = [(3d0,0d0),(2d0,0d0),(1d0,0d0),(1d0,0d0)]
call cmplx_roots_gen(roots, poly, degree, polish_roots_after, &
   use_roots_as_starting_points)
print *,roots
end program

mecej4 is quite right on getting the wrong cube root of something. When I tried coding Cardano’s method I had to do the part of the positive-discriminant case that involved a 1.0/3 power in real arithmetic, e.g. if w was complex with zero imaginary part, and its real part was rw, then w**(1.0/3) could have argument (in the complex number sense not the Fortran sense) of pi/3 but rw**(1.0/3)
would not compile if rw<0. However sign(abs(rw)**1.0/3, rw) would be what was needed.

Derek Muller from the Veritasium Youtube channel has made a great video on Cardano and how imaginary numbers were invented. Learning about the history of imaginary numbers, and how they were invented to solve the cubic equation made so much sense.