This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / erf.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.2  1997/07/01 15:31:41  mclareni
6 * Correction for FP exceptions and restriction of calculations to realistic values, from M. Schroder
7 *
8 * Revision 1.1.1.1  1996/04/01 15:01:53  mclareni
9 * Mathlib gen
10 *
11 *
12 #include "gen/pilot.h"
13 #if defined(CERNLIB_DOUBLE)
14       FUNCTION ERF(X)
15
16       LOGICAL LEF
17       DIMENSION P2(0:4),Q2(0:4)
18
19       PARAMETER(Z1 = 1, HF = Z1/2, C1 = 0.56418 958)
20
21 C     Above the value of VMAX any calculation is pointless. The value is
22 C     choosen with a big safety margin - even the double precision
23 C     version only returns 1. for V (=ABS(X)) >= 5.9
24       PARAMETER(VMAX = 7.)
25 C     The value for SWITCH is badly chosen for the single precision
26 C     version, which returns 1. already for V >= 3.9
27       PARAMETER(SWITCH = 4.)
28
29       DATA P10,Q10,P11 /+3.67678 77, +3.25845 93, -9.79704 65E-2/
30
31       DATA (P2(I),Q2(I),I=0,4)
32      +/+7.37388 83E+0, +7.37396 09E+0, +6.86501 85E+0, +1.51849 08E+1,
33      1 +3.03179 93E+0, +1.27955 30E+1, +5.63169 62E-1, +5.35421 68E+0,
34      2 +4.31877 87E-5, +1/
35
36       DATA P30,Q30,P31 /-1.24368 54E-1, +4.40917 06E-1, -9.68210 36E-2/
37
38       LEF=.TRUE.
39       GO TO 9
40
41       ENTRY ERFC(X)
42       LEF=.FALSE.
43
44     9 V=ABS(X)
45       IF(V .LT. HF) THEN
46        Y=V**2
47        H=X*(P10+P11*Y)/(Q10+Y)
48        HC=1-H
49       ELSE
50        IF(V .LT. SWITCH) THEN
51         AP=P2(4)
52         AQ=Q2(4)
53         DO 2 I = 3,0,-1
54         AP=P2(I)+V*AP
55     2   AQ=Q2(I)+V*AQ
56         HC=EXP(-V**2)*AP/AQ
57         H=1-HC
58        ELSEIF ( V .LT. VMAX) THEN
59         Y=1/V**2
60         HC=EXP(-V**2)*(C1+Y*(P30+P31*Y)/(Q30+Y))/V
61         H=1-HC
62 C     for very big values we can save us any calculation, and the
63 C     FP-exceptions we would get from EXP.
64        ELSE
65         H  = 1.
66         HC = 0.
67        ENDIF
68        IF(X .LT. 0) THEN
69         H=-H
70         HC=2-HC
71        ENDIF
72       ENDIF
73       IF(LEF) THEN
74        ERF=H
75       ELSE
76        ERFC=HC
77       ENDIF
78       RETURN
79       END
80 #endif