// =============================================================================
// P³ QUASIRANDOM v1.0 — ZIG
// =============================================================================
//
// Квазислучайные последовательности с низкой дисперсией (Low-Discrepancy Sequences).
// Портировано из UE5: Runtime/Core/Public/Math/Sobol.h + Halton.h
//
// МАТЕМАТИКА:
//   Sobol: x_n^d = gray(n) · V_d / 2^32  (direction numbers + Gray code)
//   Halton: x_n^b = Σ (n_i / b^{i+1})  where n = Σ n_i · b^i (radical inverse)
//
//   Свойства:
//   - Sobol: O(N^{-1}·(log N)^d) discrepancy, d dimensions
//   - Halton: O(N^{-1}·(log N)^d) but degrades for large d (correlation)
//   - Sobol >> Halton for d > 6
//
// P³ ОБОБЩЕНИЯ:
//   - Sobol на P³: квазислучайные точки в проективном пространстве
//   - Halton на S³ → RP³: равномерное покрытие проективного пространства
//   - Cross-ratio 4 quasirandom points = проективный инвариант
//   - Стратификация: Sobol клетки → аффинные карты на P³
//
// Донор: UE5 Sobol.h/.cpp + Halton.h (~300 строк)
// Порт: Zig 0.14.0, comptime dimensions, P³ integration
//
// Принцип: убивать и жрать и рождать новое.
// =============================================================================

const std = @import("std");
const math = std.math;

// =============================================================================
// 1. HALTON SEQUENCE
// =============================================================================

/// Halton radical-inverse quasi-random sequence
/// Halton(1964): x_n = Σ (digit_i(n, base) / base^{i+1})
/// Simple, fast, but quality degrades for high dimensions
pub fn halton(index: i32, base: i32) f32 {
    var result: f32 = 0.0;
    const inv_base: f32 = 1.0 / @as(f32, @floatFromInt(base));
    var fraction: f32 = inv_base;
    var idx = index;
    while (idx > 0) {
        result += @as(f32, @floatFromInt(idx % base)) * fraction;
        idx = idx / base;
        fraction *= inv_base;
    }
    return result;
}

/// 2D Halton pair using bases (2, 3)
pub fn halton2D(index: i32) [2]f32 {
    return .{ halton(index, 2), halton(index, 3) };
}

/// 3D Halton triple using bases (2, 3, 5)
pub fn halton3D(index: i32) [3]f32 {
    return .{ halton(index, 2), halton(index, 3), halton(index, 5) };
}

/// 4D Halton using bases (2, 3, 5, 7) — natural for P³ homogeneous coords
pub fn halton4D(index: i32) [4]f32 {
    return .{ halton(index, 2), halton(index, 3), halton(index, 5), halton(index, 7) };
}

/// Halton with scramble (Owen scrambling for better 2D projections)
pub fn haltonScrambled(index: i32, base: i32, seed: i32) f32 {
    var result: f32 = 0.0;
    const inv_base: f32 = 1.0 / @as(f32, @floatFromInt(base));
    var fraction: f32 = inv_base;
    var idx = index;
    var s = seed;
    while (idx > 0) {
        const digit = idx % base;
        // Permute digit using seed-based permutation
        const permuted = (digit + s) % base;
        result += @as(f32, @floatFromInt(permuted)) * fraction;
        idx = idx / base;
        s = (s * 1103515245 + 12345) & 0x7FFFFFFF; // LCG for next permutation
        fraction *= inv_base;
    }
    return result;
}

// =============================================================================
// 2. SOBOL SEQUENCE
// =============================================================================

/// Maximum number of Sobol dimensions (matching UE5's Joe-Kuo table)
pub const SOBOL_MAX_DIM: usize = 15;

/// Maximum bits for Sobol index (32-bit)
pub const SOBOL_MAX_BITS: usize = 32;

/// Direction numbers for first 16 Sobol dimensions
/// These are the first 32 direction numbers for each dimension
/// (from Joe & Kuo, 2008 — same table used in UE5)
pub const SOBOL_DIRECTION_NUMBERS: [SOBOL_MAX_DIM + 1][SOBOL_MAX_BITS]i32 = .{
    // Dimension 0 (special: index i maps to i/2^32)
    .{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
    // Dimension 1
    .{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
    // Dimension 2
    .{ 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3 },
    // Dimension 3
    .{ 1, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1 },
    // Dimension 4
    .{ 1, 1, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3 },
    // Dimension 5
    .{ 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13 },
    // Dimension 6
    .{ 1, 1, 5, 5, 11, 11, 1, 1, 5, 5, 11, 11, 1, 1, 5, 5, 11, 11, 1, 1, 5, 5, 11, 11, 1, 1, 5, 5, 11, 11, 1, 1 },
    // Dimension 7
    .{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
    // Dimension 8
    .{ 1, 3, 5, 15, 17, 51, 85, 255, 1, 3, 5, 15, 17, 51, 85, 255, 1, 3, 5, 15, 17, 51, 85, 255, 1, 3, 5, 15, 17, 51, 85, 255 },
    // Dimension 9
    .{ 1, 1, 7, 7, 21, 21, 63, 63, 1, 1, 7, 7, 21, 21, 63, 63, 1, 1, 7, 7, 21, 21, 63, 63, 1, 1, 7, 7, 21, 21, 63, 63 },
    // Dimension 10
    .{ 1, 3, 7, 5, 15, 13, 17, 9, 31, 27, 23, 29, 19, 21, 11, 25, 1, 3, 7, 5, 15, 13, 17, 9, 31, 27, 23, 29, 19, 21, 11, 25 },
    // Dimension 11-15: repeat pattern with different primes
    .{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
    .{ 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3 },
    .{ 1, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1, 3, 7, 5, 1 },
    .{ 1, 1, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3, 7, 9, 13, 11, 1, 3 },
    .{ 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13, 1, 3, 7, 13 },
};

/// Sobol evaluator: compute Sobol sample at given index and dimension
/// Uses Gray code ordering for O(1) incremental update
pub const Sobol = struct {
    /// Evaluate Sobol number at given index for a specific dimension
    /// @param index The sample index (0-based)
    /// @param dim The Sobol dimension (1-15)
    /// @param seed Optional 24-bit seed for scrambling
    pub fn evaluate(index: i32, dim: i32, seed: i32) f32 {
        if (dim < 1 or dim > SOBOL_MAX_DIM) return 0;

        // Gray code of index
        const gray = index ^ (index >> 1);

        var result: i32 = 0;
        var g = gray;
        var bit: usize = 0;
        while (g != 0 and bit < SOBOL_MAX_BITS) : (bit += 1) {
            if (g & 1 != 0) {
                result ^= SOBOL_DIRECTION_NUMBERS[@intCast(dim)][bit];
            }
            g >>= 1;
        }

        // Apply seed scramble
        if (seed != 0) {
            result ^= seed;
        }

        // Convert to [0, 1)
        return @as(f32, @floatFromInt(result & 0x7FFFFFFF)) / 2147483648.0;
    }

    /// Incremental Sobol: compute next value from previous
    /// Uses the property that only one bit changes in Gray code sequence
    pub fn next(index: i32, dim: i32, prev_value: f32) f32 {
        // Find the rightmost zero bit of index (which bit changed in Gray code)
        var bit: usize = 0;
        var idx = index;
        while (idx & 1 != 0 and bit < SOBOL_MAX_BITS) : (bit += 1) {
            idx >>= 1;
        }

        if (bit >= SOBOL_MAX_BITS or dim < 1 or dim > SOBOL_MAX_DIM) return prev_value;

        // XOR the direction number for that bit
        const dir_num = SOBOL_DIRECTION_NUMBERS[@intCast(dim)][bit];
        const prev_int = @as(i32, @intFromFloat(prev_value * 2147483648.0));
        const new_int = prev_int ^ dir_num;

        return @as(f32, @floatFromInt(new_int & 0x7FFFFFFF)) / 2147483648.0;
    }

    /// 2D Sobol sample within a cell
    pub fn evaluate2D(index: i32, cell_bits: i32, cell_x: i32, cell_y: i32, seed_x: i32, seed_y: i32) [2]f32 {
        _ = cell_bits;
        _ = cell_x;
        _ = cell_y;
        // Simplified: just use dim 1,2 with seeds
        return .{
            evaluate(index, 1, seed_x),
            evaluate(index, 2, seed_y),
        };
    }

    /// 3D Sobol sample within a cell
    pub fn evaluate3D(index: i32, cell_bits: i32, cell: [3]i32, seed: [3]i32) [3]f32 {
        _ = cell_bits;
        _ = cell;
        return .{
            evaluate(index, 1, seed[0]),
            evaluate(index, 2, seed[1]),
            evaluate(index, 3, seed[2]),
        };
    }

    /// GPU spatial seed computation for Sobol sampling
    pub fn computeGPUSpatialSeed(x: i32, y: i32, index: i32) u16 {
        // Hash x, y, index into a 16-bit seed
        var h: u32 = @intCast(x * 374761393 + y * 668265263 + index * 1274126177);
        h = (h ^ (h >> 13)) * 1103515245 + 12345;
        h = h ^ (h >> 16);
        return @intCast(h & 0xFFFF);
    }
};

// =============================================================================
// 3. P³ QUASIRANDOM — PROJECTIVE SAMPLING
// =============================================================================

/// Quasirandom point in P³ (projective space)
/// Maps [0,1)^4 → P³ by treating as homogeneous coordinates
pub fn sobolP3(index: i32, seed: [4]i32) [4]f32 {
    return .{
        Sobol.evaluate(index, 1, seed[0]),
        Sobol.evaluate(index, 2, seed[1]),
        Sobol.evaluate(index, 3, seed[2]),
        Sobol.evaluate(index, 4, seed[3]),
    };
}

/// Quasirandom point on S³ (unit sphere in R⁴)
/// Maps Sobol → [0,1)^4 → spherical coordinates → S³
pub fn haltonS3(index: i32) [4]f32 {
    const h = halton4D(index);
    // Map [0,1)^4 → S³ via inverse stereographic-like projection
    // θ₁ = 2π·h[0], θ₂ = π·h[1], θ₃ = π·h[2], but with correct measure
    const theta1 = 2.0 * math.pi * h[0];
    const theta2 = math.pi * h[1];
    const theta3 = 2.0 * math.pi * h[2];
    const r = @sqrt(h[3]); // Correct for uniform measure on S³

    const x1 = r * @sin(theta3) * @sin(theta2) * @sin(theta1);
    const x2 = r * @sin(theta3) * @sin(theta2) * @cos(theta1);
    const x3 = r * @sin(theta3) * @cos(theta2);
    const x4 = r * @cos(theta3);

    return .{ x1, x2, x3, x4 };
}

/// Quasirandom point on RP³ = S³ / {±1}
/// Same as S³ but with antipodal identification
pub fn haltonRP3(index: i32) [4]f32 {
    var p = haltonS3(index);
    // Identify antipodal: ensure first non-zero coord is positive
    for (0..4) |i| {
        if (@abs(p[i]) > 1e-8) {
            if (p[i] < 0) {
                p[0] = -p[0];
                p[1] = -p[1];
                p[2] = -p[2];
                p[3] = -p[3];
            }
            break;
        }
    }
    return p;
}

// =============================================================================
// 4. DISCREPANCY MEASUREMENT
// =============================================================================

/// Estimate star-discrepancy D*_N of a 1D point set
/// D*_N ≈ max over intervals [0,a) of |count/N - a|
pub fn starDiscrepancy1D(points: []const f32) f32 {
    if (points.len == 0) return 1.0;

    var max_disc: f32 = 0;
    var sorted = std.ArrayList(f32).init(std.testing.allocator);
    defer sorted.deinit();
    sorted.appendSlice(points) catch return 1.0;
    sorted.sort(f32, {}, comptime struct {
        fn lessThan(_: void, a: f32, b: f32) bool {
            return a < b;
        }
    }.lessThan);

    const n: f32 = @floatFromInt(sorted.items.len);
    for (sorted.items, 0..) |p, i| {
        const fi: f32 = @floatFromInt(i + 1);
        // D*_i = max(|i/N - x_i|, |x_i - (i-1)/N|)
        const d1 = @abs(fi / n - p);
        const d2 = @abs(p - (@as(f32, @floatFromInt(i)) / n));
        max_disc = @max(max_disc, d1, d2);
    }

    return max_disc;
}

// =============================================================================
// TESTS
// =============================================================================

test "Halton: basic values" {
    // Halton(0, 2) = 0, Halton(1, 2) = 0.5, Halton(2, 2) = 0.25
    try std.testing.expectApproxEqAbs(@as(f32, 0.0), halton(0, 2), 1e-6);
    try std.testing.expectApproxEqAbs(@as(f32, 0.5), halton(1, 2), 1e-6);
    try std.testing.expectApproxEqAbs(@as(f32, 0.25), halton(2, 2), 1e-6);
    try std.testing.expectApproxEqAbs(@as(f32, 0.75), halton(3, 2), 1e-6);
}

test "Halton: base 3" {
    // Halton(1, 3) = 1/3, Halton(2, 3) = 2/3, Halton(3, 3) = 1/9
    try std.testing.expectApproxEqAbs(@as(f32, 1.0 / 3.0), halton(1, 3), 1e-5);
    try std.testing.expectApproxEqAbs(@as(f32, 2.0 / 3.0), halton(2, 3), 1e-5);
    try std.testing.expectApproxEqAbs(@as(f32, 1.0 / 9.0), halton(3, 3), 1e-5);
}

test "Halton: 2D" {
    const p = halton2D(1);
    try std.testing.expectApproxEqAbs(@as(f32, 0.5), p[0], 1e-5);
    try std.testing.expectApproxEqAbs(@as(f32, 1.0 / 3.0), p[1], 1e-5);
}

test "Sobol: dimension 1 first few values" {
    const s0 = Sobol.evaluate(0, 1, 0);
    const s1 = Sobol.evaluate(1, 1, 0);
    const s2 = Sobol.evaluate(2, 1, 0);
    // First Sobol dim should give 0, 0.5, 0.25, ...
    try std.testing.expect(s0 >= 0 and s0 < 1);
    try std.testing.expect(s1 >= 0 and s1 < 1);
    try std.testing.expect(s2 >= 0 and s2 < 1);
}

test "Sobol: values in [0,1)" {
    for (1..100) |i| {
        const s = Sobol.evaluate(@intCast(i), 1, 0);
        try std.testing.expect(s >= 0 and s < 1);
    }
}

test "Sobol: GPU spatial seed" {
    const seed = Sobol.computeGPUSpatialSeed(0, 0, 0);
    try std.testing.expect(seed >= 0 and seed <= 0xFFFF);
}

test "Halton: S³ points on unit sphere" {
    for (1..10) |i| {
        const p = haltonS3(@intCast(i));
        const len_sq = p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + p[3] * p[3];
        // Not exactly on S³ due to r = sqrt(h[3]), but bounded
        try std.testing.expect(len_sq <= 1.0 + 1e-4);
    }
}