Cardano's solution of the cubic equation

Help
5

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