Logo Search packages:      
Sourcecode: foptions version File versions

D1-LowDiscrepancy.f


C PART I:  HALTON SEQUENCE
C PART II: SOBOL SEQUENCE


C ##############################################################################
C PART I: HALTON SEQUENCE:


C This library is free software; you can redistribute it and/or
C modify it under the terms of the GNU Library General Public
C License as published by the Free Software Foundation; either
C version 2 of the License, or (at your option) any later version.
C
C This library is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
C GNU Library General Public License for more details.
C
C You should have received a copy of the GNU Library General 
C Public License along with this library; if not, write to the 
C Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, 
C MA  02111-1307  USA


COPYRIGHT: DIETHELM WUERTZ, SEPT. 2002


C-------------------------------------------------------------------------------

C     INITHALTON (DIMEN, QUASI, BASE, OFFSET)
C     NEXTHALTON (DIMEN, QUASI, BASE, OFFSET)
C     HALTON (QN, N, DIMEN, QUASI, BASE, OFFSET, INIT, TRANSFORM)
C     REAL*8 FUNCTION HQNORM(P)

C-------------------------------------------------------------------------------


      SUBROUTINE INITHALTON(DIMEN, QUASI, BASE, OFFSET)

C     INITIALIZE THE HALTON LOW DISCREPANCY SEQUENCE.
C     THE BASE IS CALCULATED FROM PRIMES

      INTEGER DIMEN, BASE(DIMEN), ITER(DIMEN), OFFSET, DIGIT
      REAL*8 QUASI(DIMEN), HALF
      INTRINSIC MOD
 
C     INIT BASE FROM PRIMES - THIS IMPLEMENTS A SIMMPLE SIEVE:
      BASE(1) = 2
      BASE(2) = 3
      N = 3
      NC = 2
      DO WHILE(NC.LT.DIMEN)
        M = N/2
        K = 0
        IF (MOD(N,2).NE.0.AND.MOD(N,3).NE.0) THEN
           DO I = 5, M
               IF(MOD(N,I).EQ.0) K = K + 1
         ENDDO
           IF (K.EQ.0) THEN
              NC = NC + 1
              BASE(NC) = N
           ENDIF
        ENDIF
        N = N + 1
      ENDDO
      
C     NOW CREATE THE FIRST QUASI RANDOM NUMBER:
      OFFSET = 0
      DO NB = 1, DIMEN        
        ITER(NB) = OFFSET
        QUASI(NB) = 0.0D0
        HALF = 1.0D0 / BASE(NB)
        DO WHILE (ITER(NB).NE.0)
           DIGIT = MOD ( ITER(NB), BASE(NB) )
           QUASI(NB) = QUASI(NB) + DIGIT * HALF
           ITER(NB) = ( ITER(NB) - DIGIT ) / BASE(NB)
           HALF = HALF / BASE(NB)
        ENDDO 
      ENDDO

C     SET THE COUNTER:
      OFFSET = OFFSET + 1
      
      RETURN
      END
  

C-------------------------------------------------------------------------------

  
      SUBROUTINE NEXTHALTON(DIMEN, QUASI, BASE, OFFSET) 

C     GENERATE THE NEXT POINT IN HALTON





























































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































































Generated by  Doxygen 1.6.0   Back to index