MINICERN/mathlib/gen/v/rngama.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:56  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       FUNCTION RNGAMA(P)
  11       EXTERNAL RANLUX
  12
  13       DIMENSION STOR(15)
  14
  15       H=0
  16       IF(P .GT. 15) THEN
  17     1  CALL RNORMX(A,1,RANLUX)
  18        H=P*(1-1/(9*P)+A/(3*SQRT(P)))**3
  19        IF(H .LE. 0) GO TO 1
  20       ELSE
  21        M=P
  22        F=P-M
  23        IF(M .GT. 0) THEN
  24         X=1
  25         CALL RANLUX(STOR,M)
  26         DO 2 I = 1,M
  27     2   X=X*STOR(I)
  28         H=-LOG(X)
  29        ENDIF
  30        IF(F .GE. 0.00001) THEN
  31         CALL RANLUX(X,1)
  32         X1=-LOG(X)
  33         IF(F .GE. 0.9999) THEN
  34          H=H+X1
  35         ELSE
  36          CALL RANLUX(X,1)
  37     3    WLOG=LOG(X)/F
  38          IF(WLOG .GE. -100) THEN
  39           W1=EXP(WLOG)
  40           CALL RANLUX(X,1)
  41           WLOG=LOG(X)/(1-F)
  42           IF(WLOG .LT. -100) THEN
  43            H=H+X1
  44           ELSE
  45            W=W1+EXP(WLOG)
  46            IF(W .GT. 1) GO TO 3
  47            H=H+X1*W1/W
  48           ENDIF
  49          ENDIF
  50         ENDIF
  51        ENDIF
  52       ENDIF
  53       RNGAMA=H
  54       RETURN
  55       END