PROGRAM ZUFALL
**********************************************************************
*     Berechnung einer Zufallskette aus Einsen und Nullen
*********************************************************************
 
      REAL R(97)
      EXTERNAL RAN

C     IBIT(X) = 0 für X < 0.5, IBIT(X) = 1 für X > 0.5
      IBIT(X) = INT(2.*X)

      PRINT *, 
     &     'Eingabe Länge der Zufallsfolge , Anfangszustand ISEED'
      PRINT *, 
     &     ' (2 natuerliche Zahlen):'
      READ *,M,ISEED


      PRINT *, 
     &     M,' zufaellige Werte von 0 und 1:'

C     RAN Zufallszahl zwischen 0 und 1
      DO 100, I=1,M
          PRINT *, IBIT(RAN(ISEED,IX1,IX2,IX3,R))

 100  CONTINUE
      
      END       


*
*
      REAL FUNCTION RAN(IDUM,IX1,IX2,IX3,R)
**********************************************************************
*     Returns a random deviate between 0.0 and 1.0. Set IDUM to
*     any negative value to initialize or reinitialize the sequence.
*     Numerical Recipes p.196
**********************************************************************
      REAL R(97)
      INTEGER  IX1,IX2,IX3,IC1,IC2,IC3,IA1,IA2,IA3,IFF
      PARAMETER (M1=259200,IA1=7141,IC1=54733,RM1=1./M1)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)

      DATA IFF /0/

      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
         IFF=1
         IX1=MOD(IC1-IDUM,M1)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX2=MOD(IX1,M2)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX3=MOD(IX1,M3)
         DO 11, J=1,97
           IX1=MOD(IA1*IX1+IC1,M1)
           IX2=MOD(IA2*IX2+IC2,M2)
           R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
 11      CONTINUE
         IDUM=1

       ENDIF

       IX1=MOD(IA1*IX1+IC1,M1)
       IX2=MOD(IA2*IX2+IC2,M2)
       IX3=MOD(IA3*IX3+IC3,M3)
       
       J=1+(97*IX3)/M3
       
*       IF (J.GT.97.OR.J.LT.1) PAUSE

       RAN=R(J)  


       R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1

       RETURN
       END
*
*