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 )

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

2 Likes

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

I created an issue for this:

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.

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

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

1 Like

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

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

 public t_cubic_solution, solve
                              1

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

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.

1 Like

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.