//
//  CholeskySolver.swift
//  Metrika
//
//  Author:  Simon-Pierre Boucher
//  Contact: contact@spboucher.ai
//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.
//

/// Plain-Swift Cholesky solve for the k×k normal equations of the GPU
/// bootstrap. k is the regressor count (single digits in practice), so a
/// dependency-free O(k³) routine beats reaching for LAPACK — and ZQStats
/// is the only module allowed to import Accelerate (CLAUDE.md §3).
enum CholeskySolver {
    struct SingularMatrix: Error {}

    /// Solves A·x = b for symmetric positive-definite A (column-major,
    /// k×k). Throws on non-positive-definite input.
    static func solve(_ a: [Double], k: Int, rhs: [Double]) throws -> [Double] {
        precondition(a.count == k * k && rhs.count == k)

        // Lower-triangular factor L with A = L·Lᵀ.
        var l = [Double](repeating: 0, count: k * k)
        for j in 0..<k {
            var diagonal = a[j * k + j]
            for p in 0..<j {
                diagonal -= l[p * k + j] * l[p * k + j]
            }
            guard diagonal > 0 else { throw SingularMatrix() }
            let root = diagonal.squareRoot()
            l[j * k + j] = root
            for i in (j + 1)..<k {
                var value = a[j * k + i]
                for p in 0..<j {
                    value -= l[p * k + i] * l[p * k + j]
                }
                l[j * k + i] = value / root
            }
        }

        // Forward substitution L·z = b.
        var z = [Double](repeating: 0, count: k)
        for i in 0..<k {
            var value = rhs[i]
            for p in 0..<i {
                value -= l[p * k + i] * z[p]
            }
            z[i] = value / l[i * k + i]
        }

        // Back substitution Lᵀ·x = z.
        var x = [Double](repeating: 0, count: k)
        for i in stride(from: k - 1, through: 0, by: -1) {
            var value = z[i]
            for p in (i + 1)..<k {
                value -= l[i * k + p] * x[p]
            }
            x[i] = value / l[i * k + i]
        }
        return x
    }
}