dimension t(10),m(30),iexp(30) real*16 c,zr,zi,thta,crexp,delr,t dimension c(101),zr(50),zi(50),thta(50),crexp(50) 1 read(5,5) (t(i),i=1,10) 5 format(10a8) read(5,6) ija,ijb,ijc 6 format(3i2) if(ija.eq.9) go to 65 read (5,11) n1,nre,ncm,nred,ncmd,nexpr,nexpc,nzc,delr 11 format(8i2,1d40.30) if(n1.eq.0) go to 61 if(n1.gt.100) go to 62 n2=n1+1 read(5,12) (c(i),i=1,n2) 12 format(d50.0) read(5,13) ki,kf,mi,mf,mkii,mkfi,mkmii,mkmfi,lii,lfi,nn,nx 13 format(12i5) nrt=nre+ncm nexp=nexpr+nexpc zr(1)=0.0d0 if(nrt.eq.0) go to 16 read(5,14) (zr(i),zi(i),i=1,nrt) 14 format(2d40.30) if(nexp.eq.0) go to 16 read(5,15) (iexp(i),crexp(i),i=1,nexp) 15 format(i4,1d40.30) 16 write(6,20) (t(i),i=1,10) zr(1)=zr(1)-delr if(nzc.ne.0) go to 17 nzc=1 17 do 601 izc=1,nzc zr(1)=zr(1)+delr 20 format('1',10a8) write(6,21) (i,c(i),i=1,n2) 21 format('0 input series coefficients ci = c(i+1)', 1 /,' i c(i) '/(i4,1pd40.30)) if(nrt.eq.0) go to 33 nrd=nred+ncmd write(6,22) 22 format('0 the following singularities and their complex conjugates 1 have been specified'/,15x,'real part',30x,'imag part'/) if(nrt+ncm+nrd+ncmd.gt.36) go to 62 j=1 do 31 i=1,nrt write(6,23) zr(i),zi(i) 23 format(1x,2d40.30) if(i.gt.nrd) go to 25 write(6,24) 24 format('+',85x,'double root') 25 continue if(j.gt.nexp) go to 28 26 if(iexp(j).ne.i) go to 28 write(6,27) crexp(j) 27 format('0 a critical exponent at this point is specified ', 1'to be',1d40.30) j=j+1 go to 25 28 continue thta(i)=0.0d0 if(zi(i).eq.0.0d0) go to 31 thta(i)=2.0d0*atan(1.0d0)*abs(zi(i))/zi(i) if(zr(i).eq.0.0d0) go to 30 thta(i)=atan(zi(i)/zr(i)) if(zr(i).ge.0.0d0) go to 30 thta(i)=thta(i)+4*atan(1.0d0) 30 zr(i)=(zr(i)**2+zi(i)**2)**.5d0 31 continue 33 kip=ki+1 kfp=kf+1 li=lii+1 lf=lfi+1 mki=mkii+1 mkf=mkfi+1 mkmi=mkmii+1 mkmf=mkmfi+1 do 60 kp=kip,kfp k=kp-1 do 52 lp=li,lf l=lp-1 kx=kp 48 do 53 mmax=mi,mf do 49 i=1,kx 49 m(i)=mmax do 51 mkmst=mmax-1,mmax+1 m(2)=mkmst m(3)=mkmst if(nx.eq.0) go to 491 if((m(kp)+m(k)+m(k-1)).gt.nx) go to 50 if((m(kp)+m(k)+m(k-1)).lt.nn) go to 50 491 if(m(kp).lt.(nrt+ncm+nrd+ncmd)) go to 50 call mmmml(m,l,k,c,n1,zr,thta,nre,ncm,nred,ncmd, 1 ierr,ija,ijb,ijc,nexpr,nexpc,iexp,crexp) if(ierr.eq.1) go to 56 50 continue 56 continue 51 continue 53 continue 52 continue 60 continue 601 continue go to 1 61 continue stop 62 write(6,63) 63 format('0 excessive coefficients or specified singularities'/) 65 stop end subroutine mmmml(m,l,k,c,n1,zr,thta,nre,ncm,nred,ncmd, 1 ierr,ija,ijb,ijc,nexpr,nexpc,iexp,crexp) real*16 a,b,c,q,qk,p,zr,thta,rootr,rooti,zrc,zic 1 ,crexp,zijr,zijc dimension a(20000),b(101),c(101),q(250),qk(50),p(101) 1 ,zr(50),thta(50),rootr(50),rooti(50),crexp(50) dimension m(30),indx(250),indxr(250),iexp(50) kp=k+1 kpl=k+2 lp=l+1 m(kpl)=l ierr=0 nexp=nexpr+nexpc nrt=nre+ncm nrd=nred+ncmd nq=nrt+ncm+nrd+ncmd if(m(kp).lt.nq) go to 126 if(k.eq.0) go to 109 nq=nq+nrd+ncmd if(m(k)-nrd-ncmd) 126,109,109 109 continue mx=0 do 111 i=1,kpl if(m(i).lt.mx) go to 111 mx=m(i) 111 continue mp=mx+1 nkm=mp*kp nklmx=nkm+mp do 112 i=1,nklmx 112 indx(i)=0 in=0 if(k.lt.2) go to 119 do 118 i=2,k if(m(i).lt.0) go to 118 in=in+1 indx(i)=in indxr(in)=i 118 continue 119 continue do 121 j=1,mx do 121 i=1,kp if(j.gt.m(i)) go to 121 in=in+1 indx(j*kp+i)=in indxr(in)=j*kp+i 121 continue if(l.le.0) go to 123 do 122 i=2,lp in=in+1 indx(nkm+i)=in indxr(in)=nkm+i 122 continue if (l .le. 0) go to 123 in = in + 1 indx(1) = in indxr(in) = 1 123 continue nklm=in if(nklm.lt.1) go to 199 indmx=nklm**2 do 125 ind=1,indmx 125 a(ind)=0.0d0 nq=nq+nexpr+2*nexpc nc=nklm-nq if(nc.le.1) go to 126 if(nc.le.n1) go to 128 ierr=1 126 return 128 write (6,129) k,(m(i),i=1,kpl) 129 format('0 the k=',i4,' approximant [m0;m1;...;mk;l] =',10i4) do 136 n=1,nc in=0 do 133 j=1,mp in=in+1 if(j-1.gt.n) go to 133 if(indx(in).eq.0) go to 131 ind=(indx(in)-1)*nklm+n a(ind)=c(n-j+2) 131 do 132 i=1,k in=in+1 if(indx(in).eq.0) go to 132 ind=(indx(in)-1)*nklm+n a(ind)=c(n-j+2)*float(n-j+1)**i 132 continue 133 continue if(n.gt.l) go to 134 in=nkm+n+1 if(indx(in).eq.0) go to 134 ind=(indx(in)-1)*nklm+n a(ind)=-1.0d0 134 continue b(n)=-c(n+1) if(k.eq.0) go to 135 b(n)=-c(n+1)*float(n)**k 135 continue 136 continue mk=m(kp) mkp=mk+1 n=nc if(nrt.eq.0) go to 159 do 142 i=1,nrt n=n+1 nt=n b(n)=-1.0d0 if(thta(i).eq.0.0d0) go to 140 n=n+1 b(n)=0.0d0 140 continue do 141 j=1,mk in=(j+1)*kp ind=(indx(in)-1)*nklm+nt a(ind)=(zr(i)**j)*cos(thta(i)*j) if(thta(i).eq.0.0d0) go to 141 a(ind+1)=(zr(i)**j)*sin(thta(i)*j) 141 continue 142 continue if(nrd.eq.0) go to 159 do 158 i=1,nrd n=n+1 nt=n b(n)=0.0d0 if(thta(i).eq.0.0d0) go to 151 n=n+1 b(n)=0.0d0 151 continue do 152 j=1,mk in=(j+1)*kp ind=(indx(in)-1)*nklm+nt a(ind)=(zr(i)**j)*cos(thta(i)*j)*j if(thta(i).eq.0.0d0) go to 152 a(ind+1)=(zr(i)**j)*sin(thta(i)*j)*j 152 continue if(k.eq.0) go to 158 n=n+1 nt=n b(n)=0.0d0 if(thta(i).eq.0.0d0) go to 153 n=n+1 b(n)=0.0d0 153 continue mkm=m(k) if(indx(k).eq.0) go to 156 ind=(indx(k)-1)*nklm+nt a(ind)=1.0d0 156 do 157 j=1,mkm in=(j+1)*kp-1 ind=(indx(in)-1)*nklm+nt a(ind)=(zr(i)**j)*cos(thta(i)*j) if(thta(i).eq.0.0d0) go to 157 a(ind+1)=(zr(i)**j)*sin(thta(i)*j) 157 continue 158 continue 159 continue if(nexp.eq.0) go to 1612 do 1611 i=1,nexp n=n+1 nt=n b(n)=0.0d0 iex=iexp(i) if(thta(iex).eq.0.0d0) go to 1600 n=n+1 b(n)=0.0d0 1600 if(kp.le.2) go to 1602 if(iex.gt.nrd) go to 1601 ind2=(indx(kp-2)-1)*nklm+nt a(ind2)=1.0d0 go to 1602 1601 ind1=(indx(kp-1)-1)*nklm+nt a(ind1)=1.0d0 1602 continue do 1610 j=1,mx zijr=(zr(iex)**j)*cos(thta(iex)*j) in=(j+1)*kp indx0=indx(in)-1 if(indx0.eq.-1) go to 1603 ind0=indx0*nklm+nt a(ind0)=j*zijr*(crexp(i)-kp+2.0d0) 1603 indx1=indx(in-1)-1 if(indx1.eq.-1) go to 1604 ind1=indx1*nklm+nt a(ind1)=zijr 1604 if(iex.gt.nrd) go to 1615 if(indx0.eq.-1) go to 1605 a(ind0)=a(ind0)*(j-1)/2.0d0 if(kp.le.2) go to 1605 a(ind0)=a(ind0)*(crexp(i)-kp+3.0d0) 1605 if(indx1.eq.-1) go to 1616 a(ind1)=a(ind1)*j if(kp.le.2) go to 1615 a(ind1)=a(ind1)*(crexp(i)-kp+3.0d0) 1616 indx2=indx(in-2)-1 if(indx2.eq.-1) go to 1615 ind2=indx2*nklm+nt a(ind2)=zijr 1615 if(thta(iex).eq.0.0d0) go to 1610 zijc=(zr(iex)**j)*sin(thta(iex)*j) if(indx0.eq.-1) go to 1606 a(ind0+1)=j*zijc*(crexp(i)-kp+2.0d0) 1606 if(indx1.eq.-1) go to 1607 a(ind1+1)=zijc 1607 if(iex.gt.nrd) go to 1610 if(indx0.eq.-1) go to 1608 a(ind0+1)=a(ind0+1)*(j-1)/2.0d0 if(kp.le.2) go to 1608 a(ind0+1)=a(ind0+1)*(crexp(i)-kp+3.0d0) 1608 if(indx1.eq.-1) go to 1609 a(ind1+1)=a(ind1+1)*j if(kp.le.2) go to 1610 a(ind1+1)=a(ind1+1)*(crexp(i)-kp+3.0d0) 1609 if(indx2.eq.-1) go to 1610 a(ind2+1)=zijc 1610 continue 1611 continue 1612 call simq(a,b,nklm,ks) if(ks.ne.1) go to 162 write(6,161) 161 format('0 singular matrix') return 162 continue do 170 i=1,nklmx 170 q(i)=0.0d0 q(kp)=1.0d0 do 171 i=1,nklm j=indxr(i) 171 q(j)=b(i) q(nkm+1)=q(1)*c(1) do 172 i=1,mp 172 p(i)=q(nkm+i) if(ija.eq.0) go to 187 write (6,181) mx,k 181 format('0 coefficients of recurrence relation of',i4, 1 ' terms and order',i4,/ 2 ,12x,'d0',20x,'q1',20x,'...',20x,'qk',20x,'p'/) do 186 i=1,mp ns=(i-1)*kp+1 nf=i*kp if(k.gt.2) go to 184 write(6,183) i,(q(j),j=ns,nf),p(i) 183 format(i4,1p4d28.20) go to 186 184 write(6,185) i,(q(j),j=ns,nf),p(i) 185 format(i2,(1p10d13.6)) 186 continue 187 if(ijb.eq.0) go to 189 call gen(q,p,c,mx,k,nc) 189 if(ijc.eq.0) go to 195 nnf=mp*kp call rinf(q,k,nnf) 195 do 190 j=1,mkp 190 qk(j)=q(kp*j) call root(qk,mk,rootr,rooti,ier) if(ier.eq.3) go to 191 if(ier.gt.0) go to 199 191 do 198 i=1,mk zrc=rootr(i) zic=rooti(i) if(abs(rooti(i)).gt.1.0d-08) go to 197 call rind(zrc,zic,q,mp,kp,i) go to 198 197 call cind(zrc,zic,q,mp,kp,i) 198 continue 199 return end subroutine polrt(xcof,cof,m,rootr,rooti,ier) real*16 xo,yo,x,y,xpr,ypr,ux,uy,v,yt,xt,u,xt2,yt2,sumsq, 1 dx,dy,temp,alpha,fi,xcof,cof,rootr,rooti dimension xcof(50),cof(50),rootr(50),rooti(50) ifit=0 n=m ier=0 if(xcof(n+1))10,25,10 10 if(n) 15,15,32 15 ier=1 20 return 25 ier=4 go to 20 30 ier=2 go to 20 32 if(n-50) 35,35,30 35 nx=n nxx=n+1 n2=1 kj1 = n+1 do 40 l=1,kj1 mt=kj1-l+1 40 cof(mt)=xcof(l) 45 xo=.00500101 yo=0.01000101 in=0 50 x=xo xo=-10.0*yo yo=-10.0*x x=xo y=yo in=in+1 go to 59 55 ifit=1 xpr=x ypr=y 59 ict=0 60 ux=0.0 uy=0.0 v =0.0 yt=0.0 xt=1.0 u=cof(n+1) if(u) 65,130,65 65 do 70 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 xt=xt2 uy=uy-fi*yt*temp 70 yt=yt2 sumsq=ux*ux+uy*uy if(sumsq) 75,110,75 75 dx=(v*uy-u*ux)/sumsq x=x+dx dy=-(u*uy+v*ux)/sumsq y=y+dy 78 if(abs(dy)+abs(dx)-1.0d-12) 100,80,80 80 ict=ict+1 if(ict-500) 60,85,85 85 if(ifit)100,90,100 90 if(in-5) 50,95,95 95 ier=3 go to 20 100 do 105 l=1,nxx mt=kj1-l+1 temp=xcof(mt) xcof(mt)=cof(l) 105 cof(l)=temp itemp=n n=nx nx=itemp if(ifit) 120,55,120 110 if(ifit) 115,50,115 115 x=xpr y=ypr 120 ifit=0 122 if(abs(y)-1.0d-10*abs(x)) 135,125,125 125 alpha=x+x sumsq=x*x+y*y n=n-2 go to 140 130 x=0.0 nx=nx-1 nxx=nxx-1 135 y=0.0 sumsq=0.0 alpha=x n=n-1 140 cof(2)=cof(2)+alpha*cof(1) if(n.lt.2) go to 155 145 do 150 l=2,n 150 cof(l+1)=cof(l+1)+alpha*cof(l)-sumsq*cof(l-1) 155 rooti(n2)=y rootr(n2)=x n2=n2+1 if(sumsq) 160,165,160 160 y=-y sumsq=0.0 go to 155 165 if(n) 20,20,45 end subroutine gen(q,p,c,m,k,n) real*16 q,p,c,a,s,d,g dimension q(250),p(101),c(101),g(121) kp=k+1 lx=n+8 g(1)=c(1) do 20 l=1,lx s=0.0d0 do 10 j=1,m if(j.gt.l) go to 12 a=q(j*kp+1) do 9 i=1,k 9 a=a+q(j*kp+i+1)*float(l-j)**i s=s+a*g(l-j+1) 10 continue 12 d=q(1) do 13 i=1,k 13 d=d+q(i+1)*float(l)**i if(abs(d).gt.1.0d-12) go to 30 g(l+1)=c(l+1) go to 20 30 g(l+1)=-s/d if(l.gt.m) go to 20 g(l+1)=(p(l+1)-s)/d 20 continue 21 write(6,22) (i,g(i),i=1,lx) 22 format('0 coefficients generated by recurrence',/(i4,1p1d40.30)) return end subroutine simq(a,b,n,ks) real*16 a,b,tol,biga,save dimension a(11000),b(101) tol=0.0d0 ks=0 jj=-n do 65 j=1,n jy=j+1 jj=jj+n+1 biga=0 it=jj-j do 30 i=j,n ij=it+i if(abs(biga)-abs(a(ij))) 20,30,30 20 biga=a(ij) imax=i 30 continue if(abs(biga)-tol) 35,35,40 35 ks=1 return 40 i1=j+n*(j-2) it=imax-j do 50 k=j,n i1=i1+n i2=i1+it save=a(i1) a(i1)=a(i2) a(i2)=save 50 a(i1)=a(i1)/biga save=b(imax) b(imax)=b(j) b(j)=save/biga if(j-n) 55,70,55 55 iqs=n*(j-1) do 65 ix=jy,n ixj=iqs+ix it=j-ix do 60 jx=jy,n ixjx=n*(jx-1)+ix jjx=ixjx+it 60 a(ixjx)=a(ixjx)-(a(ixj)*a(jjx)) 65 b(ix)=b(ix)-(b(j)*a(ixj)) 70 ny=n-1 it=n*n do 80 j=1,ny ia=it-j ib=n-j ic=n do 80 k=1,j b(ib)=b(ib)-a(ia)*b(ic) ia=ia-n 80 ic=ic-1 return end subroutine rind(zr,zi,q,mp,kp,i) real*16 zr,zi,q,rindx,s,s1,s2,d dimension q(250) complex*16 rind1,rind2,z if(kp.le.1) go to 52 s=q(kp-1) d=0.0d0 do 10 j=2,mp l=j-1 s=s+q(j*kp-1)*zr**l d=d+q(j*kp)*l*zr**l 10 continue if(abs(d).lt.1.0d-10) go to 20 rindx=kp-2.0d0-s/d write(6,11) i,zr,zi,rindx 11 format(i3,2d20.12,d19.11) if(abs(d).lt.1.0d-02) go to 12 return 12 write(6,13) 13 format('+',82x,' suspected double root') 20 s2=0.0d0 s1=0.0d0 d=0.0d0 do 21 j=2,mp l=j-1 s2=s2+q(j*kp-2)*zr**l s1=s1+q(j*kp-1)*l*zr**l d=d+q(j*kp)*l*(l-1)*zr**l 21 continue if(abs(d).lt.1.0d-08) go to 52 s1=.5d0-s1/d if(kp.le.2) go to 43 s2=2.0d0*(s2+q(kp-2))/d z=s1**2-s2 rind1=kp-3.0d0+s1+sqrt(z) rind2=rind1-2.0d0*sqrt(z) write(6,41) i,zr,zi,rind1,rind2 41 format(i3,2d20.12,4d19.11) return 43 rindx=kp-3.0d0+2.0d0*s1 write(6,11) i,zr,zi,rindx return 52 write(6,53) i,zr,zi 53 format(i3,2d20.12,' exponents are infinite or indeterminate') return end subroutine cind(zr,zi,q,mp,kp,i) real*16 zr,zi,q, ci, cr dimension q(250) complex*16 rind1,rind2,cindx,z,zc,s,s1,s2,d if(kp.le.1) go to 52 cr=zr ci=zi zc=dcmplx(cr,ci) s=q(kp-1) d=dcmplx(0.d0,0.d0) do 10 j=2,mp l=j-1 s=s+q(j*kp-1)*zc**l d=d+q(j*kp)*l*zc**l 10 continue if(abs(d).lt.1.d-08) go to 20 cindx=kp-2.0d0-s/d write(6,11) i,zr,zi,cindx 11 format(i3,2d20.12,2d19.11) if(abs(d).lt.1.d-02) go to 12 return 12 write(6,13) 13 format('+',82x,' suspected double root') 20 s2=dcmplx(0.d0,0.d0) s1=dcmplx(0.d0,0.d0) d=dcmplx(0.d0,0.d0) do 21 j=2,mp l=j-1 s2=s2+q(j*kp-2)*zc**l s1=s1+q(j*kp-1)*l*zc**l d=d+q(j*kp)*l*(l-1)*zc**l 21 continue if(abs(d).lt.1.d-08) go to 52 s1=-s1/d+dcmplx(.5d0,0.d0) if(kp.le.2) go to 43 s2=2.0d0*(s2+q(kp-2))/d z=s1**2-s2 rind1=kp-3.0d0+s1+sqrt(z) rind2=rind1-2.0d0*sqrt(z) write(6,41) i,zr,zi,rind1,rind2 41 format(i3,2d20.12,4d19.11) return 43 cindx=kp-3.0d0+2.0d0*s1 write(6,11) i,zr,zi,cindx return 52 write(6,53) i,zr,zi 53 format(i3,2d20.12,' exponents are infinite or indeterminate') return end subroutine root(q,m,rootr,rooti,ier) REAL*16 q,rootr,rooti,w dimension q(50),rootr(50),rooti(50),w(50) call polrt(q,w,m,rootr,rooti,ier) if(ier-1)90,60,70 60 write(6,65) 65 format(//' order of polynomial less than one') go to 100 70 if(ier-3) 75,80,78 75 write(6,77) 77 format(//' order of polynomial greater than 36') go to 100 78 write(6,79) 79 format(//' higher order coefficient is zero') go to 100 80 write(6,85) 85 format(//' unable to determine root, those already found are') 90 write(6,95) 95 format('0',10x,'singularity position',21x,'critical exponent', 1 20x,'2nd critical exponent'/,10x,'real part',10x,'imag part', 2 10x,'real part',10x,'imag part',15x,'(double root case)') 100 return end subroutine rinf(q,n,nf) implicit REAL*16 (a-h,o-z) dimension q(250),rootr(50),rooti(50),r(50),s(50) nn=n+1 nfn=nf-n do 100 i = 1,nn r(i)=q(nfn+i-1) 100 continue call polrt(r,s,n,rootr,rooti,ier) write (6,150) (i,rootr(i),rooti(i),i=1,n) 150 format ('0 exponents at infinity ',/(4(i3,2d14.7))) return end