1 SUBROUTINE zmlri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)
11 DOUBLE PRECISION ACK, AK, AP, AT, AZ, BK, CKI, CKR, CNORMI,
12 * cnormr, conei, coner, fkap, fkk, flam, fnf, fnu, pti, ptr, p1i,
13 * p1r, p2i, p2r, raz, rho, rho2, rzi, rzr, scle, sti,
str, sumi,
14 * sumr, tfnf, tol, tst, yi, yr, zeroi, zeror, zi, zr,
dgamln,
16 INTEGER I, IAZ, IDUM, IFNU, INU, ITIME, K, KK, KM, KODE, M, N, NZ
18 DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
25 at = dble(float(iaz)) + 1.0d0
38 rho = ack + dsqrt(ack*ack-1.0d0)
40 tst = (rho2+rho2)/((rho2-1.0d0)*(rho-1.0d0))
49 p2r = p1r - (ckr*ptr-cki*pti)
50 p2i = p1i - (cki*ptr+ckr*pti)
56 IF (ap.GT.tst*ak*ak) go
to 20
63 IF (inu.LT.iaz) go
to 40
71 at = dble(float(inu)) + 1.0d0
82 p2r = p1r - (ckr*ptr-cki*pti)
83 p2i = p1i - (ckr*pti+cki*ptr)
89 IF (ap.LT.tst) go
to 30
90 IF (itime.EQ.2) go
to 40
92 flam = ack + dsqrt(ack*ack-1.0d0)
93 fkap = ap/xzabs(p1r,p1i)
94 rho = dmin1(flam,fkap)
95 tst = tst*dsqrt(rho/(rho*rho-1.0d0))
104 kk = max0(i+iaz,k+inu)
105 fkk = dble(float(kk))
113 fnf = fnu - dble(float(ifnu))
115 bk =
dgamln(fkk+tfnf+1.0d0,idum) -
dgamln(fkk+1.0d0,idum) -
124 p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
125 p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
128 ak = 1.0d0 - tfnf/(fkk+tfnf)
130 sumr = sumr + (ack+bk)*p1r
131 sumi = sumi + (ack+bk)*p1i
141 p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
142 p2i = p1i + (fkk+fnf)*(rzi*ptr+rzr*pti)
145 ak = 1.0d0 - tfnf/(fkk+tfnf)
147 sumr = sumr + (ack+bk)*p1r
148 sumi = sumi + (ack+bk)*p1i
156 IF (ifnu.LE.0) go
to 90
160 p2r = p1r + (fkk+fnf)*(rzr*ptr-rzi*pti)
161 p2i = p1i + (fkk+fnf)*(rzr*pti+rzi*ptr)
164 ak = 1.0d0 - tfnf/(fkk+tfnf)
166 sumr = sumr + (ack+bk)*p1r
167 sumi = sumi + (ack+bk)*p1i
174 IF (kode.EQ.2) ptr = zeror
175 CALL
xzlog(rzr, rzi,
str, sti, idum)
178 ap =
dgamln(1.0d0+fnf,idum)
194 CALL
zmlt(ckr, cki, ptr, pti, cnormr, cnormi)
196 str = yr(i)*cnormr - yi(i)*cnormi
197 yi(i) = yr(i)*cnormi + yi(i)*cnormr
subroutine xzlog(AR, AI, BR, BI, IERR)
subroutine xzexp(AR, AI, BR, BI)
double precision function d1mach(i)
double precision function dgamln(Z, IERR)
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
subroutine zmlt(AR, AI, BR, BI, CR, CI)
subroutine zmlri(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL)