MINICERN/mathlib/gen/g/gausin64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:43  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if !defined(CERNLIB_DOUBLE)
  11       FUNCTION  GAUSIN(P)
  12 #endif
  13 #if defined(CERNLIB_DOUBLE)
  14       FUNCTION DGAUSN(P)
  15 #include "gen/imp64.inc"
  16 #endif
  17 C     Computes a "Normal Deviate"
  18 C     Based on G.W. Hill & A.W. Davis, Algorithm 442 Normal Deviate
  19 C     Collected Algorithms from CACM
  20
  21       CHARACTER NAME*(*)
  22       CHARACTER*80 ERRTXT
  23 #if !defined(CERNLIB_DOUBLE)
  24       PARAMETER (NAME = 'GAUSIN')
  25 #endif
  26 #if defined(CERNLIB_DOUBLE)
  27       PARAMETER (NAME = 'DGAUSN')
  28 #endif
  29
  30       PARAMETER (C = 2.50662 82746 31000 50D0)
  31       PARAMETER (Z1 = 1, HF = Z1/2, C1 = 3*Z1/4, C2 = 7*Z1/8, C3 = Z1/3)
  32
  33       IF(P .LE. 0 .OR. P .GE. 1) THEN
  34        H=0
  35        WRITE(ERRTXT,101) P
  36        CALL MTLPRT(NAME,'G105.1',ERRTXT)
  37       ELSEIF(P .EQ. HF) THEN
  38        H=0
  39       ELSE
  40        X=P
  41        IF(P .GT. HF) X=1-P
  42        X=SQRT(-2*LOG(X))
  43        X=X-((7.47395*X+494.877)*X+1637.720)/
  44      1     (((X+117.9407)*X+908.401)*X+659.935)
  45        IF(P .LT. HF) X=-X
  46        S=X**2
  47 #if !defined(CERNLIB_DOUBLE)
  48        Z=C*(P-FREQ(X))*EXP(HF*S)
  49 #endif
  50 #if defined(CERNLIB_DOUBLE)
  51        Z=C*(P-DFREQ(X))*EXP(HF*S)
  52 #endif
  53        H=(((((C1*S+C2)*Z+X)*X+HF)*C3*Z+HF*X)*Z+1)*Z+X
  54       ENDIF
  55 #if !defined(CERNLIB_DOUBLE)
  56        GAUSIN=H
  57 #endif
  58 #if defined(CERNLIB_DOUBLE)
  59       DGAUSN=H
  60 #endif
  61       RETURN
  62   101 FORMAT('ARGUMENT P =',1P,D15.8,' NOT IN RANGE')
  63       END
  64