Hello,
I have tried the following code to calculate pi using Sobol.jl (
https://github.com/stevengj/Sobol.jl).
using Sobol
function sobol_test( n :: Int64 )
v = [ 0.0, 0.0 ]
sobo = SobolSeq( 2 )
count = 0.0
pi_mc = 0.0
for i = 1:n
next!( sobo, v )
if v[1]^2 + v[2]^2 < 1.0
count += 1.0
end
pi_mc = 4.0 * count / i
if i % 10^7 == 0
println( "i/10^7 = ", div( i, 10^6 ), " pi_mc=", pi_mc )
if abs( pi_mc - pi ) > 0.1
@show v
end
end
end
println( "pi(mc) = ", pi_mc )
println( "pi(exact) = ", pi )
end
sobol_test( 3 * 10^9 )
It converges to pi = 3.1415... very quickly, but once the number of samples
generated exceeds
2^32-1 = 4294967295, the code started to give very large values for v[1]
and v[2]:
...
i/10^7 = 2120 pi_mc=3.1415926773584903
i/10^7 = 2130 pi_mc=3.1415927136150237
i/10^7 = 2140 pi_mc=3.14159273271028
i/10^7 = 2150 pi_mc=3.1379158883720932
v = [4.8933891e7,1.503080789e9]
i/10^7 = 2160 pi_mc=3.1233885
v = [4.294403e6,1.936346197e9]
i/10^7 = 2170 pi_mc=3.108995004608295
v = [6.4634755e7,1.397200085e9]
...
(please note that a pair of Sobol numbers are used for each iteration, so
i/10^7 = 2140 roughly
corresponds to 2^32-1.) Although the webpage of Sobol.jl says that Mersenne
Twister will be used after the 2^32-1 samples,
do I need to switch manually to the built-in random number generator in
Julia? Or am I doing something wrong
about the treatment of Integer? I am using Linux 64bit and assuming that
integer literals is of Int64.
Thanks