Epidemiological simulation in Zig

sir.zig

// zig build-exe sir.zig --release-fast --strip
const std = @import("std");

inline fn randbn(n: i64, p: f64, rng: *std.rand.Xoroshiro128) i64 {
    const log_q: f64 = std.math.ln(1.000000 - p);
    var x: i64 = i64(0);
    var sum: f64 = 0.000000;
    while (1 != 0) {
        sum += std.math.ln(rng.random.float(f64))/@intToFloat(f64, n - x);
        if (sum < log_q) break;
        x += i64(1);
    }
    return x;
}

inline fn sir(u: *[4]i64, du: *[4]i64, parms: [5]f64, rng: *std.rand.Xoroshiro128) void {
    var S: i64 = u[0];
    var I: i64 = u[1];
    var R: i64 = u[2];
    var Y: i64 = u[3];
    var beta: f64 = parms[0];
    var gamma_0: f64 = parms[1];
    var iota: f64 = parms[2];
    var N: f64 = parms[3];
    var dt: f64 = parms[4];
    var lambd: f64 = beta * (@intToFloat(f64,I) + iota)/N;
    var ifrac: f64 = 1.000000 - std.math.exp((-lambd) * dt);
    var rfrac: f64 = 1.000000 - std.math.exp((-gamma_0) * dt);
    var infection: i64 = randbn(S, ifrac, rng);
    var recovery: i64 = randbn(I, rfrac, rng);
    du[0] = (S - infection);
    du[1] = (I + infection - recovery);
    du[2] = (R + recovery);
    du[3] = (Y + infection);
    return;
}

fn simulate() f64{
    const parms = [_]f64{2.62110617498984, 0.5384615384615384, 0.5, 403.0, 0.1};
    const tf: i64 = 540;
    const nsims: i64 = 1000;
    var numerator: i64 = 0;
    var seed: u64 = 123456;
    var rng = std.rand.Xoroshiro128.init(seed);
    var i: i64 = 0;
    while (i<nsims){
        var u = [_]i64{60,1,342,0};
        var du = [_]i64{0,0,0,0};
        var j: i64 = 0;
        while (j<tf){
            sir(&u,&du,parms,&rng);
            u = du;
            j += 1;
        }
        numerator += u[3];
        i += 1;
    }
    return @intToFloat(f64,numerator)/@intToFloat(f64,nsims);
}

pub fn main() void {
    var m: f64 = simulate();
    std.debug.warn("{}\n", m);
}