Source code: riso/numerical/LowDiscrepency.java
1 package riso.numerical;
2
3 /** Translated from <tt>http://www.netlib.org/toms/659</tt>,
4 * an implementation of Sobol's low-discrepency sequence generator,
5 * by P Bratley and B L Fox.
6 * Described in <it>ACM Trans. on Mathematical Software</it>, vol. 14, no. 1, pp 88--100.
7 *
8 * <p> Note that this scheme works only in 2 or more dimensions -- it cannot generate
9 * a 1-dimensional sequence.
10 *
11 * <p> This file is distributed under the terms of the ACM Software Copyright and License Agreement.
12 * A copy of the license agreement, <tt>ACM-LICENSE.html</tt>, is included with the RISO distribution.
13 */
14 public class LowDiscrepency
15 {
16 static int [ ] primes =
17 {
18 1, 2, 3, 5, 5, 7, 7, 11, 11, 11, 11, 13, 13, 17, 17, 17, 17, 19, 19,
19 23, 23, 23, 23, 29, 29, 29, 29, 29, 29, 31, 31, 37, 37, 37, 37,
20 37, 37, 41, 41, 41
21 };
22
23 /** <tt>dimen</tt> must be at least 2.
24 */
25 public static void infaur ( boolean [ ] flag, int dimen, int atmost )
26 {
27 faure.s=dimen;
28 flag [ 1-1 ] = ( faure.s > 1 && faure.s < 41 );
29 if ( ! flag [ 1-1 ] ) return;
30 faure.qs=primes [ faure.s-1 ];
31 faure.testn=faure.qs*faure.qs*faure.qs*faure.qs;
32 faure.hisum= ( int ) ( Math.log ( atmost+faure.testn ) /Math.log ( faure.qs ) );
33 flag [ 2-1 ] = ( faure.hisum < 20 );
34 if ( ! flag [ 2-1 ] ) return;
35 faure.coef [ 0 ] [ 0 ] =1;
36 for ( int j = 1 ; j <= faure.hisum ; j++ )
37 {
38 faure.coef [ j ] [ 0 ] =1;
39 faure.coef [ j ] [ j ] =1;
40 }
41 for ( int j = 1 ; j <= faure.hisum ; j++ )
42 {
43 for ( int i = j+1 ; i <= faure.hisum ; i++ )
44 {
45 faure.coef [ i ] [ j ] = ( faure.coef [ i-1 ] [ j ] +faure.coef [ i-1 ] [ j-1 ] ) % faure.qs;
46 }
47 }
48 faure.nextn=faure.testn-1;
49 faure.hisum=3;
50 faure.rqs=1.0/faure.qs;
51 }
52
53 public static void gofaur ( double [ ] quasi )
54 {
55 int [ ] ytemp = new int [ 20 ];
56 int ztemp, ktemp,ltemp,mtemp;
57 double r;
58 ktemp=faure.testn;
59 ltemp=faure.nextn;
60 for ( int i = faure.hisum ; i >= 0 ; i-- )
61 {
62 ktemp=ktemp/faure.qs;
63 mtemp=ltemp % ktemp;
64 ytemp [ i ] = ( ltemp-mtemp ) /ktemp;
65 ltemp=mtemp;
66 }
67 r=ytemp [ faure.hisum ];
68 for ( int i = faure.hisum-1 ; i >= 0 ; i-- )
69 {
70 r=ytemp [ i ] +faure.rqs*r;
71 }
72 quasi [ 1-1 ] =r*faure.rqs;
73 for ( int k = 2 ; k <= faure.s ; k++ )
74 {
75 quasi [ k-1 ] =0;
76 r=faure.rqs;
77 for ( int j = 0 ; j <= faure.hisum ; j++ )
78 {
79 ztemp=0;
80 for ( int i = j ; i <= faure.hisum ; i++ )
81 {
82 ztemp=ztemp+faure.coef [ i ] [ j ] *ytemp [ i ];
83 }
84 ytemp [ j ] =ztemp % faure.qs;
85 quasi [ k-1 ] =quasi [ k-1 ] +ytemp [ j ] *r;
86 r=r*faure.rqs;
87 }
88 }
89 faure.nextn=faure.nextn+1;
90 if ( faure.nextn == faure.testn )
91 {
92 faure.testn=faure.testn*faure.qs;
93 faure.hisum=faure.hisum+1;
94 }
95 }
96
97 public static void inhalt ( boolean [ ] flag, int dimen, int atmost, double tiny, double [ ] quasi )
98 {
99 double delta;
100 halton.s=dimen;
101 flag [ 1-1 ] = ( halton.s > 1 && halton.s < 41 );
102 if ( ! flag [ 1-1 ] ) return;
103 halton.e=0.9* ( 1.0/ ( atmost*halton.prime [ halton.s-1 ] ) -10*tiny );
104 delta=100*tiny* ( atmost+1 ) * (Math.log( atmost )/Math.log(10));
105 flag [ 2-1 ] = ( delta <= 0.09* ( halton.e-10*tiny ) );
106 if ( ! flag [ 2-1 ] ) return;
107 for ( int i = 1 ; i <= halton.s ; i++ )
108 {
109 halton.prime [ i-1 ] =1/halton.prime [ i-1 ];
110 quasi [ i-1 ] =halton.prime [ i-1 ];
111 }
112 }
113
114 public static void gohalt ( double [ ] quasi )
115 {
116 double t, f,g,h;
117 for ( int i = 1 ; i <= halton.s ; i++ )
118 {
119 t=halton.prime [ i-1 ];
120 f=1-quasi [ i-1 ];
121 g=1;
122 h=t;
123 while ( f-h < halton.e )
124 {
125 g=h;
126 h=h*t;
127 }
128 quasi [ i-1 ] =g+h-f;
129 }
130 }
131
132 public static void main( String[] args )
133 {
134 boolean omit_output = false;
135 boolean[] flag = new boolean[2];
136 int m = 3, n = 20;
137
138 for ( int i = 0; i < args.length; i++ )
139 {
140 switch (args[i].charAt(1))
141 {
142 case 'm':
143 m = Integer.parseInt( args[++i] );
144 break;
145 case 'n':
146 n = Integer.parseInt( args[++i] );
147 break;
148 case 'o':
149 omit_output = true;
150 break;
151 }
152 }
153
154 System.err.println( "m: "+m+", n: "+n );
155 double[] quasi = new double[m];
156
157 infaur( flag, m, n );
158 for ( int i = 0; i < n; i++ )
159 {
160 gofaur( quasi );
161 if ( !omit_output )
162 {
163 for ( int j = 0; j < quasi.length; j++ )
164 System.out.print( quasi[j]+" " );
165 System.out.println("");
166 }
167 }
168 }
169 }
170
171 class halton
172 {
173 static int s;
174 static double e;
175 static double [ ] prime =
176 {
177 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
178 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107,
179 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173
180 };
181 }
182
183 class faure
184 {
185 static int s, qs, nextn, testn, hisum;
186 static double rqs;
187 static int [ ] [ ] coef = new int [ 20 ] [ 20 ];
188 }