c    @(#)cubert.f	1.1
      subroutine cubert(a1,a2,a3,a4,s1,s2,s3,ier)
c
c
c subroutine to solve the equation a1*x**3 + a2*x**2 + a3*x + a4 = 0
c
c
c only real roots are provided.  ier will be positive if there are
c complex roots which are not returned.  ier less than zero indicates
c that the equation has no roots.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      complex a, b, tmp, tmp2
      parameter (two27=2./27., rnegone3=-.33333333,
     .          rnegnin27=(-9./27.), one27=1./27.)
c
c
c
      if(a1.ne.0.0)go to 300
      if(a2.ne.0.0)go to 200
      if(a3.ne.0.0)go to 100
c
      ier=-1
      return
c
c
  100 continue
      ier=0
      s1=-a4/a3
      s2=s1
      s3=s1
      return
c
c
  200 continue
      ier=-1
      r2=a3*a3-4.0*a2*a4
      if(r2.lt.0.0)return
      ier=+2
      ta2=a2+a2
      r=sqrt(r2)/ta2
      t=-a3/ta2
      s1=t+r
      s2=t-r
      s3=s2
      return
c
c
  300 continue
      p=a2/a1
      q=a3/a1
      r=a4/a1
c
      p2=p*p
      as=q+rnegone3*p2
      bs=(two27*p2+rnegnin27*q)*p + r
c
      ta=.25*bs*bs + one27*as*as*as
      b=csqrt(cmplx(ta))
      bsh=bs/2.0
c
      tmp=-cmplx(bsh)+b
      a=0.0
      if(abs(tmp).gt.0.0)a=tmp**.33333333
      if(ta.ge.0.0 .and. real(tmp) .lt. 0.0)a=
     .                         -(abs(real(tmp))**.33333333)
      tmp=-cmplx(bsh)-b
      b=0.0
      if(abs(tmp).gt.0.0)b=tmp**.3333333
      if(ta.ge.0.0 .and. real(tmp) .lt. 0.0)b=
     .                          -(abs(real(tmp))**.33333333)
c
      pone3=.3333333*p
      if(ta.gt.0.0)go to 310
      if(ta.eq.0.0)go to 330
      if(ta.lt.0.0)go to 331
c
  310 continue
      ier=1
      s1=real(a+b)
      tmp=-(a+b)/2.0
      tmp2=(a-b)*csqrt(cmplx(-3.0))/2.0
      s2=real(tmp+tmp2)
      s3=real(tmp-tmp2)
      d2=abs(s1-s2)
      d3=abs(s1-s3)
      d4=abs(s2-s3)
      if(d2.lt.d4)s1=s3
      if(d3.lt.d4)s1=s2
      s1=s1-pone3
      s2=s1
      s3=s1
      return
c
  330 continue
      ier=0
      tmp=((a-b)*csqrt(cmplx(-3.0)))/2.0
      tmp2=-(a+b)/2.0
      s2=real(tmp2+tmp-pone3)
      s3=real(tmp2-tmp-pone3)
      s1=real(a+b)-pone3
      return
c
  331 continue
      ier=0
      s1=real(a+b)-pone3
      tmp2=-(a+b)/2.0
      tmp=((a-b)*csqrt(cmplx(-3.0)))/2.0
      s2=real(tmp2 + tmp -pone3)
      s3=real(tmp2 - tmp -pone3)
      return
      end
