Télécharger calerf.eso

Retour à la liste

Numérotation des lignes :

  1. C CALERF SOURCE CHAT 05/01/12 21:45:49 5004
  2. SUBROUTINE CALERF(ARG,RESULT,JINT)
  3. C------------------------------------------------------------------
  4. C
  5. C This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x)
  6. C for a real argument x. It contains three FUNCTION type
  7. C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
  8. C and one SUBROUTINE type subprogram, CALERF. The calling
  9. C statements for the primary entries are:
  10. C
  11. C Y=ERF(X) (or Y=DERF(X)),
  12. C
  13. C Y=ERFC(X) (or Y=DERFC(X)),
  14. C and
  15. C Y=ERFCX(X) (or Y=DERFCX(X)).
  16. C
  17. C The routine CALERF is intended for internal packet use only,
  18. C all computations within the packet being concentrated in this
  19. C routine. The function subprograms invoke CALERF with the
  20. C statement
  21. C
  22. C CALL CALERF(ARG,RESULT,JINT)
  23. C
  24. C where the parameter usage is as follows
  25. C
  26. C Function Parameters for CALERF
  27. C call ARG Result JINT
  28. C
  29. C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
  30. C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1
  31. C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2
  32. C
  33. C The main computation evaluates near-minimax approximations
  34. C from "Rational Chebyshev approximations for the error function"
  35. C by W. J. Cody, Math. Comp., 1969, PP. 631-638. This
  36. C transportable program uses rational functions that theoretically
  37. C approximate erf(x) and erfc(x) to at least 18 significant
  38. C decimal digits. The accuracy achieved depends on the arithmetic
  39. C system, the compiler, the intrinsic functions, and proper
  40. C selection of the machine-dependent constants.
  41. C
  42. C*******************************************************************
  43. C*******************************************************************
  44. C
  45. C Explanation of machine-dependent constants
  46. C
  47. C XMIN = the smallest positive floating-point number.
  48. C XINF = the largest positive finite floating-point number.
  49. C XNEG = the largest negative argument acceptable to ERFCX;
  50. C the negative of the solution to the equation
  51. C 2*exp(x*x) = XINF.
  52. C XSMALL = argument below which erf(x) may be represented by
  53. C 2*x/sqrt(pi) and above which x*x will not underflow.
  54. C A conservative value is the largest machine number X
  55. C such that 1.0 + X = 1.0 to machine precision.
  56. C XBIG = largest argument acceptable to ERFC; solution to
  57. C the equation: W(x) * (1-0.5/x**2) = XMIN, where
  58. C W(x) = exp(-x*x)/[x*sqrt(pi)].
  59. C XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
  60. C machine precision. A conservative value is
  61. C 1/[2*sqrt(XSMALL)]
  62. C XMAX = largest acceptable argument to ERFCX; the minimum
  63. C of XINF and 1/[sqrt(pi)*XMIN].
  64. C
  65. C Approximate values for some important machines are:
  66. C
  67. C XMIN XINF XNEG XSMALL
  68. C
  69. C CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15
  70. C CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15
  71. C IEEE (IBM/XT,
  72. C SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8
  73. C IEEE (IBM/XT,
  74. C SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16
  75. C IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17
  76. C UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18
  77. C VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17
  78. C VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16
  79. C
  80. C
  81. C XBIG XHUGE XMAX
  82. C
  83. C CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293
  84. C CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465
  85. C IEEE (IBM/XT,
  86. C SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37
  87. C IEEE (IBM/XT,
  88. C SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307
  89. C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75
  90. C UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307
  91. C VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38
  92. C VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307
  93. C
  94. C*******************************************************************
  95. C*******************************************************************
  96. C
  97. C Error returns
  98. C
  99. C The program returns ERFC = 0 for ARG .GE. XBIG;
  100. C
  101. C ERFCX = XINF for ARG .LT. XNEG;
  102. C and
  103. C ERFCX = 0 for ARG .GE. XMAX.
  104. C
  105. C
  106. C Intrinsic functions required are:
  107. C
  108. C ABS, AINT, EXP
  109. C
  110. C
  111. C Author: W. J. Cody
  112. C Mathematics and Computer Science Division
  113. C Argonne National Laboratory
  114. C Argonne, IL 60439
  115. C
  116. C Latest modification: March 19, 1990
  117. C
  118. C------------------------------------------------------------------
  119. IMPLICIT INTEGER(I-N)
  120. IMPLICIT REAL*8(A-H,O-Z)
  121. -INC CCREEL
  122. INTEGER I,JINT
  123. REAL*8
  124. 1 A,ARG,B,C,D,DEL,FOUR,HALF,P,ONE,Q,RESULT,SIXTEN,SQRPI,
  125. 2 TWO,THRESH,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL,
  126. 3 Y,YSQ,ZERO
  127. DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5)
  128. C------------------------------------------------------------------
  129. C Mathematical constants
  130. C------------------------------------------------------------------
  131. DATA FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/,
  132. 1 SQRPI/5.6418958354775628695D-1/,THRESH/0.46875D0/,
  133. 2 SIXTEN/16.0D0/
  134. C------------------------------------------------------------------
  135. C Machine-dependent constants
  136. C------------------------------------------------------------------
  137. C DATA XINF,XNEG,XSMALL/1.79D308,-26.628D0,1.11D-16/,
  138. C 1 XBIG,XHUGE,XMAX/26.543D0,6.71D7,2.53D307/
  139. XINF=XGRAND
  140. XNEG=-26d0
  141. XSMALL=xzprec
  142. XBIG =26d0
  143. XHUGE=1d0/(2d0*sqrt(XSMALL))
  144. XMAX=min(xinf,(1d0/(sqrt(Xpi)*Xpetit)))
  145. C------------------------------------------------------------------
  146. C Coefficients for approximation to erf in first interval
  147. C------------------------------------------------------------------
  148. DATA A/3.16112374387056560D00,1.13864154151050156D02,
  149. 1 3.77485237685302021D02,3.20937758913846947D03,
  150. 2 1.85777706184603153D-1/
  151. DATA B/2.36012909523441209D01,2.44024637934444173D02,
  152. 1 1.28261652607737228D03,2.84423683343917062D03/
  153. C------------------------------------------------------------------
  154. C Coefficients for approximation to erfc in second interval
  155. C------------------------------------------------------------------
  156. DATA C/5.64188496988670089D-1,8.88314979438837594D0,
  157. 1 6.61191906371416295D01,2.98635138197400131D02,
  158. 2 8.81952221241769090D02,1.71204761263407058D03,
  159. 3 2.05107837782607147D03,1.23033935479799725D03,
  160. 4 2.15311535474403846D-8/
  161. DATA D/1.57449261107098347D01,1.17693950891312499D02,
  162. 1 5.37181101862009858D02,1.62138957456669019D03,
  163. 2 3.29079923573345963D03,4.36261909014324716D03,
  164. 3 3.43936767414372164D03,1.23033935480374942D03/
  165. C------------------------------------------------------------------
  166. C Coefficients for approximation to erfc in third interval
  167. C------------------------------------------------------------------
  168. DATA P/3.05326634961232344D-1,3.60344899949804439D-1,
  169. 1 1.25781726111229246D-1,1.60837851487422766D-2,
  170. 2 6.58749161529837803D-4,1.63153871373020978D-2/
  171. DATA Q/2.56852019228982242D00,1.87295284992346047D00,
  172. 1 5.27905102951428412D-1,6.05183413124413191D-2,
  173. 2 2.33520497626869185D-3/
  174. C------------------------------------------------------------------
  175. X = ARG
  176. Y = ABS(X)
  177. IF (Y .LE. THRESH) THEN
  178. C------------------------------------------------------------------
  179. C Evaluate erf for |X| <= 0.46875
  180. C------------------------------------------------------------------
  181. YSQ = ZERO
  182. IF (Y .GT. XSMALL) YSQ = Y * Y
  183. XNUM = A(5)*YSQ
  184. XDEN = YSQ
  185. DO 20 I = 1, 3
  186. XNUM = (XNUM + A(I)) * YSQ
  187. XDEN = (XDEN + B(I)) * YSQ
  188. 20 CONTINUE
  189. RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
  190. IF (JINT .NE. 0) RESULT = ONE - RESULT
  191. IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
  192. GO TO 800
  193. C------------------------------------------------------------------
  194. C Evaluate erfc for 0.46875 <= |X| <= 4.0
  195. C------------------------------------------------------------------
  196. ELSE IF (Y .LE. FOUR) THEN
  197. XNUM = C(9)*Y
  198. XDEN = Y
  199. DO 120 I = 1, 7
  200. XNUM = (XNUM + C(I)) * Y
  201. XDEN = (XDEN + D(I)) * Y
  202. 120 CONTINUE
  203. RESULT = (XNUM + C(8)) / (XDEN + D(8))
  204. IF (JINT .NE. 2) THEN
  205. YSQ = AINT(Y*SIXTEN)/SIXTEN
  206. DEL = (Y-YSQ)*(Y+YSQ)
  207. RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
  208. END IF
  209. C------------------------------------------------------------------
  210. C Evaluate erfc for |X| > 4.0
  211. C------------------------------------------------------------------
  212. ELSE
  213. RESULT = ZERO
  214. IF (Y .GE. XBIG) THEN
  215. IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300
  216. IF (Y .GE. XHUGE) THEN
  217. RESULT = SQRPI / Y
  218. GO TO 300
  219. END IF
  220. END IF
  221. YSQ = ONE / (Y * Y)
  222. XNUM = P(6)*YSQ
  223. XDEN = YSQ
  224. DO 240 I = 1, 4
  225. XNUM = (XNUM + P(I)) * YSQ
  226. XDEN = (XDEN + Q(I)) * YSQ
  227. 240 CONTINUE
  228. RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
  229. RESULT = (SQRPI - RESULT) / Y
  230. IF (JINT .NE. 2) THEN
  231. YSQ = AINT(Y*SIXTEN)/SIXTEN
  232. DEL = (Y-YSQ)*(Y+YSQ)
  233. RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
  234. END IF
  235. END IF
  236. C------------------------------------------------------------------
  237. C Fix up for negative argument, erf, etc.
  238. C------------------------------------------------------------------
  239. 300 IF (JINT .EQ. 0) THEN
  240. RESULT = (HALF - RESULT) + HALF
  241. IF (X .LT. ZERO) RESULT = -RESULT
  242. ELSE IF (JINT .EQ. 1) THEN
  243. IF (X .LT. ZERO) RESULT = TWO - RESULT
  244. ELSE
  245. IF (X .LT. ZERO) THEN
  246. IF (X .LT. XNEG) THEN
  247. RESULT = XINF
  248. ELSE
  249. YSQ = AINT(X*SIXTEN)/SIXTEN
  250. DEL = (X-YSQ)*(X+YSQ)
  251. Y = EXP(YSQ*YSQ) * EXP(DEL)
  252. RESULT = (Y+Y) - RESULT
  253. END IF
  254. END IF
  255. END IF
  256. 800 RETURN
  257. C---------- Last card of CALERF ----------
  258. END
  259.  
  260.  
  261.  
  262.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales