Contents of file 'wiener/wiener.c':
1 /*10:*/
2 #line 628 "./wiener.w"
3
4 /*11:*/
5 #line 653 "./wiener.w"
6
7 #include <math.h>
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include <ctype.h>
12
13 /*:11*/
14 #line 629 "./wiener.w"
15
16 /*12:*/
17 #line 679 "./wiener.w"
18
19 #define VERSION "1.0"
20 #define COPYRIGHT "Copyright (C) 2011, Fredrik Jonsson <http://jonsson.eu>"
21 #define SUCCESS (0)
22 #define FAILURE (1)
23 #define QUALITY (1009)
24 #define KK (100)
25 #define LL (37)
26 #define DEFAULT_SEED (310952)
27 #define cmd_match(s,l,c) ((!strcmp((c),(s)))||(!strcmp((c),(l))))
28 #define MODE_WIENER_PROCESS (0)
29 #define MODE_LOCKED_UNIFORM_DISTRIBUTION (1)
30 #define MODE_LOCKED_NORMAL_DISTRIBUTION (2)
31
32 /*:12*/
33 #line 630 "./wiener.w"
34
35 /*13:*/
36 #line 701 "./wiener.w"
37
38 extern char*optarg;
39 char*progname;
40 double ranf_arr_buf[QUALITY];
41 double ranf_arr_dummy= -1.0,ranf_arr_started= -1.0;
42 double*ranf_arr_ptr= &ranf_arr_dummy;
43
44 /*:13*/
45 #line 631 "./wiener.w"
46
47 /*14:*/
48 #line 712 "./wiener.w"
49
50 /*15:*/
51 #line 760 "./wiener.w"
52
53 #define mod_sum(x,y) (((x)+(y))-(int)((x)+(y)))
54 double ran_u[KK];
55
56 void ranf_array(double aa[],int n)
57 {
58 register int i,j;
59
60 for(j= 0;j<KK;j++)aa[j]= ran_u[j];
61 for(;j<n;j++)aa[j]= mod_sum(aa[j-KK],aa[j-LL]);
62 for(i= 0;i<LL;i++,j++)ran_u[i]= mod_sum(aa[j-KK],aa[j-LL]);
63 for(;i<KK;i++,j++)ran_u[i]= mod_sum(aa[j-KK],ran_u[i-LL]);
64 }
65
66 void ranf_matrix(double**aa,int mm,int dd)
67 {
68 register int i,j,col;
69
70 for(col= 0;col<dd;col++){
71 for(j= 0;j<KK;j++)aa[j][col]= ran_u[j];
72 for(;j<mm;j++)aa[j][col]= mod_sum(aa[j-KK][col],aa[j-LL][col]);
73 for(i= 0;i<LL;i++,j++)ran_u[i]= mod_sum(aa[j-KK][col],aa[j-LL][col]);
74 for(;i<KK;i++,j++)ran_u[i]= mod_sum(aa[j-KK][col],ran_u[i-LL]);
75 }
76 }
77
78 /*:15*/
79 #line 713 "./wiener.w"
80
81 /*16:*/
82 #line 798 "./wiener.w"
83
84 #define TT 70
85 #define is_odd(s) ((s)&1)
86
87 void ranf_start(long seed){
88 register int t,s,j;
89 double u[KK+KK-1];
90 double ulp= (1.0/(1L<<30))/(1L<<22);
91 double ss= 2.0*ulp*((seed&0x3fffffff)+2);
92
93 for(j= 0;j<KK;j++){
94 u[j]= ss;
95 ss+= ss;
96 if(ss>=1.0)ss-= 1.0-2*ulp;
97 }
98 u[1]+= ulp;
99 for(s= seed&0x3fffffff,t= TT-1;t;){
100 for(j= KK-1;j> 0;j--)
101 u[j+j]= u[j],u[j+j-1]= 0.0;
102 for(j= KK+KK-2;j>=KK;j--){
103 u[j-(KK-LL)]= mod_sum(u[j-(KK-LL)],u[j]);
104 u[j-KK]= mod_sum(u[j-KK],u[j]);
105 }
106 if(is_odd(s)){
107 for(j= KK;j> 0;j--)u[j]= u[j-1];
108 u[0]= u[KK];
109 u[LL]= mod_sum(u[LL],u[KK]);
110 }
111 if(s)s>>= 1;else t--;
112 }
113 for(j= 0;j<LL;j++)ran_u[j+KK-LL]= u[j];
114 for(;j<KK;j++)ran_u[j-LL]= u[j];
115 for(j= 0;j<10;j++)ranf_array(u,KK+KK-1);
116 ranf_arr_ptr= &ranf_arr_started;
117 }
118
119 #define ranf_arr_next() (*ranf_arr_ptr>=0? *ranf_arr_ptr++: ranf_arr_cycle())
120 double ranf_arr_cycle()
121 {
122 if(ranf_arr_ptr==&ranf_arr_dummy)
123 ranf_start(314159L);
124 ranf_array(ranf_arr_buf,QUALITY);
125 ranf_arr_buf[KK]= -1;
126 ranf_arr_ptr= ranf_arr_buf+1;
127 return ranf_arr_buf[0];
128 }
129
130 /*:16*/
131 #line 714 "./wiener.w"
132
133 /*17:*/
134 #line 874 "./wiener.w"
135
136 #define PI (3.14159265358979323846264338)
137 void normdist(double**aa,int mm,int dd){
138 register int j,k;
139 register double f,z;
140
141 for(j= 0;j<dd;j++){
142 for(k= 0;k<mm-1;k+= 2){
143 if(aa[k][j]> 0.0){
144 f= sqrt(-2*log(aa[k][j]));
145 z= 2.0*PI*aa[k+1][j];
146 aa[k][j]= f*cos(z);
147 aa[k+1][j]= f*sin(z);
148 }else{
149 fprintf(stderr,"%s: Error: Zero element detected!\n",progname);
150 fprintf(stderr,"%s: (row %d, column %d)\n",progname,k,j);
151 }
152 }
153 }
154 return;
155 }
156 #undef PI
157
158 /*:17*/
159 #line 715 "./wiener.w"
160
161 /*18:*/
162 #line 942 "./wiener.w"
163
164 void wiener(double**aa,int mm,int dd,int seed,short mode)
165 {
166 register int j,k;
167
168 ranf_start(seed);
169 ranf_matrix(aa,mm,dd);
170 if(mode==MODE_LOCKED_UNIFORM_DISTRIBUTION)return;
171 normdist(aa,mm,dd);
172 if(mode==MODE_LOCKED_NORMAL_DISTRIBUTION)return;
173 for(j= 0;j<dd;j++){
174 aa[0][j]= 0.0;
175 for(k= 1;k<mm;k++){
176 aa[k][j]+= aa[k-1][j];
177 }
178 }
179 }
180
181 /*:18*/
182 #line 716 "./wiener.w"
183
184 /*19:*/
185 #line 965 "./wiener.w"
186
187 double**dmatrix(long nrl,long nrh,long ncl,long nch)
188 {
189 long i,nrow= nrh-nrl+1,ncol= nch-ncl+1;
190 double**m;
191 m= (double**)malloc((size_t)((nrow+1)*sizeof(double*)));
192 if(!m){
193 fprintf(stderr,"%s: Allocation failure 1 in dmatrix()\n",progname);
194 exit(FAILURE);
195 }
196 m+= 1;
197 m-= nrl;
198 m[nrl]= (double*)malloc((size_t)((nrow*ncol+1)*sizeof(double)));
199 if(!m[nrl]){
200 fprintf(stderr,"%s: Allocation failure 2 in dmatrix()\n",progname);
201 exit(FAILURE);
202 }
203 m[nrl]+= 1;
204 m[nrl]-= ncl;
205 for(i= nrl+1;i<=nrh;i++)m[i]= m[i-1]+ncol;
206 return m;
207 }
208
209 /*:19*/
210 #line 717 "./wiener.w"
211
212 /*20:*/
213 #line 992 "./wiener.w"
214
215 void free_dmatrix(double**m,long nrl,long nrh,long ncl,long nch){
216 free((char*)(m[nrl]+ncl-1));
217 free((char*)(m+nrl-1));
218 }
219
220 /*:20*/
221 #line 718 "./wiener.w"
222
223 /*21:*/
224 #line 1000 "./wiener.w"
225
226 void display_help_message(void){
227 fprintf(stderr,"Usage: %s M [options]\n",progname);
228 fprintf(stderr,"Options:\n");
229 fprintf(stderr," -h, --help\n");
230 fprintf(stderr," Display this help message and exit clean.\n");
231 fprintf(stderr," -v, --verbose\n");
232 fprintf(stderr," Toggle verbose mode. Default: off.\n");
233 fprintf(stderr," -M, --num_samples <M>\n");
234 fprintf(stderr," Generate M samples of data. Here M should always be\n");
235 fprintf(stderr," an even number, greater than the long lag KK=%d.\n",KK);
236 fprintf(stderr," If an odd number is specified, the program will\n");
237 fprintf(stderr," automatically adjust this to the next higher\n");
238 fprintf(stderr," integer. Default: M=%d.\n",KK);
239 fprintf(stderr," -D, --dimension <D>\n");
240 fprintf(stderr," Generate D-dimensional samples of data. Default: D=1.\n");
241 fprintf(stderr," -s, --seed <seed>\n");
242 fprintf(stderr," Define a custom seed number for the initialization\n");
243 fprintf(stderr," of the random number generator. Default: seed=%d.\n",
244 DEFAULT_SEED);
245 fprintf(stderr," -u, --uniform\n");
246 fprintf(stderr," Instead of generating a sequence of D-dimensional\n");
247 fprintf(stderr," corresponding to a Wiener process, lock the program\n");
248 fprintf(stderr," to simply generate a uniform distribution of\n");
249 fprintf(stderr," D-dimensional points, with each element distributed\n");
250 fprintf(stderr," over the interval [0,1].\n");
251 fprintf(stderr," -n, --normal\n");
252 fprintf(stderr," Instead of generating a sequence of D-dimensional\n");
253 fprintf(stderr," corresponding to a Wiener process, lock the program\n");
254 fprintf(stderr," to simply generate a normal distribution of\n");
255 fprintf(stderr," D-dimensional points, with each element distributed\n");
256 fprintf(stderr," with zero expectation value and unit variance.\n");
257 }
258
259 /*:21*/
260 #line 719 "./wiener.w"
261
262 /*22:*/
263 #line 1039 "./wiener.w"
264
265 short pathcharacter(int ch){
266 return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
267 (ch=='-')||(ch=='+'));
268 }
269
270 /*:22*/
271 #line 720 "./wiener.w"
272
273 /*23:*/
274 #line 1051 "./wiener.w"
275
276 char*strip_away_path(char filename[]){
277 int j,k= 0;
278 while(pathcharacter(filename[k]))k++;
279 j= (--k);
280 while(isalnum((int)(filename[j])))j--;
281 j++;
282 return(&filename[j]);
283 }
284
285 /*:23*/
286 #line 721 "./wiener.w"
287
288
289 /*:14*/
290 #line 632 "./wiener.w"
291
292
293 int main(int argc,char*argv[])
294 {
295 /*24:*/
296 #line 1092 "./wiener.w"
297
298 double**aa;
299 unsigned long mm= KK;
300 unsigned dd= 1;
301 int seed= DEFAULT_SEED;
302 short mode= MODE_WIENER_PROCESS,verbose= 0;
303 int no_arg;
304 register int j,k;
305
306 /*:24*/
307 #line 636 "./wiener.w"
308
309 /*25:*/
310 #line 1106 "./wiener.w"
311
312 progname= strip_away_path(argv[0]);
313 no_arg= argc;
314 while(--argc){
315 if(cmd_match("-h","--help",argv[no_arg-argc])){
316 display_help_message();
317 exit(SUCCESS);
318 }else if(cmd_match("-v","--verbose",argv[no_arg-argc])){
319 verbose= (verbose?0:1);
320 }else if(cmd_match("-M","--num_samples",argv[no_arg-argc])){
321 --argc;
322 if(!sscanf(argv[no_arg-argc],"%lu",&mm)){
323 fprintf(stderr,"%s: Error detected when parsing the number of "
324 "samples.\n",progname);
325 display_help_message();
326 exit(FAILURE);
327 }
328 if(mm<KK){
329 fprintf(stderr,"%s: The M number of data points must be at least "
330 "the long lag of the generator, M >= KK = %d.\n",progname,KK);
331 exit(FAILURE);
332 }
333 if(is_odd(mm))mm++;
334 }else if(cmd_match("-D","--dimension",argv[no_arg-argc])){
335 --argc;
336 if(!sscanf(argv[no_arg-argc],"%ud",&dd)){
337 fprintf(stderr,"%s: Error detected when parsing dimension.\n",
338 progname);
339 display_help_message();
340 exit(FAILURE);
341 }
342 if(dd<1){
343 fprintf(stderr,"%s: Dimension D should be at least 1.\n",progname);
344 exit(FAILURE);
345 }
346 }else if(cmd_match("-s","--seed",argv[no_arg-argc])){
347 --argc;
348 if(!sscanf(argv[no_arg-argc],"%d",&seed)){
349 fprintf(stderr,"%s: Error detected when parsing the seed of the "
350 "initializer.\n",progname);
351 display_help_message();
352 exit(FAILURE);
353 }
354 }else if(cmd_match("-u","--uniform",argv[no_arg-argc])){
355 mode= MODE_LOCKED_UNIFORM_DISTRIBUTION;
356 }else if(cmd_match("-n","--normal",argv[no_arg-argc])){
357 mode= MODE_LOCKED_NORMAL_DISTRIBUTION;
358 }else{
359 fprintf(stderr,"%s: Sorry, I do not recognize option '%s'.\n",
360 progname,argv[no_arg-argc]);
361 display_help_message();
362 exit(FAILURE);
363 }
364 }
365 if(verbose)
366 fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
367
368 /*:25*/
369 #line 637 "./wiener.w"
370
371 /*26:*/
372 #line 1169 "./wiener.w"
373
374 aa= dmatrix(0,mm-1,0,dd-1);
375
376 /*:26*/
377 #line 638 "./wiener.w"
378
379 /*27:*/
380 #line 1174 "./wiener.w"
381
382 wiener(aa,mm,dd,seed,mode);
383
384 /*:27*/
385 #line 639 "./wiener.w"
386
387 /*28:*/
388 #line 1179 "./wiener.w"
389
390 for(k= 0;k<mm;k++){
391 for(j= 0;j<dd-1;j++)
392 printf("%.20f ",aa[k][j]);
393 printf("%.20f\n",aa[k][j]);
394 }
395
396 /*:28*/
397 #line 640 "./wiener.w"
398
399 /*29:*/
400 #line 1188 "./wiener.w"
401
402 free_dmatrix(aa,0,mm-1,0,dd-1);
403
404 /*:29*/
405 #line 641 "./wiener.w"
406
407 return(SUCCESS);
408 }
409
410 /*:10*/
411
Generated by ::viewsrc::