C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-C
C======================================================================C
C                                                                      C
C     *********************************************************        C
C     RETOURNE UN REEL ALEATOIRE SUIVANT LA LOI UNIFORME U(A,B)        C
C     *********************************************************        C
C                                                                      C
C                                                                      C
C     ARGUMENTS D'ENTREE :                                             C
C     ------------------                                               C
C                                                                      C
C     UA     : BORNE INFERIEURE                                        C
C     UB     : BORNE SUPERIEURE                                        C
C                                                                      C
C              L'OBSERVATION RETOURNEE VERIFIE : UA < LOIUAB < UB      C
C                                                                      C
C======================================================================C
C                                                                      C
C     ALGORITHME :                                                     C
C                                                                      C
C         (1)    S   :=   ( A * S  + B )   MOD 2**M                    C
C         (2)    S NORMALISE EN FLOTTANT ENTRE 0 ET 1                  C
C                ( DIVISION REELLE DE S PAR 2**M )                     C
C         (3)    PASSAGE DE [0;1] A ]UA;UB[                            C
C                                                                      C
C     VALEURS ADOPTEES DANS LA ROUTINE G05CAF DU LOGICIEL NAG:         C
C                                                                      C
C     A = 13**13 ; B = 0 ; M = 59                                      C
C     INITIALISATION : S = 123456789*(1+2**32) ; PERIODE : 2**57       C
C                                                                      C
C======================================================================C
C                                                                      C
C     TRANSPORTABILITE ( SOURCE ET RESULTATS ) :                       C
C                                                                      C
C     A , B , ET  S , SONT REPRESENTES EN BASE 2 : 1 BIT PAR ELEMENT   C
C     DE TYPE INTEGER , AVEC AU MAXIMUM N BITS , ORDONNES DU PLUS      C
C     SIGNIFICATIF AU MOINS SIGNIFICATIF . CONSEQUENCE SUR LA          C
C     REPRESENTATION DES GRANDS ENTIERS ET SUR LEUR CONVERSION         C
C     EN DOUBLE PRECISION : ELLES NE SONT PAS DEPENDANT MACHINE .      C
C     PAR CONTRE , LA PRECISION DES OBSERVATIONS DEPEND DE LA          C
C     REPRESENTATION INTERNE DE L'ARITHMETIQUE DOUBLE PRECISION .      C
C                                                                      C
C     ATTENTION : LES REJETS DES OBSERVATIONS DUS AUX BORNES PEUVENT   C
C                 DIFFERER SUIVANT LA PRECISION DE L'ARITHMETIQUE      C
C                                                                      C
C     NOTE : ON NE PEUT DONC PAS AVOIR : M > N                         C
C                                                                      C
C======================================================================C
C                                                                      C
C     AUTRES PROGRAMMATIONS DE LOIUAB :                                C
C                                                                      C
C     AFIN DE DIMINUER LE COUT CPU , D'AUTRES VERSIONS DE LOIUAB       C
C     PEUVENT EXISTER ( EN ASSEMBLEUR NOTAMMENT ) ; ELLES ONT EN       C
C     COMMUN LE NOM DU COMMON ET DU BLOCK DATA , LA DECLARATION DE     C
C     FONCTION ET DE SES ARGUMENTS , ET DOIVENT RETOURNER LES MEMES    C
C     VALEURS SUR TOUTES LES MACHINES ( VOIR CI-AVANT ) ; LE COMMON    C
C     / UABLOI / ETANT STRICTEMENT INTERNE AU MODULE LOIUAB , LES      C
C     VARIABLES S , A , B , M , PEUVENT ETRE DECLAREES DIFFEREMENT ,   C
C     ET DES VARIABLES AUXILLIAIRES PEUVENT EXISTER DANS / UABLOI / .  C
C                                                                      C
C======================================================================C
C
      BLOCK DATA BDLUAB
C
      IMPLICIT INTEGER  ( A - Z )
C
      PARAMETER         ( N = 64 )
C
      COMMON / UABLOI / S (N) , A (N) , B (N) , M
C
      DATA              S / 0,0,0,0 , 0,1,1,1  ,  0,1,0,1 , 1,0,1,1 ,
     ,                      1,1,0,0 , 1,1,0,1  ,  0,0,0,1 , 0,1,0,1 ,
     ,                      0,0,0,0 , 0,1,1,1  ,  0,1,0,1 , 1,0,1,1 ,
     ,                      1,1,0,0 , 1,1,0,1  ,  0,0,0,1 , 0,1,0,1 /
C
      DATA              A / 0,0,0,0 , 0,0,0,0  ,  0,0,0,0 , 0,0,0,1 ,
     ,                      0,0,0,1 , 0,0,1,1  ,  0,1,1,1 , 0,1,1,0 ,
     ,                      1,0,0,1 , 1,0,1,1  ,  0,0,1,0 , 0,0,1,1 ,
     ,                      1,1,0,0 , 0,1,0,1  ,  1,1,1,1 , 1,1,0,1 /
C
      DATA              B / 0,0,0,0 , 0,0,0,0  ,  0,0,0,0 , 0,0,0,0 ,
     ,                      0,0,0,0 , 0,0,0,0  ,  0,0,0,0 , 0,0,0,0 ,
     ,                      0,0,0,0 , 0,0,0,0  ,  0,0,0,0 , 0,0,0,0 ,
     ,                      0,0,0,0 , 0,0,0,0  ,  0,0,0,0 , 0,0,0,0 /
C
      DATA              M / 59 /
C
      END
C
C======================================================================C
C
      DOUBLE PRECISION FUNCTION LOIUAB ( UA , UB )
C
      IMPLICIT INTEGER ( A - Z )
C
      DOUBLE PRECISION UA , UB ,
     ,                 X , DPI , LOIU01 ,
     ,                 ZERO , UN , UNDEMI
C
C
C ... EXTERNAL FORCE LE LINKER A CHERCHER LE BLOCK DATA
C     DANS LA BIBLIOTHEQUE OBJET
C
      EXTERNAL BDLUAB
C
C
      PARAMETER         ( N = 64 , N1 = N+1 )
C
      COMMON / UABLOI / S (N) , A (N) , B (N) , M
C
      INTEGER           ASB (N)
C
C
      DATA              BASE / 2 /
C
      DATA              ZERO , UN , UNDEMI / 0.D0 , 1.D0 , .5D0 /
C
C
C
C ... ASB = B + A * S ; ON ANTICIPE MODULO 2**M AVANT RECOPIE DANS S
C     --------------------------------------------------------------
C
      N1M    = N1 - M
      NN1M   = N + N1M
C
100   DO I = N1M , N
          ASB (I) = B (I)
      END DO
C
      DO J = N , N1M , -1
          IF ( S (J) .NE. 0 ) THEN
              DO I = N , NN1M-J , -1
                  K       = I + J - N
CC                ASB (K) = ASB (K) + A (I) * S (J)
                  ASB (K) = ASB (K) + A (I)
              END DO
          ENDIF
      END DO
C
      DO I = N , N1M+1 , -1
          RETENU = ASB (I) / BASE
          IF ( RETENU .GT. 0 ) THEN
              ASB (I)   = ASB (I) - BASE * RETENU
              ASB (I-1) = ASB (I-1) + RETENU
          ENDIF
      END DO
C
      RETENU = ASB (N1M) / BASE
      IF ( RETENU .GT. 0 ) ASB (N1M) = ASB (N1M) - BASE * RETENU
C
      DO I = N1M , N
          S (I) = ASB (I)
      END DO
C
C
C
C ... DIVISION REELLE : S / 2**M
C     --------------------------
C
      LOIU01 = ZERO
      DPI    = UN
C
      DO I = N1M , N
          DPI    = DPI * UNDEMI
          IF ( S (I) .NE. 0 ) LOIU01 = LOIU01 + DPI
      END DO
C
C
C
C ... OBSERVATION DE LA LOI UNIFORME : BORNES EXCLUES
C     -----------------------------------------------
C
      X      = UA + LOIU01 * ( UB - UA )
C
      IF ( X .LE. UA  .OR.  X .GE. UB ) GOTO 100
C
      LOIUAB = X
C
      RETURN
C
      END