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)
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
charNDArray max(char d, const charNDArray &m)
charNDArray min(char d, const charNDArray &m)