      subroutine factor(b,k4,rootr,rooti,kinsid,kprint,eps)
c       sets up problem, calls dproot, 
c       and checks residual values at roots  
      implicit double precision (a-h,o-z)
      double complex z,res,jay
      double precision b(1),rootr(1),rooti(1),coe(37)
      jay=(0.d0,1.d0)
      pi=4.d0*datan2(1.d0,1.d0)
      do 550 i=1,k4
550   coe(i)=b(i)
      k4m=k4-1
      call dproot(k4m,coe,rootr,rooti,kerr,kprint,eps)
c      write(0,600)kerr
600   format(' return from dproot with kerr=',i5)
      if(kerr.gt.0)stop
      kinsid=0
      resmax=0.d0
      rmax=0.d0
      rmin=2.d0**(32)
      dist=2.d0**(32)
      amax=2.d0**(32)
      r2=amax**(1.d0/k4)
      do 701 j=1,k4m
      z=rootr(j)+jay*rooti(j)
      r=dsqrt(rootr(j)**2+rooti(j)**2)
c         skip residue calculation if root is too big
      if(r.gt.r2)goto713
      res=b(k4)
      do 705 k=2,k4
705   res=res*z+b(k4-k+1)
      partr=res
      parti=-jay*res
      resmag=dsqrt(partr**2+parti**2)
      if(resmax.le.resmag)resmax=resmag
713   if(rmax.lt.r)rmax=r
      if(rmin.gt.r)rmin=r
      if(r.lt.1.d0)kinsid=kinsid+1
      distr=dabs(r-1.d0)
      if(dist.gt.distr)dist=distr
701   continue
c     write(0,703)resmax
c     write(0,704)rmax,rmin,dist
c703   format('resmax= ',d20.10)
c704   format('rmax= ',d20.10/'rmin= ',d20.10/'dist=',d20.10)
      return
      end

      subroutine dproot(mm,a,rootr,rooti,kerr,kprint,eps)
c        mm=degree of polynomial
c        a=coefficient array, lowest to highest degree
c        kprint=1 for full printing
c        kerr=0 is normal return
      implicit double precision (a-h,o-z)
      double complex b(37),c(37),p,pp,z,w
      double complex bb(37),cc(37),jay
      double precision a(1),rootr(1),rooti(1)
      double precision save(37)
      jay=(0.d0,1.d0)
      mmp=mm+1
      m=mm
      mp=mmp
      do 700 i=1,mp
700   save(i)=a(i)
c       kount is number of iterations so far
      kount=0
c       kmax is maximum total number of iterations allowed
      kmax =20*m
c       newst is number of re-starts
      newst=0
c       ktrym is number of attempted iterations before re-starting
      ktrym=20
c       kpolm is number of attempted iterations before polishing is stopped
      kpolm=20
c       amax is the largest number we allow
      amax=2.d0**(32)
      amin=1.d0/amax
c       rr1 and rr2 are radii within which we work for polishing
      rr1=amin**(1.d0/m)
      rr2=amax**(1.d0/m)
c     eps is the tolerance for convergence
      sqteps=dsqrt(eps)
c        main loop; m is current degree
10    if(m.le.0)goto200
c        new z, a point on the unit circle
      rkount=kount
      z=dcos(rkount)+jay*dsin(rkount)
      ktry=0
c       r1 and r2 are boundaries of an expanding annulus within which we work
      r1=amin**(1.d0/m)
      r2=amax**(1.d0/m)
c        inside loop
20    partr=z
      parti=-jay*z
      size=dsqrt(partr**2+parti**2)
      if(size.lt.r1 .or. size.gt.r2)goto300
      if(ktry.ge.ktrym)goto300
      ktry=ktry+1
      if(kount.ge.kmax)goto400
      kount=kount+1
c        get value of polynomial at z, synthetic division
      b(mp)=a(mp)
      do 30 j=1,m
      k=m-j+1
30    b(k)=z*b(k+1)+a(k)
      p=b(1)
      partr=p
      parti=-jay*p
      if(dsqrt(partr**2+parti**2).gt.amax)goto300
c        get value of derivative at z, synthetic division
      c(mp)=b(mp)
      mdec=m-1
      do 60 j=1,mdec
      k=m-j+1
60    c(k)=z*c(k+1)+b(k)
      pp=c(2)
      partr=pp
      parti=-jay*pp
      if(dsqrt(partr**2+parti**2).lt.amin)goto300
c        test for convergence
      partr=p
      parti=-jay*p
      size=dsqrt(partr**2+parti**2)
      if(size.gt.eps)goto775
      nroot=mm-m+1
      if(kprint.eq.1)write(0,776)kount,nroot
776   format('kount=',i5,'  root no.=',i5)
      goto500
775   continue
      z=z-p/pp
      goto20
c        end of main loop

c        normal return
200   kerr=0
      goto600

c        new start
300   rkount=kount
      z=dcos(rkount)+jay*dsin(rkount)
      ktry=0
      newst=newst+1
      goto20

c        too many iterations
400   kerr=400
      goto600

c        root z located
c        polish z to get w
500   w=z
      kpol=0
510   partr=w
      parti=-jay*w
      size=dsqrt(partr**2+parti**2)
c       give up polishing if w is outside annulus
      if(size.lt.rr1 .or. size.gt.rr2)goto501
c       give up polishing if kpol>=kpolm
      if(kpol.ge.kpolm)goto501
      kpol=kpol+1
      if(kount.ge.kmax)goto400
      kount=kount+1
      bb(mmp)=save(mmp)
      do 530 j=1,mm
      k=mm-j+1
530   bb(k)=w*bb(k+1)+save(k)
      p=bb(1)
      partr=p
      parti=-jay*p
      if(dsqrt(partr**2+parti**2).gt.amax)goto300
      cc(mmp)=bb(mmp)
      mdec=mm-1
      do 560 j=1,mdec
      k=mm-j+1
560   cc(k)=w*cc(k+1)+bb(k)
      pp=cc(2)
      partr=pp
      parti=-jay*pp
      if(dsqrt(partr**2+parti**2).lt.amin)goto300
      partr=p
      parti=-jay*p
      size=dsqrt(partr**2+parti**2)
c       test for convergence of polishing
      if(size.le.eps)goto501
      w=w-p/pp
      goto510

c        deflate
501   b(mp)=a(mp)
      do 830 j=1,m
      k=m-j+1
830   b(k)=z*b(k+1)+a(k)
      p=b(1)
      rootr(m)=w
      rooti(m)=-jay*w
      m=m-1
      mp=mp-1
      parti=-jay*w
      if(dabs(parti).gt.sqteps)goto140
c        real root
      rooti(m+1)=0.d0
      do 100 j=1,mp
100   a(j)=b(j+1)
      goto10
c        complex root
140   partr=z
      parti=-jay*z
      z=partr-jay*parti
      c(mp)=b(mp+1)
      do 110 j=1,m
      k=m-j+1
110   c(k)=z*c(k+1)+b(k+1)
      rootr(m)=w
      rooti(m)=-(-jay*w)
      m=m-1
      mp=mp-1
      do 130 j=1,mp
130   a(j)=c(j+1)
      goto10
c        report and return
600   real1=kount
      real2=mm
      temp=real1/real2
c     write(0,150)kount,temp
150   format(' kount=',i10,' kount/root=',f15.5)
c     write(0,151)newst,kerr
151   format(' new starts=',i10,' kerr=',i10)
      return
      end
