subroutine rootp3(a0,a1,a2,a3,r1,r2,r3,tol,info)

C*******************************************************************************
C
C  rootp3 finds the real root(s) for cubic polynomial 
C      a3 x^3 + a2 x^2 + a1 x + a0 
C  in double precision. This subroutine uses Newton's method to find 
C  the maximum root, then use deflation technique to get the quadratic
C  polynomial. At last this routine uses rootp2 to find other roots.
C  
C  a0,a1,a2,a3: On entry, the coefficients of the cubic.
C  info: Output which indicates how many real roots are found. info takes
C        value 0,1,2,3, or 100. info = 100 means a3=a2=a1=0, so the input
C        coefficients are wrong. info = 0 means a3 = 0, so the cubic is
C        actually a quadratic and no real root is found. info = 2 means a3 = 0,
C        but there are two real roots.
C  r1,r2,r3, Output: The real roots found. some of r3,r2,r1 can be dumb if 
C        info < 3. 
C        
C*******************************************************************************

        implicit double precision (a-h,o-z)

	if( abs(a3) .le. 1.0d-10) then
	    call rootp2(a0,a1,a2,r1,r2,info)
	    return
	else
	    r4 = max( abs(a0/a3), abs(a1/a3)+1.0d0, abs(a2/a3)+1.0d0)
	    a = r4
	    b = -r4

C------ Bisection to get good initial guess ----------------------------------

	do k=1,10
	    rm = 0.5d0*(a+b)
	    if( f(a0,a1,a2,a3,rm)*f(a0,a1,a2,a3,a) .le. 0.0d0) then
	        b = rm
	    else
	        a = rm
	    endif
	enddo

	r4 = 0.5d0*(a+b)	

 10	    r1 = r4 - f(a0,a1,a2,a3,r4)/d(a0,a1,a2,a3,r4)
	    if( abs(r4-r1) .gt. tol .and. 
     1		  abs(f(a0,a1,a2,a3,r4)) .gt. tol) then
	       r4 = r1
	       goto 10
	    else
	       a2 = a3*r1 + a2
	       a1 = a2*r1 + a1
	       call rootp2(a1,a2,a3,r2,r3,info)
	       info = info + 1
	       return
	    endif
	endif

	end

C----------------------------------------------------------------------------

	double precision function f(a0,a1,a2,a3,x)
	implicit double precision (a-h,o-z)
	  f = a0 + x*(a1+x*(a2+a3*x))
	  return
	end

	double precision functiond(a0,a1,a2,a3,x)
        implicit double precision (a-h,o-z)
	   d = a1 + x*(2.0d0*a2 +x*3.0d0*a3)
          return
        end

C------------------------------------------------------------------------------

	subroutine rootp2(a0,a1,a2,r1,r2,info)
        implicit double precision (a-h,o-z)

	if( abs(a2) .le. 1.0d-10) then
	    if( abs(a1) .le. 1.0d-10) then
	      if( abs(a0) .le. 1.0d-10) then
	       info = 10
	      else
	       info = 100
	      endif
	      return
	    else 
	       info = 2
	       r1 = -a0/a1
	       r2 = r1
	       return
	    endif
	else
	    t = a1*a1 - 4.0d0*a0*a2
	    if ( t .ge. 0.0) then
	       info = 2
	       t = dsqrt(t)
	       if( a1 .ge. 0.0d0) then
	          r1 = ( -a1 - t)/(2.0d0*a2)
	       else
	          r1 = ( -a1 + t)/(2.0d0*a2)
	       endif
	       r2 = a0/(a2*r1)
	       return
	    else
	       info = 0
	       return
	    endif
	endif

	end