Monte Carlo methods: Difference between revisions

Line 669:

include futlib.numeric

default(f32)

fun dirvcts(): [2][30]~~int~~ =

fun dirvcts(): [2][30]i32 =

[

Line 682:

Line 684:

fun grayCode(x: ~~int~~): ~~int~~ = (x >> 1) ^ x

fun grayCode(x: i32): i32 = (x >> 1) ^ x

----------------------------------------

--- Sobol Generator

----------------------------------------

fun testBit(n: ~~int~~, ind: ~~int~~): bool =

fun testBit(n: i32, ind: i32): bool =

let t = (1 << ind) in (n & t) == t

fun xorInds(n: ~~int~~) (dir_vs: [num_bits]~~int~~): ~~int~~ =

fun xorInds(n: i32) (dir_vs: [num_bits]i32): i32 =

let reldv_vals = zipWith (fn dv i =>

let reldv_vals = zipWith (\ dv i ->

if testBit(grayCode n,i)

then dv else 0)

Line 697:

Line 699:

in reduce (^) 0 reldv_vals

fun sobolIndI (dir_vs: [m][num_bits]~~int~~, n: ~~int~~): [m]~~int~~ =

fun sobolIndI (dir_vs: [m][num_bits]i32, n: i32): [m]i32 =

map (xorInds n) dir_vs

fun sobolIndR(dir_vs: [m][num_bits]~~int~~) (n: ~~int~~ ): [m]f32 =

fun sobolIndR(dir_vs: [m][num_bits]i32) (n: i32 ): [m]f32 =

let divisor = 2.0 ** f32(num_bits)

let arri = sobolIndI( dir_vs, n )

in map (fn (x: ~~int~~): f32 => f32(x) / divisor) arri

in map (\ (x: i32): f32 -> f32(x) / divisor) arri

fun main(n: ~~int~~): f32 =

fun main(n: i32): f32 =

let rand_nums = map (sobolIndR (dirvcts())) (iota n)

let dists = map (fn xy =>

let dists = map (\xy ->

let (x,y) = (xy[0],xy[1]) in ~~sqrt32~~(x*x + y*y))

let (x,y) = (xy[0],xy[1]) in F32.sqrt(x*x + y*y))

rand_nums

let bs = map (fn d => if d <= 1.0f32 then 1 else 0) dists

let bs = map (\d -> if d <= 1.0f32 then 1 else 0) dists

let inside = reduce (+) 0 bs