1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
/* A two-dimensional random walk test based on the code of I. Vattulainen,
et al */
#include <stdio.h>
#include "tests.h"
#include <math.h>
#define sgn(a) ((a<0)?0:1)
#ifndef ANSI_ARGS
#ifdef __STDC__
#define ANSI_ARGS(args) args
#else
#define ANSI_ARGS(args) ()
#endif
#endif
void random_walk ANSI_ARGS((int walk_length, double *quadrant));
#ifdef __STDC__
main(int argc, char *argv[])
#else
main(argc, argv)
int argc;
char *argv[];
#endif
{
long ntests, walk_length, i;
double result, expected, chisum, quadrant[4];
quadrant[0] = quadrant[1] = quadrant[2] = quadrant[3] = 0;
if(argc != N_STREAM_PARAM + 2)
{
fprintf(stderr,"USAGE: %s (... %d arguments)\n",argv[0], N_STREAM_PARAM+1);
exit(1);
}
ntests = init_tests(argc,argv);
walk_length = atoi(argv[N_STREAM_PARAM+1]);
for(i=0; i<ntests; i++)
{
random_walk(walk_length, quadrant);
next_stream();
}
#if defined(SPRNG_MPI)
reduce_sum(quadrant,4);
#endif
if(proc_rank == 0)
{
chisum = 0.0;
expected = NTESTS/4.0;
for(i=0; i<4; i++)
{
chisum += (quadrant[i]-expected)*(quadrant[i]-expected)/expected;
}
result = chipercent(chisum, 3);
printf("Chi square value = %f, percent = %f\n", chisum, result*100);
}
#if defined(SPRNG_MPI)
MPI_Finalize();
#endif
}
#ifdef __STDC__
void random_walk(int walk_length, double *quadrant)
#else
void random_walk(walk_length, quadrant)
int walk_length;
double *quadrant;
#endif
{
int x, y, i, j;
double rn;
x = y = 0;
for(i=0; i<walk_length; i++)
{
rn = get_rn();
if(rn > 0.75)
x++;
else if(rn > 0.5)
x--;
else if(rn > 0.25)
y++;
else
y--;
}
if(x == 0 && y == 0)
{
quadrant[0] += 0.25;
quadrant[1] += 0.25;
quadrant[2] += 0.25;
quadrant[3] += 0.25;
}
else if(x == 0)
{
quadrant[sgn(y)] += 0.5;
quadrant[2+sgn(y)] += 0.5;
}
else if(y == 0)
{
quadrant[2*sgn(x)] += 0.5;
quadrant[2*sgn(x)+1] += 0.5;
}
else
quadrant[2*sgn(x)+sgn(y)] += 1.0;
}
|
