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

Reply via email to