c        subroutine dpolrt
c           computes the real and complex roots of a real polynomial
      subroutine dpolrt (xcof,cof,m,rootr,rooti,ier)
      dimension xcof(1), cof(1), rootr(1), rooti(1)
      double precision xo, yo, x, y, xpr, ypr, ux, uy, v, yt, xt, u, xt2
     1 , yt2, sumsq, dx, dy, temp, alpha
      double precision xcof, cof, rootr, rooti
      ifit=0
      n=m
      ier=0
      if (xcof(n+1)) 10,40,10
   10 if (n) 20,20,60
c
c        set error code to 1
c
   20 ier=1
   30 return
c
c        set error code to 4
c
   40 ier=4
      go to 30
c
c        set error code to 2
c
   50 ier=2
      go to 30
   60 if (n-36) 70,70,50
   70 nx=n
      nxx=n+1
      n2=1
      kj1=n+1
      do 80 l=1,kj1
      mt=kj1-l+1
   80 cof(mt)=xcof(l)
c
c        set initial values
c
   90 xo=.00500101d0
      yo=0.01000101d0
c
c        zero initial value counter
c
      in=0
  100 x=xo
c
c        increment initial values and counter
c
      xo=-10.d0*yo
      yo=-10.d0*x
c
c        set x and y to current value
c
      x=xo
      y=yo
      in=in+1
      go to 120
  110 ifit=1
      xpr=x
      ypr=y
c
c        evaluate polynomial and derivatives
c
  120 ict=0
  130 ux=0.0d0
      uy=0.0d0
      v=0.0d0
      yt=0.0d0
      xt=1.0d0
      u=cof(n+1)
      if (u) 140,270,140
  140 do 150 i=1,n
      l=n-i+1
      temp=cof(l)
      xt2=x*xt-y*yt
      yt2=x*yt+y*xt
      u=u+temp*xt2
      v=v+temp*yt2
      fi=i
      ux=ux+fi*xt*temp
      uy=uy-fi*yt*temp
      xt=xt2
  150 yt=yt2
      sumsq=ux*ux+uy*uy
      if (sumsq) 160,230,160
  160 dx=(v*uy-u*ux)/sumsq
      x=x+dx
      dy=-(u*uy+v*ux)/sumsq
      y=y+dy
      if (dabs(dy)+dabs(dx)-1.0d-12) 210,170,170
c
c        step iteration counter
c
  170 ict=ict+1
      if (ict-500) 130,180,180
  180 if (ifit) 210,190,210
  190 if (in-5) 100,200,200
c
c        set error code to 3
c
  200 ier=3
      go to 30
  210 do 220 l=1,nxx
      mt=kj1-l+1
      temp=xcof(mt)
      xcof(mt)=cof(l)
  220 cof(l)=temp
      itemp=n
      n=nx
      nx=itemp
      if (ifit) 250,110,250
  230 if (ifit) 240,100,240
  240 x=xpr
      y=ypr
  250 ifit=0
      if (dabs(y)-1.d-10*dabs(x)) 280,260,260
  260 alpha=x+x
      sumsq=x*x+y*y
      n=n-2
      go to 290
  270 x=0.0d0
      nx=nx-1
      nxx=nxx-1
  280 y=0.0d0
      sumsq=0.0d0
      alpha=x
      n=n-1
  290 cof(2)=cof(2)+alpha*cof(1)
      do 300 l=2,n
  300 cof(l+1)=cof(l+1)+alpha*cof(l)-sumsq*cof(l-1)
  310 rooti(n2)=y
      rootr(n2)=x
      n2=n2+1
      if (sumsq) 320,330,320
  320 y=-y
      sumsq=0.0d0
      go to 310
  330 if (n) 30,30,90
      end