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