spb/metrika Public

1//2//  PhiloxStream.swift3//  Metrika4//5//  Author:  Simon-Pierre Boucher6//  Contact: contact@spboucher.ai7//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8//910import Foundation1112/// Sequential random-variate stream over the Philox counter space, for13/// samplers whose draws are inherently ordered (MCMC chains). Unlike the14/// bootstrap's counter-addressable draws, a stream consumes positions one15/// by one — rejection samplers use a variable number — but the sequence16/// is fully determined by (seed, stream id), so chains are reproducible17/// and independent chains never overlap.18///19/// Stream ids occupy bit 62 of the word-counter space plus a 2⁴⁸-word20/// block per id, disjoint from bootstrap draws (low positions) and21/// permutation keys (bit 63).22public struct PhiloxStream: Sendable {23    private let generator: Philox4x3224    private let base: UInt6425    private var position: UInt64 = 02627    public init(seed: UInt64, stream: UInt64 = 0) {28        self.generator = Philox4x32(seed: seed)29        self.base = (UInt64(1) << 62) &+ stream &* (UInt64(1) << 48)30    }3132    /// Uniform in (0, 1) — endpoints excluded so inverse-CDF transforms33    /// and logs stay finite.34    public mutating func nextUniform() -> Double {35        position &+= 136        let u = generator.uniform(at: base &+ position)37        return min(max(u, 5e-324), 1 - 2.2e-16)38    }3940    /// Standard normal via Box–Muller (two uniforms per pair, the spare41    /// is cached).42    private var cachedNormal: Double?43    public mutating func nextNormal() -> Double {44        if let cached = cachedNormal {45            cachedNormal = nil46            return cached47        }48        let u1 = nextUniform()49        let u2 = nextUniform()50        let radius = (-2 * Foundation.log(u1)).squareRoot()51        let angle = 2 * Double.pi * u252        cachedNormal = radius * Foundation.sin(angle)53        return radius * Foundation.cos(angle)54    }5556    /// Gamma(shape, rate) via Marsaglia–Tsang squeeze (with the standard57    /// boost for shape < 1).58    public mutating func nextGamma(shape: Double, rate: Double) -> Double {59        precondition(shape > 0 && rate > 0)60        if shape < 1 {61            // Gamma(a) = Gamma(a+1) · U^(1/a)62            let boosted = nextGamma(shape: shape + 1, rate: rate)63            return boosted * Foundation.pow(nextUniform(), 1 / shape)64        }65        let d = shape - 1.0 / 3.066        let c = 1 / (9 * d).squareRoot()67        while true {68            let z = nextNormal()69            let v = (1 + c * z) * (1 + c * z) * (1 + c * z)70            guard v > 0 else { continue }71            let u = nextUniform()72            if Foundation.log(u) < 0.5 * z * z + d - d * v + d * Foundation.log(v) {73                return d * v / rate74            }75        }76    }77}78