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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
|
/***************************************************************************/
/* ____Demonstrates SPRNG use for Monte Carlo integration____ */
/* Compute pi using Monte Carlo integration. Random points are generated */
/* in a square of size 2. The value of pi is estimated as four times the */
/* the proportion of samples that fall within a circle of unit radius. */
/***************************************************************************/
#include <stdio.h>
#include <math.h>
#include <string.h>
#define SIMPLE_SPRNG /* simple interface */
#include "sprng.h"
#define EXACT_PI 3.141592653589793238462643
int gtype; /*--- */
main(argc,argv)
int argc;
char *argv[];
{
int in, n, in_old, n_old;
double pi, error, stderror, p=EXACT_PI/4.0;
char filename[80];
/*--- reading in a generator type */
#include "gen_types_menu.h"
printf("Type in a generator type (integers: 0,1,2,3,4,5): ");
scanf("%d", >ype);
initialize_function(&n, &in_old, &n_old, filename); /* read args & initialize */
in = count_in_circle(n); /* count samples in circle */
in += in_old; /* # in circle, in all runs */
n += n_old; /* # of samples, in all runs */
pi = (4.0*in)/n; /* estimate pi */
error = fabs(pi - EXACT_PI); /* determine error */
stderror = 4*sqrt(p*(1.0-p)/n); /* standard error */
printf( "pi is estimated as %18.16f from %12d samples.\n", pi, n);
printf( "\tError = %10.8g, standard error = %10.8g\n", error, stderror);
save_state(filename, in, n); /* check-point final state */
}
int count_in_circle(n) /* count # of samples in circle */
int n;
{
int i, in;
double x, y;
for (i=in=0; i<n; i++) /* find number of samples in circle */
{
x = 2*sprng() - 1.0; /* x coordinate */
y = 2*sprng() - 1.0; /* y coordinate */
if (x*x + y*y < 1.0) /* check if point (x,y) is in circle */
in++;
}
return in;
}
/* Read arguments and initialize stream */
int initialize_function(n, in_old, n_old, filename)
int *n, *in_old, *n_old;
char *filename;
{
int seed, size, new_old;
char buffer[MAX_PACKED_LENGTH];
FILE *fp;
printf("Enter 9 for a new run, or 2 for the continuation of an old run:\n");
scanf("%d", &new_old);
printf("Enter name of file to store/load state of the stream:\n");
scanf("%s", filename);
printf("Enter number of new samples:\n");
scanf("%d", n);
if(new_old == 9) /* new set of runs */
{
seed = make_sprng_seed(); /* make seed from date/time information */
init_sprng(gtype,seed,CRAYLCG); /* initialize stream */
print_sprng(); /* print information abour stream */
*in_old = 0;
*n_old = 0;
}
else /* continue from previously stored state */
{
fp = fopen(filename,"r"); /* open file */
if(fp == NULL)
{
fprintf(stderr,"ERROR opening file %s\n", filename);
exit(1);
}
fread(in_old,1,sizeof(int),fp); /* cumulative # in circle previously */
fread(n_old,1,sizeof(int),fp); /* cumulative # of samples previously */
fread(&size,1,sizeof(int),fp); /* size of stored stream state */
fread(buffer,1,size,fp); /* copy stream state to buffer */
unpack_sprng(buffer); /* retrieve state of the stream */
fclose(fp); /* close file */
}
return 0;
}
int save_state(filename, in, n) /* store the state */
char *filename;
int in, n;
{
FILE *fp;
char *bytes;
int size;
fp = fopen(filename,"w"); /* open file to store state */
if(fp == NULL)
{
fprintf(stderr,"Could not open file %s for writing\nCheck path or permissions\n", filename);
exit(1);
}
fwrite(&in,1,sizeof(int),fp); /* store # in circle in all runs */
fwrite(&n,1,sizeof(int),fp); /* store # of samples in all runs */
size = pack_sprng(&bytes); /* pack stream state into an array */
fwrite(&size,1,sizeof(int),fp); /* store # of bytes required for storage */
fwrite(bytes,1,size,fp); /* store stream state */
fclose(fp);
free(bytes); /* free memory needed to store stream state*/
return 0;
}
|
