OpenEuphoria

Re: Phix Random Numbers

• Posted by andreasWagner 1 month ago
• 544 views

Hallo,

I just delegated that to Gemini for you. The period should be enough for at least 100k in value.
It came right back like that, just a few minor syntax errors related to "constant".

-- Custom LCG Pseudorandom Number Generator for Monte Carlo simulations 
 
-- Global variable for the current generator state (the seed) 
atom lcg_state = 123456789  -- Any starting seed (!= 0) 
 
-- Parameters according to Park & Miller "Minimal Standard" 
-- The period is 2^31 - 1 (over 2 billion values), ideal for 100k runs 
constant lcg_m = 2147483647  -- 2^31 - 1 (Mersenne prime) 
constant lcg_a = 48271     -- Multiplier 
constant lcg_q = 44488     -- m / a (used in Schrage's method to prevent overflow) 
constant lcg_r = 3399      -- m % a 
 
-- Function to generate the next pseudorandom float between 0.0 and 1.0 
function lcg_next_float() 
    -- Schrage's method to avoid 64-bit overflows during multiplication 
    atom hi = floor(lcg_state / lcg_q) 
    atom lo = remainder(lcg_state, lcg_q) 
     
    lcg_state = lcg_a * lo - lcg_r * hi 
     
    if lcg_state <= 0 then 
        lcg_state += lcg_m 
    end if 
     
    -- Normalize to the range [0.0, 1.0) 
    return lcg_state / lcg_m 
end function 
 
-- === MONTE CARLO EXPERIMENT: Estimating Pi === 
-- We drop "points" into a square and count how many land inside the quarter circle. 
 
integer n_hits = 0 
integer n_total = 10000000 
 
for i = 1 to n_total do 
    atom x = lcg_next_float() 
    atom y = lcg_next_float() 
     
    -- Pythagorean theorem: Check distance to the origin (x*x + y*y <= 1) 
    if (x * x + y * y) <= 1.0 then 
        n_hits += 1 
    end if 
end for 
 
-- Pi is calculated as: 4 * (hits / total) 
atom pi_estimate = 4.0 * (n_hits / n_total) 
 
printf(1, "Monte Carlo simulation with %d values:\n", {n_total}) 
printf(1, "Estimated value for Pi = %1.5f\n", {pi_estimate}) 
wait_key()