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

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.