C test cube root function:
C
C
      double complex x1, x2, x3, x4
      double complex theroots(4)

      print *,' '
      print *,'-------------------------------------'
      x1=dcmplx(0.112800000000000d0,0.0d0)
      x2=dcmplx(62.515954868409032d0, -19.496062188142883d0)
      x3=dcmplx(-1.746805372434688d+02 , 9.761793888231617d+01)
      x4=dcmplx(-3.013141689125188d+02 , 2.079155627385262e+02)


      call getcubicroots(x1,x2,x3,x4,theroots)
      
      print *,theroots


      print *,' '
      print *,'-------------------------------------'
      print *,' testing matlabs roots:'
      ! see if matlab's roots work:
      x1=dcmplx(0.112800000000000d0,0.0d0)
      x2=dcmplx(62.515954868409032d0, -19.496062188142883d0)
      x3=dcmplx(-1.746805372434688d+02 , 9.761793888231617d+01)
      x4=dcmplx(-3.013141689125188d+02 , 2.079155627385262e+02)

      theroots(1)=dcmplx(-5.571844753211568d2,1.734645600157649d2)
      theroots(2)=dcmplx(4.244379082127986d0, -0.781521989536298d0)
      theroots(3)=dcmplx(-1.279361814242566d0 , 0.154392719718983d0)

      print *,'xs=',x1,x2,x3,x4
      print *,'roots=',theroots
      call checkroots(x1,x2,x3,x4,theroots)


      print *,' '
      print *,'-------------------------------------'
      print *,' testing SINGLE-PRECISION matlabs roots:'
      ! see if matlab's roots work:
      x1=dcmplx(0.1128d0,0.0d0)
      x2=dcmplx(62.51595d0, -19.4960d0)
      x3=dcmplx(-1.7468d+02 , 9.76179d+01)
      x4=dcmplx(-3.013141d+02 , 2.0792e+02)

      theroots(1)=dcmplx(-5.571844753211568d2,1.734645600157649d2)
      theroots(2)=dcmplx(4.244379082127986d0, -0.781521989536298d0)
      theroots(3)=dcmplx(-1.279361814242566d0 , 0.154392719718983d0)

      print *,'xs=',x1,x2,x3,x4
      print *,'roots=',theroots
      call checkroots(x1,x2,x3,x4,theroots)


      print *,' '
      print *,'-------------------------------------'
      x1 = dcmplx(1.,0.)
      x2 = dcmplx(2.,-10.)
      x3 = dcmplx(-4.,1.)
      x4 = dcmplx(5.,2.)

      call getcubicroots(x1,x2,x3,x4,theroots)
      
      print *,theroots

      stop



      print *,' '
      print *,'-------------------------------------'
      x1 = dcmplx(1.,0.)
      x2 = dcmplx(0.,0.)
      x3 = dcmplx(0.,0.)
      x4 = dcmplx(-8.,0.)

      call getcubicroots(x1,x2,x3,x4,theroots)
      
      print *,theroots


      print *,' '
      print *,'-------------------------------------'
      x1 = dcmplx(1.,0.)
      x2 = dcmplx(-6.,0.)
      x3 = dcmplx(-72.,0.)
      x4 = dcmplx(-27.,0.)

      call getcubicroots(x1,x2,x3,x4,theroots)
      
      print *,theroots



      STOP
      END

      !*********************************************************
      SUBROUTINE getCubicRoots(x1,x2,x3,x4,theroots)

      double COMPLEX x1, x2, x3, x4
      double COMPLEX theroots(4)
    ! coefficients represent the equation:
    ! X1*x^3 + X2*x^2 + X3*x + X4 =0
    ! everything is complex.

      double COMPLEX p,q,r
      double COMPLEX a,b

      double COMPLEX a1, a2, a3, aa, bb
      double COMPLEX apb, amb, scale
      double COMPLEX ya, yb
      double COMPLEX xx1, xx2, xx3

      double complex t1, t2
      
      p = x2/x1
      q = x3/x1
      r = x4/x1

      print *,'p=',p
      print *,'q=',q
      print *,'r=',r

      a = q-p*p/3.
      b = (1./27.)*(2.*p*p*p - 9.*p*q + 27.*r)

      print *,'a=',a,'  b=',b

      !print *,'a=',q,-p*p/3.
      !print *,'b=1/27of:',2.*p*p*p, -9.*p*q, 27.*r
      !print *,'b=1/27of:',2.*p**3, -9.*p*q, 27.*r

      a1 = 0.25*b*b + a*a*a/27.
      a2 = cdsqrt(a1)

      !OK
      !print *,'1 does:',a2*a2,' = ',a1
      
      a3 = -0.5*b
      aa = (a3 + a2)**(1./3.)
      bb = (a3 - a2)**(1./3.)

      print *,' A=',aa,'  B=',bb

      print *,'a3=',a3,'   a2=',a2

      ! debug only:
      t1 = a3+a2
      t2 = a3-a2
      !OK
      !print *,'4 does:',aa**3,' = ', t1
      !OK
      !print *,'5 does:',bb**3,' = ', t2
      

      ! NO
      print *,'2 does:',-3.*aa*bb,' = ',a

      ! OK
      !print *,'3 does:',-aa*aa*aa-bb*bb*bb,' = ', b

      apb = aa+bb
      amb = aa-bb

      scale = CMPLX(0., 0.5*sqrt(3.))

      ya = 0.5*apb
      yb = scale*amb

      xx1 = -p/3. + apb
      xx2 = -p/3. - ya + yb
      xx3 = -p/3. - ya - yb

      theroots(1) = xx1
      theroots(2) = xx2
      theroots(3) = xx3

      CALL checkroots(x1,x2,x3,x4, theroots)
      
      return
      end




      
      SUBROUTINE checkroots(x1,x2,x3,x4,theroots)
      double complex x1,x2,x3,x4,theroots(4)

      double complex y
      double complex aa, a1, a2, a3, a4
      integer i
      do i = 1,3
        y = theroots(i)
        a1 = x1*(y**3)
        a2 = x2*(y**2)
        a3 = x3*y
        a4 =x4
        aa = x1*(y**3) + x2*(y**2) + x3*y +x4
        !print *,'root',i,'=',aa,' as=',a1,a2,a3,a4
        print *,'root',i,'=',aa
      ENDDO

      RETURN
      END
      
