trop/tri/determinant.tri

// ---
// tags: trop, trident
// crystal-type: circuit
// crystal-domain: comp
// ---

//  Tropical determinant (minimum-weight perfect matching).
//
//  trop_det(A) = min over all permutations sigma of sum_i A[i][sigma(i)].
//
//  Naive O(n!) enumeration via Heap's algorithm.
//  Bounded to n <= 8 in Trident (8! = 40320 iterations, feasible in circuit).

module trop.determinant

use trop.element.{INF, Tropical, add, mul, is_inf, from_u32}
use trop.matrix.{MAX_N, TropMatrix, get}

/// Maximum dimension for naive determinant (8! = 40320).
const MAX_DET_DIM: U32 = 8

/// Evaluate the cost of a single permutation and return tropical min with best.
fn eval_perm(a: TropMatrix, perm: [U32; 8], n: U32, best: U32) -> U32 {
    let mut cost: Tropical = Tropical { val: 0 }
    for i in 0..8 {
        if i < n {
            let entry: Tropical = get(a, i, perm[i])
            cost = mul(cost, entry)
            if is_inf(cost) {
                return best
            }
        }
    }
    let cur_best: Tropical = from_u32(best)
    add(cur_best, cost).val
}

/// Compute the tropical determinant of a square matrix.
///
/// trop_det(A) = min_{sigma in S_n} sum_i A[i][sigma(i)]
///
/// n must be <= 8.
pub fn determinant(a: TropMatrix) -> U32 {
    let n: U32 = a.n
    if n == 0 {
        return 0
    }

    // Initialize permutation [0, 1, 2, ..., 7]
    let mut perm: [U32; 8] = [0; 8]
    perm[0] = 0
    perm[1] = 1
    perm[2] = 2
    perm[3] = 3
    perm[4] = 4
    perm[5] = 5
    perm[6] = 6
    perm[7] = 7

    // Heap's algorithm control array
    let mut c: [U32; 8] = [0; 8]

    let mut best: U32 = INF

    // Evaluate identity permutation
    best = eval_perm(a, perm, n, best)

    // Generate remaining permutations via iterative Heap's algorithm.
    // Total permutations for n elements: n!
    // Max iterations: 8! = 40320
    let mut i: U32 = 0
    for iter in 0..40320 {
        if i < n {
            if c[i] < i {
                // Perform swap
                let even: U32 = i & 1
                if even == 0 {
                    // Swap perm[0] and perm[i]
                    let tmp: U32 = perm[0]
                    perm[0] = perm[i]
                    perm[i] = tmp
                } else {
                    // Swap perm[c[i]] and perm[i]
                    let ci: U32 = c[i]
                    let tmp: U32 = perm[ci]
                    perm[ci] = perm[i]
                    perm[i] = tmp
                }

                best = eval_perm(a, perm, n, best)

                c[i] = c[i] + 1
                i = 0
            } else {
                c[i] = 0
                i = i + 1
            }
        }
    }

    best
}