2 SUBROUTINE psifn (X, N, KODE, M, ANS, NZ, IERR)
107 INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ
109 REAL ANS, ARG, B, DEN, ELIM, EPS, FLN, FN, FNP, FNS, FX, RLN,
110 * rxsq, r1m4, r1m5, s, slope, t, ta, tk, tol, tols, trm, trmr,
111 * tss, tst, tt, t1, t2, wdtol, x, xdmln, xdmy, xinc, xln, xm,
114 dimension b(22), trm(22), trmr(100), ans(*)
120 DATA b(1), b(2), b(3), b(4), b(5), b(6), b(7), b(8), b(9), b(10),
121 * b(11), b(12), b(13), b(14), b(15), b(16), b(17), b(18), b(19),
122 * b(20), b(21), b(22) /1.00000000000000000e+00,
123 * -5.00000000000000000e-01,1.66666666666666667e-01,
124 * -3.33333333333333333e-02,2.38095238095238095e-02,
125 * -3.33333333333333333e-02,7.57575757575757576e-02,
126 * -2.53113553113553114e-01,1.16666666666666667e+00,
127 * -7.09215686274509804e+00,5.49711779448621554e+01,
128 * -5.29124242424242424e+02,6.19212318840579710e+03,
129 * -8.65802531135531136e+04,1.42551716666666667e+06,
130 * -2.72982310678160920e+07,6.01580873900642368e+08,
131 * -1.51163157670921569e+10,4.29614643061166667e+11,
132 * -1.37116552050883328e+13,4.88332318973593167e+14,
133 * -1.92965793419400681e+16/
138 IF (x.LE.0.0e0) ierr=1
140 IF (kode.LT.1 .OR. kode.GT.2) ierr=1
142 IF (ierr.NE.0)
RETURN
144 nx =
min(-i1mach(12),i1mach(13))
146 r1m4 = r1mach(4)*0.5e0
147 wdtol =
max(r1m4,0.5e-18)
151 elim = 2.302e0*(nx*r1m5-3.0e0)
161 IF (
abs(t).GT.elim) go
to 290
162 IF (x.LT.wdtol) go
to 260
166 rln = r1m5*i1mach(11)
167 rln =
min(rln,18.06e0)
168 fln =
max(rln,3.0e0) - 3.0e0
169 yint = 3.50e0 + 0.40e0*fln
170 slope = 0.21e0 + fln*(0.0006038e0*fln+0.008677e0)
175 xm = -2.302e0*rln -
min(0.0e0,xln)
181 IF (
abs(arg).LT.1.0e-3) xm = -arg
184 IF (xm.GT.7.0e0 .AND. fln.LT.15.0e0) go
to 200
189 IF (x.GE.xmin) go
to 60
202 IF (tk.GT.elim) go
to 380
207 IF (nn.NE.0) t1 = tt + 1.0e0/fn
208 rxsq = 1.0e0/(xdmy*xdmy)
212 IF (
abs(s).LT.tst) go
to 80
215 t = t*((tk+fn+1.0e0)/(tk+1.0e0))*((tk+fn)/(tk+2.0e0))*rxsq
217 IF (
abs(trm(k)).LT.tst) go
to 80
223 IF (xinc.EQ.0.0e0) go
to 100
229 IF (nx.GT.nmax) go
to 390
230 IF (nn.EQ.0) go
to 160
244 IF (fn.EQ.0.0e0) go
to 180
254 IF (fn.NE.0.0e0) t1 = tt + 1.0e0/fn
257 IF (
abs(s).LT.tst) go
to 120
260 trm(k) = trm(k)*fnp/tk
261 IF (
abs(trm(k)).LT.tst) go
to 120
267 IF (xinc.EQ.0.0e0) go
to 140
268 IF (fn.EQ.0.0e0) go
to 160
280 IF (fn.EQ.0.0e0) go
to 180
288 s = s + 1.0e0/(x+nx-i)
291 IF (kode.EQ.2) go
to 190
295 IF (xdmy.EQ.x)
RETURN
315 IF (n.NE.0) go
to 220
316 IF (kode.EQ.2) ans(1) = s + xln
332 IF (trm(i).LT.tols) go
to 240
343 IF (mm.EQ.1) go
to 280
351 IF (kode.EQ.2) ans(1) = ans(1) + xln
354 IF (t.GT.0.0e0) go
to 380
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
subroutine psifn(X, N, KODE, M, ANS, NZ, 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)
charNDArray max(char d, const charNDArray &m)
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
charNDArray min(char d, const charNDArray &m)