!23456789 123456789 123456789 123456789 123456789 123456789 123456789 12

      subroutine Lodovico_Ferrari(a,b,c,d,e,  x)

! tested with gnu f95
! finds analitically solutions of quartic ax^4+bx^2+cx^3+dx+e=0
! from  http://en.wikipedia.org/wiki/Quartic_equation

      implicit complex(8) (A-H,O-Z)
      real a,b,c,d,e
      complex(8) x(1:4) ! x is a vector to host four complex solutions

      alpha = -(3.*b**2)/(8.*a**2) + c/a
      beta = b**3./(8.*a**3)-(b*c)/(2.*a**2)+d/a
      gamma = -(3.*b**4)/(256.*a**4) + (c*b**2)/(16.*a**3) -
     &           (b*d)/(4.*a**2) + e/a
      print *, "beta = ",beta
      if(beta.eq.0.) then
!      sign(t)=+
        x(1)=-b/(4*a)+ sqrt( .5*(-alpha + sqrt(alpha**2-4*gamma)) ) 
        x(2)=-b/(4*a)- sqrt( .5*(-alpha + sqrt(alpha**2-4*gamma)) )     
!      sign(t)=-
        x(3)=-b/(4*a)+ sqrt( .5*(-alpha - sqrt(alpha**2-4*gamma)) ) 
        x(4)=-b/(4*a)- sqrt( .5*(-alpha - sqrt(alpha**2-4*gamma)) ) 
       return
      else
       P = -(alpha**2)/12. - gamma
       Q = -(alpha**3)/108.+alpha*gamma/3.-(beta**2)/8.
       R = -Q/2.+sqrt((Q**2)/4.+(P**3)/27.)    
       U = R**(1./3.) 
       if(U.ne.0.) y= -5.*alpha/6. +U -P/(3.*U)
       if(U.eq.0.) y= -5.*alpha/6. -Q**(1./3.)
       W = sqrt(alpha+2.*y)
!      sign(s)=+
        x(1)=-b/(4*a)+.5*(+W+sqrt(-(3*alpha+2.*y+2*beta/W)))
        x(2)=-b/(4*a)+.5*(+W-sqrt(-(3*alpha+2.*y+2*beta/W)))
!      sign(s)=-
        x(3)=-b/(4*a)+.5*(-W+sqrt(-(3*alpha+2.*y-2*beta/W)))
        x(4)=-b/(4*a)+.5*(-W-sqrt(-(3*alpha+2.*y-2*beta/W)))
      endif
      return
      end



      program testFerrari

      logical EOF /.false./
      complex(8) x(1:4) ! x is a vector to host four complex solutions

      do while (.not.EOF)
       print *, "--------------------------------------------------"
       print *, "enter a,b,c,d,e"
       READ (*, *, IOSTAT=IST) a,b,c,d,e
       EOF = IST.LT.0
       IF (.NOT.EOF) THEN
        print *, "a,b,c,d,e: ", a,b,c,d,e
        call Lodovico_Ferrari(a,b,c,d,e,  x)
        print*,"Ferrari returned: ",  (x(i),i=1,4)
       ENDIF
      end do
      end

