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
|
C/***************************************************************************/
C/* ____Demonstrates SPRNG use for Monte Carlo integration____ */
C/* Compute pi using Monte Carlo integration. Random points are generated */
C/* in a square of size 2. The value of pi is estimated as four times */
C/* the proportion of samples that fall within a circle of unit radius. */
C/***************************************************************************/
program pif_simple
implicit none
#define SIMPLE_SPRNG ! simple interface
#include "sprng_f.h"
integer in, n, in_old, n_old, count_in_circle
real*8 pi, error, stderror, p, EXACT_PI
character filename*7
C--- reading in a generator type
integer gtype
#include "genf_types_menu.h"
print *,'Type in a generator type (integers: 0,1,2,3,4,5): '
read *, gtype
C---
EXACT_PI = 3.141592653589793238462643
p = EXACT_PI/4.0
C-- initialization: to initialize n, in_old, n_old, and filename
call initialize_function(gtype, n, in_old, n_old, filename)
in = count_in_circle(n) ! count samples in circle
in = in + in_old ! # in circle, in all runs
n = n + n_old ! # of samples, in all runs
pi = (4.0*in)/n
error = abs(pi - EXACT_PI)
stderror = 4*sqrt(p*(1.0-p)/n)
write(*,114) pi, n
write(*, 115) error, stderror
114 format('pi is estimated as ', f18.16, ' from ', i12, ' samples.')
115 format('Error = ', g18.8, ' standard error = ', g18.8)
call save_state(filename, in, n) !check-point final state
end
C-- count # of samples in cirlce
integer function count_in_circle(n)
implicit none
#include "sprng_f.h"
integer n, i, in
real*8 x, y
in = 0
do 10 i=1, n
x = 2*sprng() - 1.0 ! x coordinate
y = 2*sprng() - 1.0 ! y coordinate
if (x*x + y*y .lt. 1.0) then !check if point (x,y) is in circle
in = in +1
endif
10 continue
count_in_circle = in
return
end
C-- initialization --
subroutine initialize_function(gtype, n, in_old, n_old, filename)
implicit none
#include "sprng_f.h"
C---
integer gtype
C---
integer n, in_old, n_old, seed, size, junk
SPRNG_POINTER junkPtr
character filename*7, firstChar*1, buffer(MAX_PACKED_LENGTH)
print*, 'Enter 9 for a new run, or 2 to continue an old run'
read(*,113) firstChar
print*, 'Enter file name(length 7) to store final state:'
read(*, 111) filename ! 7 characters please
print*, 'Enter number of new samples:'
read(*,*)n
111 format(A7)
113 format(A1)
if (firstChar .eq. '9') then ! new set of runs
seed = make_sprng_seed() ! make seed from date/time information
junkPtr = init_sprng(gtype, seed, CRAYLCG) ! initialize
! stream
junk = print_sprng()
in_old = 0
n_old = 0
else ! continue from previously stored state
open(31, file = filename, status = 'unknown',
& form = 'unformatted')
read(31) in_old, n_old, size, buffer !read previous run info.
junkPtr = unpack_sprng(buffer)
close (31)
endif
return
end
C-- save the final state of the default stream into a file
subroutine save_state(filename, in, n)
implicit none
#include "sprng_f.h"
integer in, n, size
character filename*7, buffer(MAX_PACKED_LENGTH)
open(31, file = filename, status = 'unknown',
& form = 'unformatted')
size = pack_sprng(buffer)
write(31)in,n,size,buffer !write the current run info. for future use
close(31)
return
end
|
