#!/usr/bin/env julia

"""
Exact verifier for the frequency-22 radial Class-I cubature theorem.

Every mathematical equality and sign decision is made in Q(sqrt(5)).
BigFloat is used only to print short orientation decimals.

Verified statements:

1. The fixed Class-I grid X_22 supports strictly positive H3-orbit-constant
   cubature through algebraic degree 69.
2. The unique degree-70/71 invariant solution has a negative mass on orbit
   (21,1,0).
3. An exact Farkas/cardinal dual has the same negative spherical mean.

By symmetrization, statement 2 excludes every nonnegative degree-70 rule on
the fixed node set, even when the original node weights are not orbit-constant.
"""

using SHA
using Printf

const QQ = Rational{BigInt}

rat(n::Integer) = BigInt(n) // BigInt(1)

struct K5
    a::QQ
    b::QQ
end

K5() = K5(rat(0), rat(0))
K5(a::Integer) = K5(rat(a), rat(0))
K5(a::QQ) = K5(a, rat(0))
K5(a::Integer, b::Integer) = K5(rat(a), rat(b))

Base.zero(::Type{K5}) = K5()
Base.one(::Type{K5}) = K5(1)
Base.iszero(x::K5) = iszero(x.a) && iszero(x.b)
Base.:(==)(x::K5, y::K5) = x.a == y.a && x.b == y.b
Base.:+(x::K5, y::K5) = K5(x.a + y.a, x.b + y.b)
Base.:-(x::K5) = K5(-x.a, -x.b)
Base.:-(x::K5, y::K5) = x + (-y)
Base.:*(x::K5, y::K5) =
    K5(x.a * y.a + 5 * x.b * y.b, x.a * y.b + x.b * y.a)

function Base.inv(x::K5)
    d = x.a * x.a - 5 * x.b * x.b
    iszero(d) && throw(DivideError())
    K5(x.a / d, -x.b / d)
end

Base.:/(x::K5, y::K5) = x * inv(y)

function Base.:^(x::K5, n::Integer)
    n < 0 && return inv(x) ^ (-n)
    result = K5(1)
    base = x
    exponent = n
    while exponent > 0
        isodd(exponent) && (result *= base)
        base *= base
        exponent >>= 1
    end
    result
end

"""Exact sign under the real embedding sqrt(5)>0."""
function exact_sign(x::K5)
    iszero(x.b) && return sign(x.a)
    x.a >= 0 && x.b >= 0 && return 1
    x.a <= 0 && x.b <= 0 && return -1
    aa = x.a * x.a
    fivebb = 5 * x.b * x.b
    aa == fivebb && return 0
    x.a > 0 ? (aa > fivebb ? 1 : -1) : (aa > fivebb ? -1 : 1)
end

function orientation_decimal(x::K5)
    setprecision(BigFloat, 128) do
        a = BigFloat(numerator(x.a)) / BigFloat(denominator(x.a))
        b = BigFloat(numerator(x.b)) / BigFloat(denominator(x.b))
        a + b * sqrt(BigFloat(5))
    end
end

const SQRT5 = K5(0, 1)
const PHI = (K5(1) + SQRT5) / K5(2)
const R2 = PHI + K5(2)

const Point3 = NTuple{3,K5}
const Exponent3 = NTuple{3,Int}
const Polynomial = Dict{Exponent3,K5}

const VERTICES = Point3[]
for s1 in (1, -1), s2 in (1, -1)
    push!(VERTICES, (K5(0), K5(s1), K5(s2) * PHI))
    push!(VERTICES, (K5(s1), K5(s2) * PHI, K5(0)))
    push!(VERTICES, (K5(s1) * PHI, K5(0), K5(s2)))
end

const FACE = (
    (K5(0), K5(1), PHI),
    (K5(-1), PHI, K5(0)),
    (K5(1), PHI, K5(0)),
)

function dot3(x::Point3, y::Point3)
    total = K5()
    for i in 1:3
        total += x[i] * y[i]
    end
    total
end

function poly_add(left::Polynomial, right::Polynomial)
    result = copy(left)
    for (exponent, coefficient) in right
        value = get(result, exponent, K5()) + coefficient
        iszero(value) ? delete!(result, exponent) : (result[exponent] = value)
    end
    result
end

function poly_scale(poly::Polynomial, scalar::K5)
    result = Polynomial()
    for (exponent, coefficient) in poly
        value = scalar * coefficient
        !iszero(value) && (result[exponent] = value)
    end
    result
end

function poly_multiply(left::Polynomial, right::Polynomial)
    result = Polynomial()
    for ((i, j, k), a) in left
        for ((u, v, w), b) in right
            exponent = (i + u, j + v, k + w)
            result[exponent] = get(result, exponent, K5()) + a * b
        end
    end
    filter!(pair -> !iszero(last(pair)), result)
    result
end

function multinomial3(n::Int, i::Int, j::Int, k::Int)
    factorial(big(n)) ÷
        (factorial(big(i)) * factorial(big(j)) * factorial(big(k)))
end

function linear_power(vector::Point3, n::Int)
    result = Polynomial()
    vx, vy, vz = vector
    for i in 0:n
        for j in 0:(n - i)
            k = n - i - j
            coefficient =
                K5(multinomial3(n, i, j, k)) * vx^i * vy^j * vz^k
            !iszero(coefficient) && (result[(i, j, k)] = coefficient)
        end
    end
    result
end

function invariant_polynomial(n::Int)
    result = Polynomial()
    for vertex in VERTICES
        result = poly_add(result, linear_power(vertex, n))
    end
    poly_scale(result, R2 ^ (-(n ÷ 2)))
end

const U6_POLY = invariant_polynomial(6)
const U10_POLY = invariant_polynomial(10)

function odd_double_factorial(n::Int)
    result = BigInt(1)
    for value in n:-2:1
        result *= value
    end
    result
end

function sphere_monomial_mean(i::Int, j::Int, k::Int)
    (isodd(i) || isodd(j) || isodd(k)) && return rat(0)
    a, b, c = i ÷ 2, j ÷ 2, k ÷ 2
    numerator_value =
        odd_double_factorial(2 * a - 1) *
        odd_double_factorial(2 * b - 1) *
        odd_double_factorial(2 * c - 1)
    denominator_value = odd_double_factorial(2(a + b + c) + 1)
    numerator_value // denominator_value
end

function sphere_integral(poly::Polynomial)
    result = K5()
    for ((i, j, k), coefficient) in poly
        result += coefficient * K5(sphere_monomial_mean(i, j, k))
    end
    result
end

const POWERS6 = Dict{Int,Polynomial}(
    0 => Polynomial((0, 0, 0) => K5(1)),
    1 => U6_POLY,
)
const POWERS10 = Dict{Int,Polynomial}(
    0 => Polynomial((0, 0, 0) => K5(1)),
    1 => U10_POLY,
)

function invariant_moment(r::Int, s::Int)
    while maximum(keys(POWERS6)) < r
        n = maximum(keys(POWERS6)) + 1
        POWERS6[n] = poly_multiply(POWERS6[n - 1], U6_POLY)
    end
    while maximum(keys(POWERS10)) < s
        n = maximum(keys(POWERS10)) + 1
        POWERS10[n] = poly_multiply(POWERS10[n - 1], U10_POLY)
    end
    sphere_integral(poly_multiply(POWERS6[r], POWERS10[s]))
end

function basis_pairs(degree::Int)
    pairs = Tuple{Int,Int}[]
    for s in 0:(degree ÷ 10)
        for r in 0:((degree - 10 * s) ÷ 6)
            push!(pairs, (r, s))
        end
    end
    sort!(pairs, by = pair -> (6 * pair[1] + 10 * pair[2], pair[2], pair[1]))
end

function partitions3(total::Int)
    result = NTuple{3,Int}[]
    for a in total:-1:0
        for b in min(a, total - a):-1:0
            c = total - a - b
            0 <= c <= b && push!(result, (a, b, c))
        end
    end
    result
end

function orbit_invariants(representative::NTuple{3,Int})
    a, b, c = representative
    point = ntuple(
        i -> K5(a) * FACE[1][i] + K5(b) * FACE[2][i] + K5(c) * FACE[3][i],
        3,
    )
    norm_squared = dot3(point, point)
    sum6 = K5()
    sum10 = K5()
    for vertex in VERTICES
        product = dot3(point, vertex)
        sum6 += product^6
        sum10 += product^10
    end
    u6 = sum6 / (R2^3 * norm_squared^3)
    u10 = sum10 / (R2^5 * norm_squared^5)
    (u6, u10)
end

function exact_rref(matrix::Matrix{K5}, right_hand_side::Vector{K5})
    row_count, column_count = size(matrix)
    augmented = Matrix{K5}(undef, row_count, column_count + 1)
    augmented[:, 1:column_count] .= matrix
    augmented[:, column_count + 1] .= right_hand_side

    pivot_row = 1
    pivots = Int[]
    for column in 1:column_count
        selected_row = 0
        for row in pivot_row:row_count
            if !iszero(augmented[row, column])
                selected_row = row
                break
            end
        end
        selected_row == 0 && continue
        if selected_row != pivot_row
            for j in 1:(column_count + 1)
                augmented[pivot_row, j], augmented[selected_row, j] =
                    augmented[selected_row, j], augmented[pivot_row, j]
            end
        end

        inverse_pivot = inv(augmented[pivot_row, column])
        for j in 1:(column_count + 1)
            augmented[pivot_row, j] *= inverse_pivot
        end

        for row in 1:row_count
            row == pivot_row && continue
            factor = augmented[row, column]
            iszero(factor) && continue
            for j in 1:(column_count + 1)
                augmented[row, j] -= factor * augmented[pivot_row, j]
            end
        end

        push!(pivots, column)
        pivot_row += 1
        pivot_row > row_count && break
    end

    for row in pivot_row:row_count
        all_zero = all(iszero(augmented[row, column]) for column in 1:column_count)
        all_zero && !iszero(augmented[row, column_count + 1]) &&
            error("inconsistent exact moment system")
    end
    augmented, pivots
end

function choose_solution(
    column_count::Int,
    reduced::Matrix{K5},
    pivots::Vector{Int},
    free_values::Vector{K5},
)
    pivot_set = Set(pivots)
    free_columns = [column for column in 1:column_count if !(column in pivot_set)]
    length(free_columns) == length(free_values) ||
        error("wrong number of free values")

    solution = fill(K5(), column_count)
    for (column, value) in zip(free_columns, free_values)
        solution[column] = value
    end
    for (row, column) in enumerate(pivots)
        value = reduced[row, end]
        for free_column in free_columns
            value -= reduced[row, free_column] * solution[free_column]
        end
        solution[column] = value
    end
    solution, free_columns
end

function verify_system(
    matrix::Matrix{K5},
    right_hand_side::Vector{K5},
    solution::Vector{K5},
)
    for row in axes(matrix, 1)
        lhs = K5()
        for column in axes(matrix, 2)
            lhs += matrix[row, column] * solution[column]
        end
        lhs == right_hand_side[row] || return false
    end
    true
end

function moment_system(frequency::Int, degree::Int)
    representatives = partitions3(frequency)
    pairs = basis_pairs(degree)
    values = [orbit_invariants(rep) for rep in representatives]
    matrix = Matrix{K5}(undef, length(pairs), length(representatives))
    for (row, (r, s)) in enumerate(pairs)
        for column in eachindex(representatives)
            u6, u10 = values[column]
            matrix[row, column] = u6^r * u10^s
        end
    end
    rhs = [invariant_moment(r, s) for (r, s) in pairs]
    reduced, pivots = exact_rref(matrix, rhs)
    representatives, pairs, values, matrix, rhs, reduced, pivots
end

canonical_fraction(value::QQ) =
    string(numerator(value), "/", denominator(value))
canonical_k5(value::K5) =
    string(canonical_fraction(value.a), ";", canonical_fraction(value.b))

function solution_digest(representatives, weights)
    records = String[]
    for ((a, b, c), weight) in zip(representatives, weights)
        push!(records, string(a, ",", b, ",", c, ":", canonical_k5(weight)))
    end
    bytes2hex(sha256(codeunits(join(records, "\n"))))
end

const POSITIVE_FREE_VALUES = K5[
    K5(BigInt(11607197) // BigInt(394254842)),
    K5(BigInt(2312954) // BigInt(161188803)),
    K5(BigInt(7156114) // BigInt(472522843)),
]

const EXPECTED_DEGREE69_SHA256 =
    "0865572f76c7b8c9816ce937a79361c8926abe785ac3f41b06e4f289665bf340"
const EXPECTED_DEGREE71_SHA256 =
    "c1e620fc9cf6d99623423cf82450540a5991679e8ded6632bd959914068949f2"

function main()
    started = time()
    frequency = 22

    reps, pairs69, values, matrix69, rhs69, reduced69, pivots69 =
        moment_system(frequency, 69)
    weights69, free_columns =
        choose_solution(length(reps), reduced69, pivots69, POSITIVE_FREE_VALUES)

    @assert length(reps) == 52
    @assert length(pairs69) == length(pivots69) == 49
    @assert free_columns == [50, 51, 52]
    @assert reps[free_columns] == [(9, 7, 6), (8, 8, 6), (8, 7, 7)]
    @assert verify_system(matrix69, rhs69, weights69)
    @assert foldl(+, weights69; init = K5()) == K5(1)
    @assert all(exact_sign(weight) > 0 for weight in weights69)

    minimum_index = 1
    for index in 2:length(weights69)
        exact_sign(weights69[index] - weights69[minimum_index]) < 0 &&
            (minimum_index = index)
    end
    @assert reps[minimum_index] == (11, 11, 0)
    digest69 = solution_digest(reps, weights69)
    @assert digest69 == EXPECTED_DEGREE69_SHA256

    pairs71 = basis_pairs(71)
    old_pairs = Set(pairs69)
    new_pairs = [pair for pair in pairs71 if !(pair in old_pairs)]
    @assert new_pairs == [(10, 1), (5, 4), (0, 7)]

    extra_matrix = Matrix{K5}(undef, length(new_pairs), length(reps))
    for (row, (r, s)) in enumerate(new_pairs)
        for column in eachindex(reps)
            u6, u10 = values[column]
            extra_matrix[row, column] = u6^r * u10^s
        end
    end
    extra_rhs = [invariant_moment(r, s) for (r, s) in new_pairs]

    constants = fill(K5(), length(reps))
    coefficients = fill(K5(), length(reps), 3)
    for (k, column) in enumerate(free_columns)
        coefficients[column, k] = K5(1)
    end
    for (row, column) in enumerate(pivots69)
        constants[column] = reduced69[row, end]
        for (k, free_column) in enumerate(free_columns)
            coefficients[column, k] = -reduced69[row, free_column]
        end
    end

    three_by_three = Matrix{K5}(undef, 3, 3)
    three_rhs = fill(K5(), 3)
    for row in 1:3
        for k in 1:3
            value = K5()
            for column in eachindex(reps)
                value += extra_matrix[row, column] * coefficients[column, k]
            end
            three_by_three[row, k] = value
        end
        constant_part = K5()
        for column in eachindex(reps)
            constant_part += extra_matrix[row, column] * constants[column]
        end
        three_rhs[row] = extra_rhs[row] - constant_part
    end

    reduced3, pivots3 = exact_rref(three_by_three, three_rhs)
    @assert pivots3 == [1, 2, 3]
    forced_free, no_free = choose_solution(3, reduced3, pivots3, K5[])
    @assert isempty(no_free)

    weights71 = fill(K5(), length(reps))
    for column in eachindex(reps)
        value = constants[column]
        for k in 1:3
            value += coefficients[column, k] * forced_free[k]
        end
        weights71[column] = value
    end

    matrix71 = vcat(matrix69, extra_matrix)
    rhs71 = vcat(rhs69, extra_rhs)
    @assert size(matrix71) == (52, 52)
    @assert verify_system(matrix71, rhs71, weights71)
    @assert foldl(+, weights71; init = K5()) == K5(1)

    negative_indices =
        [index for index in eachindex(weights71) if exact_sign(weights71[index]) < 0]
    @assert length(negative_indices) == 1
    bad_index = only(negative_indices)
    bad_rep = reps[bad_index]
    bad_weight = weights71[bad_index]
    @assert bad_rep == (21, 1, 0)

    # permutedims reorders the matrix without attempting scalar transpose on
    # K5, which is a field element but intentionally is not a subtype of Number.
    transpose_matrix = permutedims(matrix71)
    selector = [index == bad_index ? K5(1) : K5() for index in eachindex(reps)]
    reduced_dual, pivots_dual = exact_rref(transpose_matrix, selector)
    @assert length(pivots_dual) == 52
    dual_coefficients, dual_free =
        choose_solution(52, reduced_dual, pivots_dual, K5[])
    @assert isempty(dual_free)
    @assert verify_system(transpose_matrix, selector, dual_coefficients)

    dual_mean = K5()
    for index in eachindex(rhs71)
        dual_mean += dual_coefficients[index] * rhs71[index]
    end
    @assert dual_mean == bad_weight
    @assert exact_sign(dual_mean) < 0

    digest71 = solution_digest(reps, weights71)
    @assert digest71 == EXPECTED_DEGREE71_SHA256

    println("OMEGA NAUGHT - EXACT JULIA CLASS-I f=22 CERTIFICATE")
    println("Arithmetic: Q(sqrt(5)); equality/sign checks: exact")
    println("nodes=$(10 * frequency^2 + 2)  H3_orbits=$(length(reps))")
    println(
        "degree_69_conditions=$(length(pairs69))  rank=$(length(pivots69))  " *
        "affine_dimension=$(length(free_columns))",
    )
    println("free_orbits=$(reps[free_columns])")
    println("free_masses=$(POSITIVE_FREE_VALUES)")
    @printf(
        "degree_69: EXACT and STRICTLY POSITIVE; minimum_orbit_mass=%.16g at %s\n",
        Float64(orientation_decimal(weights69[minimum_index])),
        string(reps[minimum_index]),
    )
    println("degree_69_solution_sha256=$digest69")
    println("new_degree_70_invariants=$new_pairs")
    @printf(
        "degree_70/71: UNIQUE exact solution; negative_orbit_mass=%.16g at %s\n",
        Float64(orientation_decimal(bad_weight)),
        string(bad_rep),
    )
    println("degree_71_solution_sha256=$digest71")
    println("Farkas/cardinal dual: EXACT PASS; mean equals negative forced mass")
    println("THEOREM: maximum nonnegative degree on fixed X_22 is 69")
    @printf("CERTIFICATE PASS (%.2f seconds)\n", time() - started)
    true
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end