SUBROUTINE harmadfoku(a1, b1, c1, d1, x1, x2, x3)
!Elofeltetel: a1 nem nulla
 IMPLICIT NONE

 DOUBLE PRECISION, INTENT(IN) :: a1, b1, c1, d1 !A bemeno egyutthatok
 DOUBLE PRECISION :: a, b, c !A lenormalt egyutthatok
 DOUBLE PRECISION ::  a2, b2, theta !segedvaltozok
 COMPLEX :: q, r !Segedvaltozok
 DOUBLE PRECISION, PARAMETER :: pi=3.1415926535897932 !A pi szam
 COMPLEX, INTENT(OUT) :: x1, x2, x3 !Gyokok

!Az egyenlet egyutthatoinak atszamitasa, hogy a foegyutthato 1 legyen
    a=b1/a1
    b=c1/a1
    c=d1/a1
!Segedvaltozok
    q=(a**2-3*b)/9
    r=(2*a**3-9*a*b+27*c)/54
   IF ((real(r)**2) .LT. (real(q)**3)) THEN 
!Casus irreducibilis, a valos gyokok csak komplex szamolassal erhetok el
       theta=acos(real(r)/(sqrt((real(q))**3)))
!A komplex gyokvonasnak megfeleloen harom gyok adodik
       x1=-2.0*sqrt(q)*cos(theta/3)-a/3
       x2=-2.0*sqrt(q)*cos((theta+2*pi)/3)-a/3
       x3=-2.0*sqrt(q)*cos((theta-2*pi)/3)-a/3
   ELSE 
!Casus reducibilis, egy valos es ket konjugalt komplex gyok
       a2=0
       IF (real(r) .LT. 0) THEN
          a2=(abs(r)+sqrt(r**2-q**3))**(1.0/3)
       ELSE
          a2=-1*(abs(r)+sqrt(r**2-q**3))**(1.0/3)
       ENDIF
       b2=0
       IF (a2 .NE. 0) THEN
          b2=q/a2
       ENDIF
!A valos gyok
       x1=(a2+b2)-a/3
!A komplex gyokok
       x2=cmplx((-1.0/2)*(a2+b2)-a/3, sqrt(3.0)/2*(a2-b2))
       x3=cmplx((-1.0/2)*(a2+b2)-a/3, -1*sqrt(3.0)/2*(a2-b2))
    ENDIF
 RETURN
 STOP
 END SUBROUTINE harmadfoku