mirror of
https://github.com/logos-storage/plonky2.git
synced 2026-01-04 06:43:07 +00:00
u32 division (#517)
* First draft for division. * `eval_division` work * Division * Minor: outdated fixme * Tests and better column names * Minor lints * Remove redundant constraint * Make division proof more formal * Minor proof and comments Co-authored-by: Hamish Ivey-Law <hamish@ivey-law.name>
This commit is contained in:
Jakub Nabaglo
GitHub
parent
2cedd1b02a
commit
06fef55bfb
93
system_zero/src/alu/div-constraints-proof.txt
Normal file
93
system_zero/src/alu/div-constraints-proof.txt
Normal file
@ -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.
|
||||
@ -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<F: PrimeField64>(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<F: Field, P: PackedField<Scalar = F>>(
|
||||
local_values: &[P; NUM_COLUMNS],
|
||||
lv: &[P; NUM_COLUMNS],
|
||||
yield_constr: &mut ConstraintConsumer<P>,
|
||||
) {
|
||||
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<F: RichField + Extendable<D>, const D: usize>(
|
||||
builder: &mut CircuitBuilder<F, D>,
|
||||
local_values: &[ExtensionTarget<D>; NUM_COLUMNS],
|
||||
lv: &[ExtensionTarget<D>; NUM_COLUMNS],
|
||||
yield_constr: &mut RecursiveConstraintConsumer<F, D>,
|
||||
) {
|
||||
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::<u32>());
|
||||
values[COL_DIV_INPUT_DIVISOR] = F::from_canonical_u32(rng.gen::<u32>());
|
||||
|
||||
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.
|
||||
}
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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 {
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user