zyga 0.5.1

ZYGA zero-knowledge proof system - CLI and library for generating ZK proofs
Documentation

/// Evaluate Lagrange polynomial for constraint matrix at point x
/// This interpolates the constraint matrix M to create a polynomial that passes through
/// all constraint points (1, M[0]), (2, M[1]), ..., (n, M[n-1])
#[allow(clippy::needless_range_loop)] // Indices are used for mathematical calculations
pub fn lagrange_interpolate(matrix: &[Vec<f64>], x: f64) -> Vec<f64> {
    let n_constraints = matrix.len();
    let n_vars = if n_constraints > 0 {
        matrix[0].len()
    } else {
        0
    };

    let mut result = vec![0.0; n_vars];

    // For each variable j
    for j in 0..n_vars {
        let mut p = 0.0;

        // Lagrange interpolation: sum over all constraints
        for i in 0..n_constraints {
            let mut lagrange_basis = 1.0;

            // Compute Lagrange basis polynomial L_i(x)
            // L_i(x) = product of (x - c)/(i+1 - c) for all c != i+1
            for c in 1..=n_constraints {
                if c != i + 1 {
                    lagrange_basis *= (x - c as f64) / ((i + 1) as f64 - c as f64);
                }
            }

            // Add contribution of this constraint
            p += lagrange_basis * matrix[i][j];
        }

        result[j] = p;
    }

    result
}

/// Create a polynomial function that interpolates the constraint matrix
pub struct LagrangePolynomial {
    matrix: Vec<Vec<f64>>,
}

impl LagrangePolynomial {
    pub fn new(matrix: Vec<Vec<f64>>) -> Self {
        LagrangePolynomial { matrix }
    }

    /// Evaluate the polynomial at point x
    pub fn evaluate(&self, x: f64) -> Vec<f64> {
        lagrange_interpolate(&self.matrix, x)
    }
}

/// Compute the vanishing polynomial Z_n(x) = (x-1)(x-2)...(x-n)
pub fn vanishing_polynomial(n_constraints: usize, x: f64) -> f64 {
    let mut result = 1.0;
    for i in 1..=n_constraints {
        result *= x - i as f64;
    }
    result
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_lagrange_interpolation() {
        // Test with a simple 2x3 matrix
        let matrix = vec![vec![1.0, 0.0, 0.0], vec![0.0, 1.0, 0.0]];

        // At x=1, should get first row
        let result = lagrange_interpolate(&matrix, 1.0);
        assert!((result[0] - 1.0).abs() < 1e-10);
        assert!((result[1] - 0.0).abs() < 1e-10);

        // At x=2, should get second row
        let result = lagrange_interpolate(&matrix, 2.0);
        assert!((result[0] - 0.0).abs() < 1e-10);
        assert!((result[1] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_lagrange_matches_python() {
        // Test with same 3x3 identity matrix as Python
        let matrix = vec![
            vec![1.0, 0.0, 0.0],
            vec![0.0, 1.0, 0.0],
            vec![0.0, 0.0, 1.0],
        ];

        // Test at x=4 (should match Python: [1.0, -3.0, 3.0])
        let result = lagrange_interpolate(&matrix, 4.0);
        println!("Rust LP(M, 4) = {:?}", result);
        assert!(
            (result[0] - 1.0).abs() < 1e-10,
            "Expected result[0] = 1.0, got {}",
            result[0]
        );
        assert!(
            (result[1] - (-3.0)).abs() < 1e-10,
            "Expected result[1] = -3.0, got {}",
            result[1]
        );
        assert!(
            (result[2] - 3.0).abs() < 1e-10,
            "Expected result[2] = 3.0, got {}",
            result[2]
        );
    }

    #[test]
    fn test_vanishing_polynomial() {
        // Z_3(x) = (x-1)(x-2)(x-3)
        let z = vanishing_polynomial(3, 4.0);
        assert!((z - 6.0).abs() < 1e-10); // (4-1)(4-2)(4-3) = 3*2*1 = 6

        // Should be zero at constraint points
        assert!(vanishing_polynomial(3, 1.0).abs() < 1e-10);
        assert!(vanishing_polynomial(3, 2.0).abs() < 1e-10);
        assert!(vanishing_polynomial(3, 3.0).abs() < 1e-10);
    }
}