diff --git a/system_zero/src/alu/div-constraints-proof.txt b/system_zero/src/alu/div-constraints-proof.txt new file mode 100644 index 00000000..bd0b687c --- /dev/null +++ b/system_zero/src/alu/div-constraints-proof.txt @@ -0,0 +1,93 @@ +Constraints A (implemented in code): +A1. dividend ∈ {0, …, u32::MAX} +A2. divisor ∈ {0, …, u32::MAX} +A3. quotient ∈ {0, …, u32::MAX} +A4. remainder ∈ {0, …, u32::MAX} +A5. divisor_rem_diff_m1 ∈ {0, …, u32::MAX} +A6. divisor * div_inverse = div_div_inverse +A7. (div_div_inverse - 1) * (remainder - quotient - u32::MAX) = 0 +A8. divisor * (div_div_inverse - 1) = 0 +A9. div_inverse * dividend = quotient + remainder * div_inverse +A10. divisor * (divisor - remainder - 1 - divisor_rem_diff_m1) = 0 + +Constraints B (intuitive division): +B1. dividend ∈ {0, …, u32::MAX} +B2. divisor ∈ {0, …, u32::MAX} +B3. divisor = 0 => quotient = 0 +B4. divisor = 0 => remainder = u32::MAX +B5. divisor ≠ 0 => dividend = quotient * divisor + remainder +B6. divisor ≠ 0 => quotient ∈ {0, …, u32::MAX} +B7. divisor ≠ 0 => remainder ∈ {0, …, divisor - 1} + + + +Assume we meet constraints A for some dividend, divisor, quotient, remainder, divisor_rem_diff_m1, div_inverse, and div_div_inverse. We want to show that constrants B are met. + +B1. Trivial by A1. + +B2. Trivial by A2. + +B3. Assume divisor = 0. Then div_div_inverse = 0 by A6. div_div_inverse - 1 ≠ 0, so remainder - quotient = u32::MAX by A7. +quotient ∈ {0, …, u32::MAX} by A3 and remainder ∈ {0, …, u32::MAX} by A4. Then remainder - quotient ∈ {-quotient, …, u32::MAX - quotient}. +If quotient ≠ 0, then quotient ∈ {1, …, u32::MAX} and remainder - quotient ∈ {-u32::MAX, …, u32::MAX - 1}, which does not include u32::MAX, contradicting A7. + +B4. Assume divisor = 0. Then div_div_inverse = 0 by A6. div_div_inverse - 1 ≠ 0, so remainder - quotient = u32::MAX by A7. +quotient ∈ {0, …, u32::MAX} by A3 and remainder ∈ {0, …, u32::MAX} by A4. Then remainder - quotient ∈ {remainder - u32::MAX, …, remainder}. +If remainder ≠ u32::MAX, then remainder ∈ {0, …, u32::MAX - 1} and remainder - quotient ∈ {-u32::MAX, …, u32::MAX - 1} which does not include u32::MAX, contradicting A7. + +B5. Assume divisor ≠ 0. By A8, div_div_inverse = 1. By A6, div_inverse = divisor^-1. Multiplying both sides of A9 by divisor, dividend = quotient * divisor + remainder. + +B6. Follows from A3. + +B7. remainder ∈ {0, …, u32::MAX} by A4. Assume divisor ≠ 0. Then divisor_rem_diff_m1 = divisor - remainder - 1 by A10. divisor ∈ {1, …, u32::MAX} by A2. If remainder ∈ {divisor, …, u32::MAX}, then divisor - remainder - 1 ∈ {-u32::MAX, …, u32::MAX - divisor} ⊆ {-u32::MAX, …, u32::MAX - 1}, contradicting A5. Hence, remainder ∈ {0, …, divisor - 1}. + + + +Assume we meet constraints B for some dividend, divisor, quotient, and remainder. We want to show +that there exist divisor_rem_diff_m1, div_inverse, div_div_inverse, such that constrants A are met. + +If divisor = 0, set divisor_rem_diff_m1 = 0, div_inverse = 0, div_div_inverse = 0. +Otherwise, set divisor_rem_diff_m1 = divisor - remainder - 1, div_inverse = divisor^-1, div_div_inverse = 1. + +A1. Trivial by B1. + +A2. Trivial by B2. + +The remainder is by cases: + +(divisor = 0) + + A3. Follows from B3. + + A4. Follows from B4. + + A5. Follows from our choice of divisor_rem_diff_m1 = 0. + + A6. Follows from our choice of div_div_inverse = 0. + + A7. quotient = 0 by B3. remainder = u32::MAX by B4. Then remainder - quotient = u32::MAX. + + A8. Trivial since divisor = 0. + + A9. By our choice, div_inverse = 0. quotient = 0 by B3. + + A10. Trivial since divisor = 0. + + +(divisor ≠ 0) + + A3. Follows from B6. + + A4. By B7, remainder ∈ {0, …, divisor - 1}, and by B2, divisor ∈ {0, …, u32::MAX}, implying that remainder ∈ {0, …, u32::MAX - 1}. + + A5. We've set divisor_rem_diff_m1 = divisor - remainder - 1. remainder ∈ {0, …, divisor - 1}, so divisor - remainder ∈ {1, …, divisor} and divisor - remainder - 1 = divisor_rem_diff_m1 ∈ {0, …, divisor - 1}. From B2, divisor ∈ {0, …, u32::MAX}, so divisor_rem_diff_m1 ∈ {0, …, u32::MAX - 1} as desired. + + A6. div_inverse = divisor^-1 by choice, so divisor * div_inverse = 1. div_div_inverse = 1 by choice. + + A7. div_div_inverse = 1 by choice, so div_div_inverse - 1 = 0. + + A8. div_div_inverse = 1 by choice, so div_div_inverse - 1 = 0. + + A9. From B5, dividend = quotient * divisor + remainder. Since divisor ≠ 0, div_inverse = divisor^-1 by choice. Multiplying both sides by div_inverse, dividend * div_inverse = quotient * divisor * div_inverse + remainder * div_inverse = quotient + remainder * div_inverse. + + A10. By our choice of divisor_rem_diff_m1 = divisor - remainder - 1. diff --git a/system_zero/src/alu/division.rs b/system_zero/src/alu/division.rs index f0d645f1..eda44442 100644 --- a/system_zero/src/alu/division.rs +++ b/system_zero/src/alu/division.rs @@ -9,23 +9,215 @@ use starky::constraint_consumer::{ConstraintConsumer, RecursiveConstraintConsume use crate::registers::alu::*; use crate::registers::NUM_COLUMNS; +/// Division instruction of a u32 divisor N into a u32 dividend D, +/// with u32 quotient Q and u32 remainder R. If D is not zero, then +/// the values will satisfy N = Q*D + R with 0 <= R < D. If D is +/// zero, then the remainder is set to the special value u32::MAX = +/// 2^32 - 1 (which is not a valid remainder for any nonzero D) and +/// the quotient is set to zero. In particular, no overflow is +/// possible. + pub(crate) fn generate_division(values: &mut [F; NUM_COLUMNS]) { - // TODO + let dividend = values[COL_DIV_INPUT_DIVIDEND].to_canonical_u64() as u32; + let divisor = values[COL_DIV_INPUT_DIVISOR].to_canonical_u64() as u32; + + // `COL_DIV_INVDIVISOR` is `divisor^-1` if `divisor != 0` and `0` otherwise. + // `COL_DIV_NONZERO_DIVISOR` is `1` if `divisor != 0` and `0` otherwise. + + // `COL_DIV_RANGE_CHECKED_TMP` is set to `divisor - rem - 1` if `divisor != 0` and `0` + // otherwise. This is used to ensure that `rem < divisor` when `divisor != 0`. + + if divisor == 0 { + // Outputs + values[COL_DIV_OUTPUT_QUOT_0] = F::ZERO; + values[COL_DIV_OUTPUT_QUOT_1] = F::ZERO; + values[COL_DIV_OUTPUT_REM_0] = F::from_canonical_u16(u16::MAX); + values[COL_DIV_OUTPUT_REM_1] = F::from_canonical_u16(u16::MAX); + + // Temporaries + values[COL_DIV_RANGE_CHECKED_TMP_0] = F::ZERO; + values[COL_DIV_RANGE_CHECKED_TMP_1] = F::ZERO; + values[COL_DIV_INVDIVISOR] = F::ZERO; + values[COL_DIV_NONZERO_DIVISOR] = F::ZERO; + } else { + let quo = dividend / divisor; + let rem = dividend % divisor; + + let div_rem_diff_m1 = divisor - rem - 1; + + // Outputs + values[COL_DIV_OUTPUT_QUOT_0] = F::from_canonical_u16(quo as u16); + values[COL_DIV_OUTPUT_QUOT_1] = F::from_canonical_u16((quo >> 16) as u16); + values[COL_DIV_OUTPUT_REM_0] = F::from_canonical_u16(rem as u16); + values[COL_DIV_OUTPUT_REM_1] = F::from_canonical_u16((rem >> 16) as u16); + + // Temporaries + values[COL_DIV_RANGE_CHECKED_TMP_0] = F::from_canonical_u16(div_rem_diff_m1 as u16); + values[COL_DIV_RANGE_CHECKED_TMP_1] = F::from_canonical_u16((div_rem_diff_m1 >> 16) as u16); + values[COL_DIV_INVDIVISOR] = F::from_canonical_u32(divisor).inverse(); + values[COL_DIV_NONZERO_DIVISOR] = F::ONE; + } } pub(crate) fn eval_division>( - local_values: &[P; NUM_COLUMNS], + lv: &[P; NUM_COLUMNS], yield_constr: &mut ConstraintConsumer

, ) { - let is_div = local_values[IS_DIV]; - // TODO + let base = F::from_canonical_u64(1 << 16); + let u32_max = P::from(F::from_canonical_u32(u32::MAX)); + + // Filter + let is_div = lv[IS_DIV]; + + // Inputs + let dividend = lv[COL_DIV_INPUT_DIVIDEND]; + let divisor = lv[COL_DIV_INPUT_DIVISOR]; + + // Outputs + let quotient = lv[COL_DIV_OUTPUT_QUOT_0] + lv[COL_DIV_OUTPUT_QUOT_1] * base; + let remainder = lv[COL_DIV_OUTPUT_REM_0] + lv[COL_DIV_OUTPUT_REM_1] * base; + + // Temporaries + let divinv = lv[COL_DIV_INVDIVISOR]; + let div_divinv = lv[COL_DIV_NONZERO_DIVISOR]; + let div_rem_diff_m1 = lv[COL_DIV_RANGE_CHECKED_TMP_0] + lv[COL_DIV_RANGE_CHECKED_TMP_1] * base; + + // Constraints + yield_constr.constraint(is_div * (divisor * divinv - div_divinv)); + yield_constr.constraint(is_div * (div_divinv - F::ONE) * (remainder - quotient - u32_max)); + yield_constr.constraint(is_div * divisor * (div_divinv - F::ONE)); + yield_constr.constraint(is_div * (quotient + remainder * divinv - divinv * dividend)); + yield_constr.constraint(is_div * divisor * (divisor - remainder - F::ONE - div_rem_diff_m1)); } pub(crate) fn eval_division_recursively, const D: usize>( builder: &mut CircuitBuilder, - local_values: &[ExtensionTarget; NUM_COLUMNS], + lv: &[ExtensionTarget; NUM_COLUMNS], yield_constr: &mut RecursiveConstraintConsumer, ) { - let is_div = local_values[IS_DIV]; - // TODO + let base = builder.constant_extension(F::Extension::from_canonical_u64(1 << 16)); + let u32_max = builder.constant_extension(F::Extension::from_canonical_u32(u32::MAX)); + let one = builder.constant_extension(F::Extension::ONE); + + // Filter + let is_div = lv[IS_DIV]; + + // Inputs + let dividend = lv[COL_DIV_INPUT_DIVIDEND]; + let divisor = lv[COL_DIV_INPUT_DIVISOR]; + + // Outputs + let quotient = + builder.mul_add_extension(lv[COL_DIV_OUTPUT_QUOT_1], base, lv[COL_DIV_OUTPUT_QUOT_0]); + let remainder = + builder.mul_add_extension(lv[COL_DIV_OUTPUT_REM_1], base, lv[COL_DIV_OUTPUT_REM_0]); + + // Temporaries + let divinv = lv[COL_DIV_INVDIVISOR]; + let div_divinv = lv[COL_DIV_NONZERO_DIVISOR]; + let div_rem_diff_m1 = builder.mul_add_extension( + lv[COL_DIV_RANGE_CHECKED_TMP_1], + base, + lv[COL_DIV_RANGE_CHECKED_TMP_0], + ); + + // Constraints + let constr6 = builder.mul_sub_extension(divisor, divinv, div_divinv); + let constr7 = { + let t = builder.sub_extension(div_divinv, one); + let u = builder.sub_extension(remainder, quotient); + let v = builder.sub_extension(u, u32_max); + builder.mul_extension(t, v) + }; + let constr8 = { + let t = builder.sub_extension(div_divinv, one); + builder.mul_extension(divisor, t) + }; + let constr9 = { + let t = builder.sub_extension(remainder, dividend); + builder.mul_add_extension(t, divinv, quotient) + }; + let constr10 = { + let t = builder.sub_extension(divisor, remainder); + let u = builder.add_extension(one, div_rem_diff_m1); + let v = builder.sub_extension(t, u); + builder.mul_extension(divisor, v) + }; + + let constr6 = builder.mul_extension(is_div, constr6); + let constr7 = builder.mul_extension(is_div, constr7); + let constr8 = builder.mul_extension(is_div, constr8); + let constr9 = builder.mul_extension(is_div, constr9); + let constr10 = builder.mul_extension(is_div, constr10); + + yield_constr.constraint(builder, constr6); + yield_constr.constraint(builder, constr7); + yield_constr.constraint(builder, constr8); + yield_constr.constraint(builder, constr9); + yield_constr.constraint(builder, constr10); +} + +#[cfg(test)] +mod tests { + use plonky2::field::field_types::Field; + use plonky2::field::goldilocks_field::GoldilocksField; + use rand::{Rng, SeedableRng}; + use rand_chacha::ChaCha8Rng; + use starky::constraint_consumer::ConstraintConsumer; + + use super::*; + use crate::registers::NUM_COLUMNS; + + #[test] + fn generate_eval_consistency_not_div() { + const D: usize = 1; + type F = GoldilocksField; + + let mut rng = ChaCha8Rng::seed_from_u64(0x6feb51b7ec230f25); + let mut values = [F::default(); NUM_COLUMNS].map(|_| F::rand_from_rng(&mut rng)); + + // if `IS_DIV == 0`, then the constraints should be met even if all values are garbage. + values[IS_DIV] = F::ZERO; + + let mut constrant_consumer = ConstraintConsumer::new( + vec![GoldilocksField(2), GoldilocksField(3), GoldilocksField(5)], + GoldilocksField::ONE, + GoldilocksField::ONE, + GoldilocksField::ONE, + ); + eval_division(&values, &mut constrant_consumer); + for &acc in &constrant_consumer.constraint_accs { + assert_eq!(acc, GoldilocksField::ZERO); + } + } + + #[test] + fn generate_eval_consistency_div() { + const D: usize = 1; + type F = GoldilocksField; + + let mut rng = ChaCha8Rng::seed_from_u64(0x6feb51b7ec230f25); + let mut values = [F::default(); NUM_COLUMNS].map(|_| F::rand_from_rng(&mut rng)); + + // set `IS_DIV == 1` and ensure all constraints are satisfied. + values[IS_DIV] = F::ONE; + // set `DIVIDEND` and `DIVISOR` to `u32`s + values[COL_DIV_INPUT_DIVIDEND] = F::from_canonical_u32(rng.gen::()); + values[COL_DIV_INPUT_DIVISOR] = F::from_canonical_u32(rng.gen::()); + + generate_division(&mut values); + + let mut constrant_consumer = ConstraintConsumer::new( + vec![GoldilocksField(2), GoldilocksField(3), GoldilocksField(5)], + GoldilocksField::ONE, + GoldilocksField::ONE, + GoldilocksField::ONE, + ); + eval_division(&values, &mut constrant_consumer); + for &acc in &constrant_consumer.constraint_accs { + assert_eq!(acc, GoldilocksField::ZERO); + } + } + + // TODO: test eval_division_recursively. } diff --git a/system_zero/src/registers/alu.rs b/system_zero/src/registers/alu.rs index 585ecab1..7b89b479 100644 --- a/system_zero/src/registers/alu.rs +++ b/system_zero/src/registers/alu.rs @@ -66,4 +66,27 @@ pub(crate) const COL_MUL_ADD_OUTPUT_2: usize = super::range_check_16::col_rc_16_ /// The fourth 16-bit chunk of the output, based on little-endian ordering. pub(crate) const COL_MUL_ADD_OUTPUT_3: usize = super::range_check_16::col_rc_16_input(3); -pub(super) const END: usize = super::START_ALU + NUM_SHARED_COLS; +/// Dividend for division, as an unsigned u32 +pub(crate) const COL_DIV_INPUT_DIVIDEND: usize = shared_col(0); +/// Divisor for division, as an unsigned u32 +pub(crate) const COL_DIV_INPUT_DIVISOR: usize = shared_col(1); +/// Inverse of divisor, if one exists, and 0 otherwise +pub(crate) const COL_DIV_INVDIVISOR: usize = shared_col(2); +/// 1 if divisor is nonzero and 0 otherwise +pub(crate) const COL_DIV_NONZERO_DIVISOR: usize = shared_col(3); + +/// The first 16-bit chunk of the quotient, based on little-endian ordering. +pub(crate) const COL_DIV_OUTPUT_QUOT_0: usize = super::range_check_16::col_rc_16_input(0); +/// The second 16-bit chunk of the quotient, based on little-endian ordering. +pub(crate) const COL_DIV_OUTPUT_QUOT_1: usize = super::range_check_16::col_rc_16_input(1); +/// The first 16-bit chunk of the remainder, based on little-endian ordering. +pub(crate) const COL_DIV_OUTPUT_REM_0: usize = super::range_check_16::col_rc_16_input(2); +/// The second 16-bit chunk of the remainder, based on little-endian ordering. +pub(crate) const COL_DIV_OUTPUT_REM_1: usize = super::range_check_16::col_rc_16_input(3); + +/// The first 16-bit chunk of a temporary value (divisor - remainder - 1). +pub(crate) const COL_DIV_RANGE_CHECKED_TMP_0: usize = super::range_check_16::col_rc_16_input(4); +/// The second 16-bit chunk of a temporary value (divisor - remainder - 1). +pub(crate) const COL_DIV_RANGE_CHECKED_TMP_1: usize = super::range_check_16::col_rc_16_input(5); + +pub(super) const END: usize = START_SHARED_COLS + NUM_SHARED_COLS; diff --git a/system_zero/src/registers/range_check_16.rs b/system_zero/src/registers/range_check_16.rs index 674df302..cdb11065 100644 --- a/system_zero/src/registers/range_check_16.rs +++ b/system_zero/src/registers/range_check_16.rs @@ -1,6 +1,6 @@ //! Range check unit which checks that values are in `[0, 2^16)`. -pub(crate) const NUM_RANGE_CHECKS: usize = 5; +pub(super) const NUM_RANGE_CHECKS: usize = 6; /// The input of the `i`th range check, i.e. the value being range checked. pub(crate) const fn col_rc_16_input(i: usize) -> usize {