implement TIP 0001: Contiguity Argument for Memory Consistency

Fix #77

Co-authored-by: Alexander Lemmens <alexanderlemmens@hotmail.com>

triton-vm/src/table/base_matrix.rs

@@ -1,13 +1,18 @@
 use std::fmt::{Display, Formatter};
+use std::ops::Mul;
 
 use itertools::Itertools;
-use num_traits::Zero;
+use num_traits::{One, Zero};
 use twenty_first::shared_math::b_field_element::BFieldElement;
+use twenty_first::shared_math::polynomial::Polynomial;
 use twenty_first::shared_math::traits::FiniteField;
 use twenty_first::shared_math::x_field_element::XFieldElement;
 
 use crate::instruction::AnInstruction::*;
 use crate::table::table_column::ProcessorExtTableColumn::*;
+use crate::table::table_column::RamBaseTableColumn::{
+    BezoutCoefficientPolynomialCoefficient0, BezoutCoefficientPolynomialCoefficient1, RAMP,
+};
 use crate::table::table_column::{
     InstructionBaseTableColumn, OpStackBaseTableColumn, ProcessorBaseTableColumn,
     ProgramBaseTableColumn, RamBaseTableColumn,
@@ -173,6 +178,75 @@ impl BaseMatrices {
         {
             ram_matrix[idx][usize::from(RamBaseTableColumn::InverseOfRampDifference)] = inverse;
         }
+
+        // compute Bézout coefficient polynomials
+        let all_ramps = ram_matrix
+            .iter()
+            .map(|row| row[usize::from(RAMP)])
+            .dedup()
+            .collect_vec();
+        let num_of_different_ramps = all_ramps.len();
+        let polynomial_with_ramps_as_roots = all_ramps
+            .into_iter()
+            .map(|ramp| Polynomial::new(vec![-ramp, BFieldElement::one()]))
+            .fold(
+                Polynomial::from_constant(BFieldElement::one()),
+                Polynomial::mul,
+            );
+
+        let formal_derivative_coefficients = polynomial_with_ramps_as_roots
+            .coefficients
+            .iter()
+            .enumerate()
+            .map(|(i, &coefficient)| BFieldElement::new(i as u64) * coefficient)
+            .skip(1)
+            .collect_vec();
+        let formal_derivative = Polynomial::new(formal_derivative_coefficients);
+
+        let (gcd, bezout_0, bezout_1) =
+            XFieldElement::xgcd(polynomial_with_ramps_as_roots, formal_derivative);
+        assert_eq!(Polynomial::from_constant(BFieldElement::one()), gcd);
+        assert!(
+            bezout_0.degree() < num_of_different_ramps as isize,
+            "The Bézout coefficient must be of degree at most {}.",
+            num_of_different_ramps - 1
+        );
+        assert!(
+            bezout_1.degree() <= num_of_different_ramps as isize,
+            "The Bézout coefficient must be of degree at most {num_of_different_ramps}."
+        );
+
+        let mut bezout_coefficient_polynomial_coefficient_0 = bezout_0.coefficients;
+        bezout_coefficient_polynomial_coefficient_0
+            .resize(num_of_different_ramps, BFieldElement::zero());
+
+        let mut bezout_coefficient_polynomial_coefficient_1 = bezout_1.coefficients;
+        bezout_coefficient_polynomial_coefficient_1
+            .resize(num_of_different_ramps, BFieldElement::zero());
+
+        // take care of first row
+        let mut current_bcpc_0 = bezout_coefficient_polynomial_coefficient_0.pop().unwrap();
+        let mut current_bcpc_1 = bezout_coefficient_polynomial_coefficient_1.pop().unwrap();
+        ram_matrix.first_mut().unwrap()[usize::from(BezoutCoefficientPolynomialCoefficient0)] =
+            current_bcpc_0;
+        ram_matrix.first_mut().unwrap()[usize::from(BezoutCoefficientPolynomialCoefficient1)] =
+            current_bcpc_1;
+
+        // fill in rest of table
+        let mut previous_ramp = ram_matrix.first().unwrap()[usize::from(RAMP)];
+        for row in ram_matrix.iter_mut().skip(1) {
+            if previous_ramp != row[usize::from(RAMP)] {
+                current_bcpc_0 = bezout_coefficient_polynomial_coefficient_0.pop().unwrap();
+                current_bcpc_1 = bezout_coefficient_polynomial_coefficient_1.pop().unwrap();
+                previous_ramp = row[usize::from(RAMP)]
+            }
+            row[usize::from(BezoutCoefficientPolynomialCoefficient0)] = current_bcpc_0;
+            row[usize::from(BezoutCoefficientPolynomialCoefficient1)] = current_bcpc_1;
+        }
+
+        assert_eq!(0, bezout_coefficient_polynomial_coefficient_0.len());
+        assert_eq!(0, bezout_coefficient_polynomial_coefficient_1.len());
+
         ram_matrix
     }
 

triton-vm/src/table/challenges.rs

@@ -26,7 +26,7 @@ pub struct AllChallenges {
 }
 
 impl AllChallenges {
-    pub const TOTAL_CHALLENGES: usize = 127;
+    pub const TOTAL_CHALLENGES: usize = 128;
 
     pub fn create_challenges(mut weights: Vec<XFieldElement>) -> Self {
         let processor_table_challenges = ProcessorTableChallenges {
@@ -106,6 +106,7 @@ impl AllChallenges {
         };
 
         let ram_table_challenges = RamTableChallenges {
+            bezout_relation_sample_point: weights.pop().unwrap(),
             processor_perm_row_weight: processor_table_challenges.ram_perm_row_weight,
             clk_weight: processor_table_challenges.ram_table_clk_weight,
             ramv_weight: processor_table_challenges.ram_table_ramv_weight,

triton-vm/src/table/ram_table.rs

@@ -1,5 +1,5 @@
 use itertools::Itertools;
-use num_traits::Zero;
+use num_traits::{One, Zero};
 use strum::EnumCount;
 use twenty_first::shared_math::b_field_element::BFieldElement;
 use twenty_first::shared_math::mpolynomial::{Degree, MPolynomial};
@@ -103,34 +103,61 @@ impl RamTable {
         interpolant_degree: Degree,
     ) -> ExtRamTable {
         let mut extension_matrix: Vec<Vec<XFieldElement>> = Vec::with_capacity(self.data().len());
-        let mut running_product = PermArg::default_initial();
-
+        let mut running_product_for_perm_arg = PermArg::default_initial();
+
+        // initialize columns establishing Bézout relation
+        let mut running_product_of_ramp = challenges.bezout_relation_sample_point
+            - XFieldElement::new_const(self.data().first().unwrap()[usize::from(RAMP)]);
+        let mut formal_derivative = XFieldElement::one();
+        let mut bezout_coefficient_0 = XFieldElement::zero();
+        let mut bezout_coefficient_1 = self.data().first().unwrap()
+            [usize::from(BezoutCoefficientPolynomialCoefficient1)]
+        .lift();
+
+        let mut previous_row: Option<Vec<BFieldElement>> = None;
         for row in self.data().iter() {
             let mut extension_row = [0.into(); FULL_WIDTH];
             extension_row[..BASE_WIDTH]
                 .copy_from_slice(&row.iter().map(|elem| elem.lift()).collect_vec());
 
-            let (clk, ramp, ramv) = (
-                extension_row[usize::from(CLK)],
-                extension_row[usize::from(RAMP)],
-                extension_row[usize::from(RAMV)],
-            );
-
-            let (clk_w, ramp_w, ramv_w) = (
-                challenges.clk_weight,
-                challenges.ramp_weight,
-                challenges.ramv_weight,
-            );
+            let clk = extension_row[usize::from(CLK)];
+            let ramp = extension_row[usize::from(RAMP)];
+            let ramv = extension_row[usize::from(RAMV)];
+
+            if let Some(prow) = previous_row {
+                if prow[usize::from(RAMP)] != row[usize::from(RAMP)] {
+                    let bcpc0 = extension_row[usize::from(BezoutCoefficientPolynomialCoefficient0)];
+                    let bcpc1 = extension_row[usize::from(BezoutCoefficientPolynomialCoefficient1)];
+                    let bezout_challenge = challenges.bezout_relation_sample_point;
+
+                    formal_derivative =
+                        (bezout_challenge - ramp) * formal_derivative + running_product_of_ramp;
+                    running_product_of_ramp *= bezout_challenge - ramp;
+                    bezout_coefficient_0 = bezout_coefficient_0 * bezout_challenge + bcpc0;
+                    bezout_coefficient_1 = bezout_coefficient_1 * bezout_challenge + bcpc1;
+                }
+            }
+
+            extension_row[usize::from(RunningProductOfRAMP)] = running_product_of_ramp;
+            extension_row[usize::from(FormalDerivative)] = formal_derivative;
+            extension_row[usize::from(BezoutCoefficient0)] = bezout_coefficient_0;
+            extension_row[usize::from(BezoutCoefficient1)] = bezout_coefficient_1;
+
+            // permutation argument to Processor Table
+            let clk_w = challenges.clk_weight;
+            let ramp_w = challenges.ramp_weight;
+            let ramv_w = challenges.ramv_weight;
 
             // compress multiple values within one row so they become one value
             let compressed_row_for_permutation_argument =
                 clk * clk_w + ramp * ramp_w + ramv * ramv_w;
 
-            // compute the running *product* of the compressed column (for permutation argument)
-            running_product *=
+            // compute the running product of the compressed column for permutation argument
+            running_product_for_perm_arg *=
                 challenges.processor_perm_row_weight - compressed_row_for_permutation_argument;
-            extension_row[usize::from(RunningProductPermArg)] = running_product;
+            extension_row[usize::from(RunningProductPermArg)] = running_product_for_perm_arg;
 
+            previous_row = Some(row.clone());
             extension_matrix.push(extension_row.to_vec());
         }
 
@@ -251,27 +278,43 @@ impl Extendable for RamTable {
 impl TableLike<XFieldElement> for ExtRamTable {}
 
 impl ExtRamTable {
-    fn ext_initial_constraints(
-        _challenges: &RamTableChallenges,
-    ) -> Vec<MPolynomial<XFieldElement>> {
+    fn ext_initial_constraints(challenges: &RamTableChallenges) -> Vec<MPolynomial<XFieldElement>> {
         use RamBaseTableColumn::*;
 
-        let variables: Vec<MPolynomial<XFieldElement>> =
-            MPolynomial::variables(FULL_WIDTH, 1.into());
+        let one = MPolynomial::from_constant(1.into(), FULL_WIDTH);
+        let bezout_challenge =
+            MPolynomial::from_constant(challenges.bezout_relation_sample_point, FULL_WIDTH);
+        let variables = MPolynomial::variables(FULL_WIDTH, 1.into());
+
         let clk = variables[usize::from(CLK)].clone();
         let ramp = variables[usize::from(RAMP)].clone();
         let ramv = variables[usize::from(RAMV)].clone();
+        let bcpc0 = variables[usize::from(BezoutCoefficientPolynomialCoefficient0)].clone();
+        let bcpc1 = variables[usize::from(BezoutCoefficientPolynomialCoefficient1)].clone();
+        let rp = variables[usize::from(RunningProductOfRAMP)].clone();
+        let fd = variables[usize::from(FormalDerivative)].clone();
+        let bc0 = variables[usize::from(BezoutCoefficient0)].clone();
+        let bc1 = variables[usize::from(BezoutCoefficient1)].clone();
 
-        // Cycle count clk is 0.
         let clk_is_0 = clk;
-
-        // RAM pointer ramp is 0.
-        let ramp_is_0 = ramp;
-
-        // RAM value ramv is 0.
+        let ramp_is_0 = ramp.clone();
         let ramv_is_0 = ramv;
+        let bezout_coefficient_polynomial_coefficient_0_is_0 = bcpc0;
+        let bezout_coefficient_0_is_0 = bc0;
+        let bezout_coefficient_1_is_bezout_coefficient_polynomial_coefficient_1 = bc1 - bcpc1;
+        let formal_derivative_is_1 = fd - one;
+        let running_product_is_initialized_correctly = rp - (bezout_challenge - ramp);
 
-        vec![clk_is_0, ramp_is_0, ramv_is_0]
+        vec![
+            clk_is_0,
+            ramp_is_0,
+            ramv_is_0,
+            bezout_coefficient_polynomial_coefficient_0_is_0,
+            bezout_coefficient_0_is_0,
+            bezout_coefficient_1_is_bezout_coefficient_polynomial_coefficient_1,
+            formal_derivative_is_1,
+            running_product_is_initialized_correctly,
+        ]
     }
 
     fn ext_consistency_constraints(
@@ -282,66 +325,117 @@ impl ExtRamTable {
     }
 
     fn ext_transition_constraints(
-        _challenges: &RamTableChallenges,
+        challenges: &RamTableChallenges,
     ) -> Vec<MPolynomial<XFieldElement>> {
         use RamBaseTableColumn::*;
 
         let variables: Vec<MPolynomial<XFieldElement>> =
             MPolynomial::variables(2 * FULL_WIDTH, 1.into());
         let one = MPolynomial::from_constant(1.into(), 2 * FULL_WIDTH);
 
+        let bezout_challenge =
+            MPolynomial::from_constant(challenges.bezout_relation_sample_point, FULL_WIDTH);
+
         let clk = variables[usize::from(CLK)].clone();
         let ramp = variables[usize::from(RAMP)].clone();
         let ramv = variables[usize::from(RAMV)].clone();
-        let hv6 = variables[usize::from(InverseOfRampDifference)].clone();
+        let iord = variables[usize::from(InverseOfRampDifference)].clone();
+        let bcpc0 = variables[usize::from(BezoutCoefficientPolynomialCoefficient0)].clone();
+        let bcpc1 = variables[usize::from(BezoutCoefficientPolynomialCoefficient1)].clone();
+        let rp = variables[usize::from(RunningProductOfRAMP)].clone();
+        let fd = variables[usize::from(FormalDerivative)].clone();
+        let bc0 = variables[usize::from(BezoutCoefficient0)].clone();
+        let bc1 = variables[usize::from(BezoutCoefficient1)].clone();
 
         let clk_next = variables[FULL_WIDTH + usize::from(CLK)].clone();
         let ramp_next = variables[FULL_WIDTH + usize::from(RAMP)].clone();
         let ramv_next = variables[FULL_WIDTH + usize::from(RAMV)].clone();
+        let bcpc0_next =
+            variables[FULL_WIDTH + usize::from(BezoutCoefficientPolynomialCoefficient0)].clone();
+        let bcpc1_next =
+            variables[FULL_WIDTH + usize::from(BezoutCoefficientPolynomialCoefficient1)].clone();
+        let rp_next = variables[FULL_WIDTH + usize::from(RunningProductOfRAMP)].clone();
+        let fd_next = variables[FULL_WIDTH + usize::from(FormalDerivative)].clone();
+        let bc0_next = variables[FULL_WIDTH + usize::from(BezoutCoefficient0)].clone();
+        let bc1_next = variables[FULL_WIDTH + usize::from(BezoutCoefficient1)].clone();
+
+        let ramp_diff = ramp_next.clone() - ramp;
+        let ramp_changes = ramp_diff.clone() * iord.clone();
+
+        // iord is 0 or iord is the inverse of (ramp' - ramp)
+        let iord_is_0_or_iord_is_inverse_of_ramp_diff = iord * (ramp_changes.clone() - one.clone());
+
+        // (ramp' - ramp) is zero or iord is the inverse of (ramp' - ramp)
+        let ramp_diff_is_0_or_iord_is_inverse_of_ramp_diff =
+            ramp_diff.clone() * (ramp_changes.clone() - one.clone());
+
+        // The ramp does not change or the new ramv is 0
+        let ramp_does_not_change_or_ramv_becomes_0 = ramp_diff.clone() * ramv_next.clone();
 
-        let ramp_diff = ramp_next - ramp;
+        // The ramp does change or the ramv does not change or the clk increases by 1
+        let ramp_does_not_change_or_ramv_does_not_change_or_clk_increases_by_1 =
+            (ramp_changes.clone() - one.clone())
+                * (ramv_next - ramv)
+                * (clk_next - (clk + one.clone()));
 
-        // hv6 is 0 or hv6 is the inverse of (ramp' - ramp).
-        //
-        // $ hv6·(hv6·(ramp' - ramp) - 1) = 0 $
-        let hv6_is_0_or_hv6_is_inverse_of_ramp_diff =
-            hv6.clone() * (hv6.clone() * ramp_diff.clone() - one.clone());
+        let bcbp0_only_changes_if_ramp_changes =
+            (one.clone() - ramp_changes.clone()) * (bcpc0_next.clone() - bcpc0);
 
-        // (ramp' - ramp) is zero or hv6 is the inverse of (ramp' - ramp).
-        //
-        // $ (ramp' - ramp)·(hv6·(ramp' - ramp) - 1) = 0 $
-        let ramp_diff_is_0_or_hv6_is_inverse_of_ramp_diff =
-            ramp_diff.clone() * (hv6.clone() * ramp_diff.clone() - one.clone());
+        let bcbp1_only_changes_if_ramp_changes =
+            (one.clone() - ramp_changes.clone()) * (bcpc1_next.clone() - bcpc1);
 
-        // The ramp does not change or the new ramv is 0.
-        //
-        // (ramp' - ramp)·ramv'
-        let ramp_does_not_change_or_ramv_becomes_0 = ramp_diff.clone() * ramv_next.clone();
+        let running_product_ramp_updates_correctly = ramp_diff.clone()
+            * (rp_next.clone() - rp.clone() * (bezout_challenge.clone() - ramp_next.clone()))
+            + (one.clone() - ramp_changes.clone()) * (rp_next - rp.clone());
 
-        // The ramp does change or the ramv does not change or the clk increases by 1.
-        //
-        // $ (hv6·(ramp' - ramp) - 1)·(ramv' - ramv)·(clk' - (clk + 1)) = 0 $
-        let ramp_does_not_change_or_ramv_does_not_change_or_clk_increases_by_1 =
-            (hv6 * ramp_diff - one.clone()) * (ramv_next - ramv) * (clk_next - (clk + one));
+        let formal_derivative_updates_correctly = ramp_diff.clone()
+            * (fd_next.clone() - rp - (bezout_challenge.clone() - ramp_next) * fd.clone())
+            + (one.clone() - ramp_changes.clone()) * (fd_next - fd);
+
+        let bezout_coefficient_0_is_constructed_correctly = ramp_diff.clone()
+            * (bc0_next.clone() - bezout_challenge.clone() * bc0.clone() - bcpc0_next)
+            + (one.clone() - ramp_changes.clone()) * (bc0_next - bc0);
+
+        let bezout_coefficient_1_is_constructed_correctly = ramp_diff
+            * (bc1_next.clone() - bezout_challenge * bc1.clone() - bcpc1_next)
+            + (one - ramp_changes) * (bc1_next - bc1);
 
         vec![
-            hv6_is_0_or_hv6_is_inverse_of_ramp_diff,
-            ramp_diff_is_0_or_hv6_is_inverse_of_ramp_diff,
+            iord_is_0_or_iord_is_inverse_of_ramp_diff,
+            ramp_diff_is_0_or_iord_is_inverse_of_ramp_diff,
             ramp_does_not_change_or_ramv_becomes_0,
             ramp_does_not_change_or_ramv_does_not_change_or_clk_increases_by_1,
+            bcbp0_only_changes_if_ramp_changes,
+            bcbp1_only_changes_if_ramp_changes,
+            running_product_ramp_updates_correctly,
+            formal_derivative_updates_correctly,
+            bezout_coefficient_0_is_constructed_correctly,
+            bezout_coefficient_1_is_constructed_correctly,
         ]
     }
 
     fn ext_terminal_constraints(
         _challenges: &RamTableChallenges,
     ) -> Vec<MPolynomial<XFieldElement>> {
-        // no further constraints
-        vec![]
+        let one = MPolynomial::from_constant(1.into(), FULL_WIDTH);
+        let variables = MPolynomial::variables(FULL_WIDTH, 1.into());
+
+        let rp = variables[usize::from(RunningProductOfRAMP)].clone();
+        let fd = variables[usize::from(FormalDerivative)].clone();
+        let bc0 = variables[usize::from(BezoutCoefficient0)].clone();
+        let bc1 = variables[usize::from(BezoutCoefficient1)].clone();
+
+        let bezout_relation_holds = bc0 * rp + bc1 * fd - one;
+
+        vec![bezout_relation_holds]
     }
 }
 
 #[derive(Debug, Clone)]
 pub struct RamTableChallenges {
+    /// The point in which the Bézout relation establishing contiguous memory regions is queried.
+    pub bezout_relation_sample_point: XFieldElement,
+
     /// The weight that combines two consecutive rows in the
     /// permutation/evaluation column of the op-stack table.
     pub processor_perm_row_weight: XFieldElement,