41 parameter(expmax=87.49823)
44 parameter(infnty=1.0e38)
48 parameter(minlog=1.0e-37)
54 REAL a,alpha,b,beta,delta,gamma,k1,k2,olda,oldb,r,s,t,u1,u2,v,w,y,
67 SAVE olda,oldb,alpha,beta,gamma,k1,k2,a,b
71 DATA olda,oldb/-1.0e37,-1.0e37/
74 qsame = (olda.EQ.aa) .AND. (oldb.EQ.bb)
77 IF (.NOT. (aa.LT.minlog.OR.bb.LT.minlog)) go
to 10
78 WRITE (*,*)
' AA or BB < ',minlog,
' in GENBET - Abort!'
79 WRITE (*,*)
' AA: ',aa,
' BB ',bb
80 CALL xstopx(
' AA or BB too small in GENBET - Abort!')
84 20
IF (.NOT. (
min(aa,bb).GT.1.0)) go
to 100
96 beta =
sqrt((alpha-2.0)/ (2.0*a*b-alpha))
104 v = beta*
log(u1/ (1.0-u1))
106 IF (v.GT.expmax) go
to 55
110 IF (w.GT.infnty/a) go
to 55
116 r = gamma*v - 1.3862944
121 IF ((s+2.609438).GE. (5.0*z)) go
to 70
135 IF (alpha/(b+w).LT.minlog) go
to 40
137 IF ((r+alpha*
log(alpha/ (b+w))).LT.t) go
to 40
141 70
IF (.NOT. (aa.EQ.a)) go
to 80
154 100
IF (qsame) go
to 110
160 k1 = delta* (0.0138889+0.0416667*b)/ (a*beta-0.777778)
161 k2 = 0.25 + (0.5+0.25/delta)*b
168 IF (u1.GE.0.5) go
to 130
174 IF ((0.25*u2+z-y).GE.k1) go
to 120
180 IF (.NOT. (z.LE.0.25)) go
to 160
181 v = beta*
log(u1/ (1.0-u1))
186 IF (a.GT.1.0) go
to 135
188 IF (v.GT.expmax) go
to 132
192 IF (w.GT.expmax) go
to 140
197 135
IF (v.GT.expmax) go
to 140
199 IF (w.GT.infnty/a) go
to 140
205 160
IF (z.GE.k2) go
to 120
212 170 v = beta*
log(u1/ (1.0-u1))
215 IF (a.GT.1.0) go
to 175
217 IF (v.GT.expmax) go
to 172
221 IF (w.GT.expmax) go
to 180
226 175
IF (v.GT.expmax) go
to 180
228 IF (w.GT.infnty/a) go
to 180
236 190
IF (alpha/(b+w).LT.minlog) go
to 120
237 IF ((alpha* (
log(alpha/ (b+w))+v)-1.3862944).LT.
log(z)) go
to 120
241 200
IF (.NOT. (a.EQ.aa)) go
to 210
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)
real function genbet(aa, bb)
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
charNDArray max(char d, const charNDArray &m)
octave_value sqrt(void) const
charNDArray min(char d, const charNDArray &m)