← Back to branch summary

~maddevelopers/mg5amcnlo/WWW_tg

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_242195/PROC_242195/Source/htuple.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_242195/PROC_242195/Source/htuple.f
 
1
      subroutine ntuple(x,a,b,ii,jconfig)
 
2
C--------------------------------------------------------------
 
3
c
 
4
c
 
5
c     This is a modified version for use with multi-pole integrations
 
6
c     it allows you to get the same set of random numbers several times
 
7
c     (1 for each configuration) jj tells it what configuration to use, 
 
8
c     so if you always put jj=1, you basically get the old version of 
 
9
c     ntuple out. It is currently configued to have maxconfig=25.
 
10
c     Modified by Tim Stelzer July 5 1995
 
11
c
 
12
c
 
13
c       Abstract:
 
14
c
 
15
c      ntuple - returns sequences of quasi-uniform random numbers
 
16
c               in the interval (a,b) using scrambled radical inverse
 
17
c               sequences. These numbers are designed for use in
 
18
c               Monte Carlo integration routines.
 
19
c
 
20
c  Author: Bill Long, UW-Madison Phenomenology Institute, 19-MAR-1991
 
21
c               Based on the algorithm for an earlier routine, htuple,
 
22
c               written by George Weller.
 
23
c
 
24
c               Original reference for this algorithm:
 
25
c
 
26
c         E. Braaten and G. Weller, J. Comp. Phys. 33,249-258 (1979)
 
27
c
 
28
c  References for radical inverse sequences (and Halton sequences)
 
29
c               are ck +^Z7ed in the above article.
 
30
c
 
31
c    Usage Notes:
 
32
c
 
33
c 1) Sequence numbers, i, range from 1 to 25 and correspond
 
34
c   to scrambled radical inverse sequences based on the first
 
35
c   25 primes (2..97).  In general, when performing a multi-dimensional
 
36
c   integral, a separate value of i should be used for the values
 
37
c   along each axis.
 
38
c
 
39
c 2) The basic algorithm generates values in the range (0,1). These
 
40
c    are rescaled to the range (a,b) in the final statement, so
 
41
c    it is not necessary to have a < b.  If a=b, the returned value, x,
 
42
c    will always be equal to a.
 
43
c
 
44
c 3) The sequences for different values of i are independent. They
 
45
c    cycle with different periods ranging from ~ 4M numbers 
 
46
c    through 147M numbers, covering numbers of precision 23-27 bits
 
47
c    appropriate for the mantissa of a single precision real value.
 
48
c
 
49
c 4) Ntuple differs from htuple primarily in that ntuple sequences
 
50
c    cycle with periods that are roughly equal for each value of i.
 
51
c    Htuple sequences cycled too quickly for small i, and too
 
52
c    slowly for large i. Ntuple is also written in a more modern
 
53
c    style, and uses considerably less memory, helping execution
 
54
c    speed on cache-sensitive machines.
 
55
c
 
56
c
 
57
c--Argument Declarations
 
58
c
 
59
      real*8  x    !  OUT - quasi-random value returned
 
60
                   !        x is in the range (a,b)
 
61
      real*8  a    !  IN  - Lower bound of interval for x
 
62
      real*8  b    !  IN  - Upper bound of interval for x
 
63
      integer ii   !  IN  - Sequence number, restricted to 1 <= i <= MaxDim
 
64
      integer jconfig!IN  - Pole number, restricted to1<=jconfig<= MaxConfigs
 
65
c
 
66
c     Constants
 
67
c
 
68
      include "genps.inc"
 
69
      include 'maxconfigs.inc'
 
70
 
 
71
      integer ndim,kdim,mdim,maxconfig,ktot
 
72
      parameter (ndim = maxdim, kdim = 181, mdim = 1060)
 
73
c      parameter (maxconfig=maxconfigs, ktot=kdim*maxconfig)
 
74
      parameter (maxconfig=lmaxconfigs, ktot=kdim*maxconfig)
 
75
 
 
76
c
 
77
c--Local Variable Declarations
 
78
c
 
79
      integer base_minus1(ndim),mix(mdim),k(kdim,maxconfig)
 
80
      double precision accum(kdim,maxconfig)
 
81
      double precision pbase(kdim)
 
82
      integer offset,koffset(ndim),mix_offset(ndim),maxj(ndim)
 
83
      integer jj
 
84
      logical first_time
 
85
 
 
86
c
 
87
c--Fixed Data Initializations
 
88
c
 
89
      data (pbase(i), i = 1, 146)
 
90
     .     /2d0,4d0,8d0,16d0,32d0,64d0,128d0,256d0,512d0,1024d0,
 
91
     .     2048d0,4096d0,8192d0,16384d0,32768d0,65536d0,131072d0,
 
92
     .     262144d0,524288d0,1048576d0,2097152d0,4194304d0,8388608d0,
 
93
     .     16777216d0,0d0,3d0,9d0,27d0,81d0,243d0,729d0,2187d0,6561d0,
 
94
     .     19683d0,59049d0,177147d0,531441d0,1594323d0,4782969d0,
 
95
     .     14348907d0,0d0,5d0,25d0,125d0,625d0,3125d0,15625d0,78125d0,
 
96
     .     390625d0,1953125d0,9765625d0,0d0,7d0,49d0,343d0,2401d0,
 
97
     .     16807d0,117649d0,823543d0,5764801d0,0d0,11d0,121d0,1331d0,
 
98
     .     14641d0,161051d0,1771561d0,19487171d0,0d0,13d0,169d0,2197d0,
 
99
     .     28561d0,371293d0,4826809d0,0d0,17d0,289d0,4913d0,83521d0,
 
100
     .     1419857d0,24137569d0,0d0,19d0,361d0,6859d0,130321d0,
 
101
     .     2476099d0,47045881d0,0d0,23d0,529d0,12167d0,279841d0,
 
102
     .     6436343d0,0d0,29d0,841d0,24389d0,707281d0,20511149d0,0d0,
 
103
     .     31d0,961d0,29791d0,923521d0,28629151d0,0d0,37d0,1369d0,
 
104
     .     50653d0,1874161d0,69343957d0,0d0,41d0,1681d0,68921d0,
 
105
     .     2825761d0,115856201d0,0d0,43d0,1849d0,79507d0,3418801d0,
 
106
     .     147008443d0,0d0,47d0,2209d0,103823d0,4879681d0,0d0,53d0,
 
107
     .     2809d0,148877d0,7890481d0,0d0,59d0,3481d0,205379d0,
 
108
     .     12117361d0,0d0,61d0,3721d0,226981d0,13845841d0,0d0/
 
109
      data (pbase(i), i=147, 181)
 
110
     .     /67d0,4489d0,300763d0,20151121d0,0d0,71d0,5041d0,357911d0,
 
111
     .     25411681d0,0d0,73d0,5329d0,389017d0,28398241d0,0d0,79d0,
 
112
     .     6241d0,493039d0,38950081d0,0d0,83d0,6889d0,571787d0,
 
113
     .     47458321d0,0d0,89d0,7921d0,704969d0,62742241d0,0d0,97d0,
 
114
     .     9409d0,912673d0,88529281d0, 0d0/
 
115
 
 
116
      data base_minus1/
 
117
     .      1,    2,    4,    6,   10,   12,   16,   18,   22,   28,
 
118
     .     30,   36,   40,   42,   46,   52,   58,   60,   66,   70,
 
119
     .     72,   78,   82,   88,   96/
 
120
 
 
121
      data maxj/ 24, 15, 10, 8,    7,    6,    6,    6,    5,    5,
 
122
     .     5,    5,    5,    5,    4,    4,    4,    4,    4,    4,
 
123
     .     4,    4,    4,    4,    4/
 
124
 
 
125
      data koffset/
 
126
     .     0,   25,   41,   52,   61,   69,   76,   83,   90,   96,
 
127
     .     102,  108,  114,  120,  126,  131,  136,  141,  146,  151,
 
128
     .     156,  161,  166,  171,  176/  
 
129
 
 
130
      data mix_offset/
 
131
     .     0,    2,    5,   10,   17,   28,   41,   58,   77,  100,
 
132
     .     129,  160,  197,  238,  281,  328,  381,  440,  501,  568,
 
133
     .     639,  712,  791,  874,  963/
 
134
 
 
135
      data (mix(i), i = 1, 412)
 
136
     .     /1,0, 1,2,0, 3,1,4,2,0, 4,2,6,1,5,3,0, 5,8,2,10,3,6,1,9,
 
137
     .     7,4,0, 6,10,2,8,4,12,1,9,5,11,3,7,0, 8,13,3,11,5,16,1,10,7,
 
138
     .     14,4,12,2,15,6,9,0, 9,14,3,17,6,11,1,15,7,12,4,18,8,2,16,10,
 
139
     .     5,13,0, 11,17,4,20,7,13,2,22,9,15,5,18,1,14,10,21,6,16,3,19,
 
140
     .     8,12,0, 15,7,24,11,20,2,27,9,18,4,22,13,26,5,16,10,23,1,19,
 
141
     .     28,6,14,17,3,25,12,8,21,0,  15,23,5,27,9,18,2,29,12,20,7,25,
 
142
     .     11,17,3,30,14,22,1,21,8,26,10,16,28,4,19,6,24,13,0, 18,28,6,
 
143
     .     23,11,34,3,25,14,31,8,20,36,1,16,27,10,22,13,32,4,29,17,7,
 
144
     .     35,19,2,26,12,30,9,24,15,33,5,21,0, 20,31,7,26,12,38,3,23,
 
145
     .     34,14,17,29,5,40,10,24,1,35,18,28,9,33,15,21,4,37,13,30,8,
 
146
     .     39,22,2,27,16,32,11,25,6,36,19,0, 21,32,7,38,13,25,3,35,17,
 
147
     .     28,10,41,5,23,30,15,37,1,19,33,11,26,42,8,18,29,4,39,14,22,
 
148
     .     34,6,24,12,40,2,31,20,27,9,36,16,0, 24,12,39,6,33,20,44,3,
 
149
     .     29,16,36,10,42,22,8,31,26,14,46,1,35,18,28,5,40,19,37,11,25,
 
150
     .     43,4,30,15,34,9,45,21,2,32,17,41,13,27,7,38,23,0,
 
151
     .     26,40,9,33,16,49,4,36,21,45,12,29,6,51,23,38,14,43,1,30,19,
 
152
     .     47,10,34,24,42,3,27,52,15,18,39,7,46,31,11,35,20,48,2,28,41,
 
153
     .     8,22,50,13,32,17,44,5,37,25,0, 29,44,10,52,18,34,4,48,23,38,
 
154
     .     13,57,7,32,41,20,54,2,26,46,15,36,24,50,8,40,16,56,5,30,43/
 
155
      data (mix(i), i = 413, 803)
 
156
     .     /21,51,11,33,1,58,27,37,14,47,19,28,45,6,53,12,35,22,42,3,
 
157
     .     55,25,31,9,49,17,39,0,30,46,10,38,18,56,4,42,24,52,14,33,21,
 
158
     .     59,6,40,27,49,2,35,16,54,12,44,26,50,8,32,58,19,1,41,29,48,
 
159
     .     13,36,22,60,7,45,23,53,9,34,17,55,3,39,28,47,15,37,20,57,5,
 
160
     .     43,25,51,11,31,0, 33,50,11,59,20,39,5,54,26,44,15,64,23,36,
 
161
     .     2,57,30,47,9,62,18,41,13,52,28,37,4,66,24,46,8,55,31,17,60,
 
162
     .     34,1,48,21,43,63,12,38,25,53,7,49,16,58,29,6,42,65,19,35,10,
 
163
     .     51,27,56,3,40,32,61,14,45,22,0, 35,53,12,62,21,41,5,67,28,
 
164
     .     46,16,56,25,8,50,38,65,2,32,59,19,44,14,70,30,48,7,39,58,22,
 
165
     .     10,63,33,26,52,1,55,18,43,68,13,36,47,4,61,24,40,29,66,9,51,
 
166
     .     17,57,23,37,3,69,31,45,15,60,11,49,34,20,64,6,54,27,42,0,
 
167
     .     36,55,12,46,22,67,5,41,61,18,30,52,8,70,27,43,15,59,33,2,64,
 
168
     .     38,24,50,10,72,20,48,31,57,4,63,25,40,14,54,35,68,7,45,17,
 
169
     .     60,28,1,66,39,21,51,11,71,32,47,13,56,26,44,3,65,34,19,58,9,
 
170
     .     49,37,69,16,29,53,6,62,23,42,0, 39,59,13,69,24,46,6,74,31,
 
171
     .     51,18,63,9,42,55,27,77,2,35,65,21,48,15,71,33,53,4,61,29,43,
 
172
     .     17,75,37,10,67,49,22,57,7,72,26,40,56,1,64,30,45,14,78,20,
 
173
     .     52,34,11,68,41,60,5,36,73,23,50,16,62,28,3,76,44,25,58,12,
 
174
     .     66,38,19,54,32,70,8,47,0,41,62,14,73,25,48,6,67,32,54,19,80/
 
175
      data (mix(i), i = 804, 1060)
 
176
     .     /10,44,58,29,76,2,37,64,22,51,16,71,35,56,8,82,27,46,12,69,
 
177
     .     39,60,4,50,24,78,31,65,17,42,74,1,53,21,61,34,11,79,43,28,
 
178
     .     68,7,55,38,75,15,47,20,70,5,57,33,81,26,49,9,63,36,66,18,45,
 
179
     .     3,77,30,59,23,52,13,72,40,0, 44,67,15,56,27,82,6,50,74,22,
 
180
     .     36,63,10,86,33,53,18,77,40,2,70,47,29,80,12,60,38,65,20,88,
 
181
     .     4,51,31,72,24,58,8,78,42,46,16,84,34,62,1,69,26,55,19,76,41,
 
182
     .     11,83,49,30,66,7,59,37,87,14,54,25,73,21,68,43,3,79,35,57,
 
183
     .     13,81,45,28,64,5,75,32,52,17,85,39,9,61,71,23,48,0, 48,73,
 
184
     .     16,61,29,89,7,55,81,34,22,69,41,94,3,52,77,19,38,85,12,64,
 
185
     .     44,26,91,58,9,71,32,79,14,50,66,24,96,1,46,83,36,59,18,75,
 
186
     .     30,87,5,54,42,68,21,92,10,62,39,80,27,56,6,86,47,72,15,35,
 
187
     .     93,43,65,2,76,25,53,84,17,37,67,11,90,49,31,74,20,60,95,4,
 
188
     .     45,63,28,82,13,57,40,78,8,88,33,51,23,70,0/
 
189
c
 
190
c--Variable Data Initializations
 
191
c
 
192
      data accum /ktot*0d0/
 
193
 
 
194
      data k /ktot*0d0/
 
195
      data first_time /.true./
 
196
c
 
197
c--Code:
 
198
c
 
199
      if (first_time) then
 
200
         write(*,*) 'Warning htuple modified for 1 configuration only'
 
201
         write(*,*) 'Using htuple configuration ',jconfig
 
202
         first_time=.false.
 
203
c     
 
204
c        to use multiple configurations need to use line
 
205
        jj = jconfig
 
206
c
 
207
      endif
 
208
      jj = jconfig
 
209
c      jj=1      !use jj=jconfig for multiconfiguraion mode
 
210
      if (jj .lt. 1 .or. jj .gt. maxconfig) then
 
211
         print*,'Error in ntuple.  Invalid pole choice',jj
 
212
         stop
 
213
      endif
 
214
      i=ii
 
215
      i = ii+jj                   !This keeps us from generating same ran #'s
 
216
      do while (i .gt. ndim)
 
217
         i=i-ndim
 
218
      enddo
 
219
c      if (i .gt. ndim) i=i-ndim   !For different configurations
 
220
      j = 1
 
221
      offset = koffset(i)
 
222
      do while (k(offset+j,jj) .eq. base_minus1(i))
 
223
         k(offset+j,jj)=0
 
224
         j=j+1
 
225
         if (j.gt.maxj(i)) then
 
226
            do j=1,maxj(i)
 
227
               k(offset+j,jj)=0
 
228
               accum(offset+j,jj)=0d0
 
229
            end do
 
230
            j=1
 
231
         end if
 
232
      end do
 
233
 
 
234
      k(offset+j,jj)=k(offset+j,jj)+1
 
235
      accum(offset+j,jj) = accum(offset+j+1,jj) +
 
236
     .     mix(mix_offset(i)+k(offset+j,jj))/pbase(offset+j)
 
237
      do  jjj=2,j-1
 
238
         accum(offset+jjj,jj) = accum(offset+j,jj)
 
239
      end do
 
240
      x = a + (b-a) * accum(offset+j,jj)
 
241
c      write(*,'(2i6,1f15.8)') jj,i,x
 
242
      end
 
243