# Cardano's solution of the cubic equation **URL:** \ **Category:** Help\ **Created:** [June 12, 2020, 9:30am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111 "2020-06-12T09:30:32Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![ELNS](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/elns/32/3338_2.png) [@ELNS](https://fortran-lang.discourse.group/u/ELNS)\ **Post date:** [June 12, 2020, 9:30am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/1 "2020-06-12T09:30:32Z")
Hi all I am going to implement exercise below in fortran: ![Exercise_14](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/7e4b290a0d4022985ef29df31e0de8e6d9c545b4.png) 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 ) ---
**Author:** ![stavros](https://avatars.discourse-cdn.com/v4/letter/s/51bf81/32.png) [@stavros](https://fortran-lang.discourse.group/u/stavros)\ **Post date:** [June 12, 2020, 12:07pm UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/2 "2020-06-12T12:07:33Z")
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)` ? ---
**Author:** ![ELNS](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/elns/32/3338_2.png) [@ELNS](https://fortran-lang.discourse.group/u/ELNS)\ **Post date:** [June 13, 2020, 8:26am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/3 "2020-06-13T08:26:04Z")
With sign(1, R), It also prints wrong answer. ---
**Author:** ![ELNS](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/elns/32/3338_2.png) [@ELNS](https://fortran-lang.discourse.group/u/ELNS)\ **Post date:** [June 13, 2020, 9:00am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/4 "2020-06-13T09:00:30Z")
I created an issue for this: > > > 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... ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [June 13, 2020, 10:26am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/5 "2020-06-13T10:26:37Z")
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](http://web.mit.edu/18.06/www/Spring17/Eigenvalue-Polynomials.pdf) or [here](https://www.math.utah.edu/~gustafso/s2016/2270/labs/lab7-polyroot-qrmethod.pdf)) and the standard eigenvalue solver from LAPACK: ```auto 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: ```auto $ 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. ---
**Author:** ![gardhor](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/gardhor/32/88_2.png) [@gardhor](https://fortran-lang.discourse.group/u/gardhor)\ **Post date:** [June 14, 2020, 10:04am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/6 "2020-06-14T10:04:09Z")
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 Q **3-R** 2\<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) ---
**Author:** ![stavros](https://avatars.discourse-cdn.com/v4/letter/s/51bf81/32.png) [@stavros](https://fortran-lang.discourse.group/u/stavros)\ **Post date:** [June 17, 2020, 1:53pm UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/7 "2020-06-17T13:53:36Z")
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: ```auto 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 ``` ```auto 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. ---
**Author:** ![ELNS](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/elns/32/3338_2.png) [@ELNS](https://fortran-lang.discourse.group/u/ELNS)\ **Post date:** [June 18, 2020, 8:46am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/8 "2020-06-18T08:46:40Z")
Thank you for helping me 🙏 🙏 Can i use this in my repository? ---
**Author:** ![ELNS](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/elns/32/3338_2.png) [@ELNS](https://fortran-lang.discourse.group/u/ELNS)\ **Post date:** [June 18, 2020, 9:35am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/9 "2020-06-18T09:35:54Z")
I got this error: m\_cardano.f90:14:34: ``` public t_cubic_solution, solve 1 ``` Error: Symbol ‘solve’ at (1) has no IMPLICIT type ---
**Author:** ![stavros](https://avatars.discourse-cdn.com/v4/letter/s/51bf81/32.png) [@stavros](https://fortran-lang.discourse.group/u/stavros)\ **Post date:** [June 18, 2020, 10:43am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/10 "2020-06-18T10:43:13Z")
Hi @ELNS, it works for me, check here [https://godbolt.org/z/guE64H](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. ---
**Author:** ![ELNS](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/elns/32/3338_2.png) [@ELNS](https://fortran-lang.discourse.group/u/ELNS)\ **Post date:** [June 18, 2020, 1:44pm UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/11 "2020-06-18T13:44:53Z")
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. ---
**Author:** ![Ana91](https://avatars.discourse-cdn.com/v4/letter/a/e9bcb4/32.png) [@Ana91](https://fortran-lang.discourse.group/u/Ana91)\ **Post date:** [May 20, 2022, 2:47am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/12 "2022-05-20T02:47:15Z")
programar en ([https://www.mycompiler.io/](https://www.mycompiler.io/)) ```auto 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 ``` ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [May 20, 2022, 8:51am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/13 "2022-05-20T08:51:30Z")
@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 🙂 ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [May 20, 2022, 10:15am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/14 "2022-05-20T10:15:51Z")
@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](https://user-images.githubusercontent.com/57490926/84569260-689e6000-ad9a-11ea-9f2f-ca803ec006ca.png) that I have verified. ---
**Author:** ![macneacail](https://avatars.discourse-cdn.com/v4/letter/m/a698b9/32.png) [@macneacail](https://fortran-lang.discourse.group/u/macneacail)\ **Post date:** [May 20, 2022, 3:38pm UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/15 "2022-05-20T15:38:09Z")
```auto ! 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. ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [May 21, 2022, 3:38am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/16 "2022-05-21T03:38:13Z")
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](https://link.springer.com/book/10.1007/978-981-13-7114-1), 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!” ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [May 21, 2022, 4:47am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/17 "2022-05-21T04:47:11Z")
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](https://math.stackexchange.com/questions/1306447/what-is-the-principal-cubic-root-of-8). ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [May 21, 2022, 1:03pm UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/18 "2022-05-21T13:03:05Z")
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](https://arxiv.org/pdf/1203.1034.pdf), 2012. Their complex polynomial solver, [cmplx\_roots\_sg](http://www.astrouw.edu.pl/~jskowron/cmplx_roots_sg/), can be applied to the polynomial example of this thread using the following driver: ```auto 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 ``` ---
**Author:** ![Harper](https://avatars.discourse-cdn.com/v4/letter/h/b5ac83/32.png) [@Harper](https://fortran-lang.discourse.group/u/Harper)\ **Post date:** [May 22, 2022, 9:06am UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/19 "2022-05-22T09:06:56Z")
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. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [May 22, 2022, 4:13pm UTC](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/20 "2022-05-22T16:13:22Z")
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. [![](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/b/bad45aa62eeee39f1f80414362197c07d561f9fa.jpeg "How Imaginary Numbers Were Invented") ](https://www.youtube.com/watch?v=cUzklzVXJwo) [Next page](https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111.md?page=2)