C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS(M)
C                    Uses RANDd() + 5 random numbers with correction
C                    Old code used in many simulations
C   N=1e8  sigma= 0.99556     mean= 0.10542E-03
c  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1589     0.1589     0.1587     15894145 15889594
c  1.50 0.6704E-01 0.6702E-01 0.6681E-01  6703736  6702430
c  2.00 0.2258E-01 0.2256E-01 0.2275E-01  2258038  2256030
c  2.50 0.5841E-02 0.5843E-02 0.6210E-02   584148   584264
c  3.00 0.1031E-02 0.1028E-02 0.1350E-02   103141   102783
c  4.00 0.5900E-06 0.6100E-06 0.3170E-04       59       61
C
C--------------------------------------
      X=0.
      DO  I=1,5
         X=X+RANDd(M)
      EndDo
      X2   =1.5491933*(X-2.5)
      GAUSS=X2*(1.-0.01*(3.-X2*X2))
      RETURN
      END
C--------------------------------------
C		normal random numbers
      FUNCTION GAUSS2(M)
C                        Uses ranndd for homogenous rand numbers
C                          sum of 12 random numbers + corrections
C                          It is  slow: very (10 times than GAUSS3)
C                          Quality is very good: upto 4 sigma
c  sigma Frac(>sig) Frac(<-sig)   True       n(>)    n(<)
c  1.00 0.1586     0.1586     0.1587     15857116 15857711
c  1.50 0.6678E-01 0.6680E-01 0.6681E-01  6678110  6679725
c  2.00 0.2277E-01 0.2275E-01 0.2275E-01  2277430  2274868
c  2.50 0.6197E-02 0.6189E-02 0.6210E-02   619732   618895
c  3.00 0.1327E-02 0.1319E-02 0.1350E-02   132713   131929
c  4.00 0.2664E-04 0.2727E-04 0.3170E-04     2664     2727
C                   
c-------------------------------------------------------------------- 
      X=0.
      DO  I=1,12
         X=X+Randd(M)
      EndDo
      X2 = X-6.0
      GAUSS2   =X2 *(1.-0.0045*(3.-X2*X2))
      RETURN
      END

C------------------------------------------------
C                                              random number generator
      FUNCTION RANDd(M)
C------------------------------------------------
      DATA LC,AM,KI,K1,K2,K3,K4,L1,L2,L3,L4
     +	/453815927,2147483648.,2147483647,536870912,131072,256,
     +	 16777216,4,16384,8388608,128/
      ML=M/K1*K1
      M1=(M-ML)*L1
      ML=M/K2*K2
      M2=(M-ML)*L2
      ML=M/K3*K3
      M3=(M-ML)*L3
      ML=M/K4*K4
      M4=(M-ML)*L4
      M5=KI-M
      IF(M1.GE.M5)M1=M1-KI-1
      ML=M+M1
      M5=KI-ML
      IF(M2.GE.M5)M2=M2-KI-1
      ML=ML+M2
      M5=KI-ML
      IF(M3.GE.M5)M3=M3-KI-1
      ML=ML+M3
      M5=KI-ML
      IF(M4.GE.M5)M4=M4-KI-1
      ML=ML+M4
      M5=KI-ML
      IF(LC.GE.M5)ML=ML-KI-1
      M=ML+LC
      RANDd=M/AM
      RETURN
      END