fix: poly_reduce and poly_divide_by_cyclo

This commit is contained in:
Enrico Bottazzi
2023-10-05 10:05:24 +02:00
parent 61a1e8e81c
commit e4ad7a6657
3 changed files with 138 additions and 184 deletions

View File

@@ -1,2 +1,3 @@
pub mod poly_distribution;
pub mod poly_operations;
pub mod utils;

View File

@@ -1,3 +1,4 @@
use crate::chips::utils::{div_euclid, vec_assigned_to_vec_u64};
use halo2_base::gates::GateChip;
use halo2_base::gates::GateInstructions;
use halo2_base::safe_types::RangeChip;
@@ -6,12 +7,11 @@ use halo2_base::utils::ScalarField;
use halo2_base::AssignedValue;
use halo2_base::Context;
use halo2_base::QuantumCell;
use zk_fhe::chips::utils::div_euclid;
/// Build the sum of the polynomials a and b as sum of the coefficients
/// DEG is the degree of the input polynomials
/// Input polynomials are parsed as a vector of assigned coefficients [a_DEG, a_DEG-1, ..., a_1, a_0] where a_0 is the constant term
pub fn poly_add<const DEG: u64, F: ScalarField>(
pub fn poly_add<const DEG: usize, F: ScalarField>(
ctx: &mut Context<F>,
a: Vec<AssignedValue<F>>,
b: Vec<AssignedValue<F>>,
@@ -19,17 +19,17 @@ pub fn poly_add<const DEG: u64, F: ScalarField>(
) -> Vec<AssignedValue<F>> {
// assert that the input polynomials have the same degree and this is equal to DEG
assert_eq!(a.len() - 1, b.len() - 1);
assert_eq!(a.len() - 1, DEG as usize);
assert_eq!(a.len() - 1, DEG);
let c: Vec<AssignedValue<F>> = a
.iter()
.zip(b.iter())
.take(2 * (DEG as usize) - 1)
.take(2 * (DEG) - 1)
.map(|(&a, &b)| gate.add(ctx, a, b))
.collect();
// assert that the sum polynomial has degree DEG
assert_eq!(c.len() - 1, DEG as usize);
assert_eq!(c.len() - 1, DEG);
c
}
@@ -123,7 +123,7 @@ pub fn poly_mul_diff_deg<F: ScalarField>(
/// Build the scalar multiplication of the polynomials a and the scalar k as scalar multiplication of the coefficients of a and k
/// DEG is the degree of the polynomial
/// Input polynomials is parsed as a vector of assigned coefficients [a_DEG, a_DEG-1, ..., a_1, a_0] where a_0 is the constant term
/// Input polynomial is parsed as a vector of assigned coefficients [a_DEG, a_DEG-1, ..., a_1, a_0] where a_0 is the constant term
pub fn poly_scalar_mul<const DEG: u64, F: ScalarField>(
ctx: &mut Context<F>,
a: Vec<AssignedValue<F>>,
@@ -141,111 +141,96 @@ pub fn poly_scalar_mul<const DEG: u64, F: ScalarField>(
c
}
// Takes a polynomial represented by its coefficients in a vector and output a new polynomial reduced by applying modulo Q
// Q is the modulus
// N is the degree of the polynomials
pub fn poly_reduce<const Q: u64, const N: u64, F: ScalarField>(
/// Takes a polynomial represented by its coefficients in a vector and output a new polynomial reduced by applying modulo Q to each coefficient
/// DEG is the degree of the polynomial
/// Input polynomial is parsed as a vector of assigned coefficients [a_DEG, a_DEG-1, ..., a_1, a_0] where a_0 is the constant term
/// It assumes that the coefficients of the input polynomial can be expressed in at most num_bits bits
pub fn poly_reduce<const Q: u64, const DEG: usize, F: ScalarField>(
ctx: &mut Context<F>,
input: Vec<AssignedValue<F>>,
range: &RangeChip<F>,
num_bits: usize,
) -> Vec<AssignedValue<F>> {
// Assert that degree is equal to the constant N
assert_eq!(input.len() - 1, N as usize);
// Assign the input polynomials to the circuit
let in_assigned: Vec<AssignedValue<F>> = input
.iter()
.map(|x| {
let result = F::from(*x as u64);
ctx.load_witness(result)})
.collect();
// Assert that degree of input polynomial is equal to the constant DEG
assert_eq!(input.len() - 1, DEG);
// Enforce that in_assigned[i] % Q = rem_assigned[i]
// coefficients of input polynomials are guaranteed to be at most 16 bits by assumption
let rem_assigned: Vec<AssignedValue<F>> =
in_assigned.iter().take(2 * N - 1).map(|&x| range.div_mod(ctx, x, Q, 16).1).collect();
let rem_assigned: Vec<AssignedValue<F>> = input
.iter()
.take(2 * (DEG) - 1)
.map(|&x| range.div_mod(ctx, x, Q, num_bits).1)
.collect();
// assert that the reduced polynomial has degree N
assert_eq!(rem_assigned.len() - 1, N as usize);
// assert that the reduced polynomial has degree DEG
assert_eq!(rem_assigned.len() - 1, DEG);
rem_assigned
}
// Takes a polynomial represented by its coefficients in a vector
// and output a remainder polynomial divided by divisor polynomial f(x)=x^m+1 where m is a power of 2 (public output)
// N is the degree of the dividend polynomial
// M is the degree of the divisor polynomial
// Q is the modulus
pub fn poly_divide_by_cyclo<const N: u64, const M: u64, const Q: u64, F: ScalarField>(
/// Takes a polynomial `divisor` represented by its coefficients in a vector
/// Takes a cyclotomic polynomial `dividend` f(x)=x^m+1 (m is a power of 2) of the form represented by its coefficients in a vector
/// Output the remainder of the division of `dividend` by `dividend` as a vector of coefficients
/// DEG_DVD is the degree of the `dividend` polynomial
/// DEG_DVS is the degree of the `divisor` polynomial
/// Q is the modulus of the Ring
/// Input polynomials is parsed as a vector of assigned coefficients [a_DEG, a_DEG-1, ..., a_1, a_0] where a_0 is the constant term
/// Assumes that the coefficients of `dividend` are in the range [0, Q - 1]
/// Assumes that divisor is a cyclotomic polynomial with coefficients either 0 or 1
pub fn poly_divide_by_cyclo<
const DEG_DVD: usize,
const DEG_DVS: usize,
const Q: u64,
F: ScalarField,
>(
ctx: &mut Context<F>,
dividend: Vec<AssignedValue<F>>,
divisor: Vec<AssignedValue<F>>,
range: &RangeChip<F>,
) -> Vec<AssignedValue<F>> {
// Assert that degree of dividend poly is equal to the constant N
assert_eq!(dividend.len() - 1, N as usize);
// Assert that degree of divisor poly is equal to the constant M
assert_eq!(divisor.len() - 1, M as usize);
// Assert that degree of dividend polynomial is equal to the constant DEG_DVD
assert_eq!(dividend.len() - 1, DEG_DVD);
// Assert that degree of divisor poly is equal to the constant DEG_DVS
assert_eq!(divisor.len() - 1, DEG_DVS);
// M must be less than N
assert!(M < N);
let mut dividend_assigned = vec![];
let mut divisor_assigned = vec![];
for i in 0..N + 1 {
let val = F::from(dividend[i]);
let assigned_val = ctx.load_witness(val);
dividend_assigned.push(assigned_val);
}
for i in 0..M + 1 {
let val = F::from(divisor[i]);
let assigned_val = ctx.load_witness(val);
divisor_assigned.push(assigned_val);
}
// assert the correct length of the assigned polynomails
assert_eq!(dividend_assigned.len() - 1, N as usize);
assert_eq!(divisor_assigned.len() - 1, M as usize);
// DEG_DVS must be strictly less than DEG_DVD
assert!(DEG_DVS < DEG_DVD);
// long division operation performed outside the circuit
let (quotient, remainder) = div_euclid(&dividend, &divisor);
// Need to convert the dividend and divisor into a vector of u64
let dividend_to_u64 = vec_assigned_to_vec_u64(&dividend);
let divisor_to_u64 = vec_assigned_to_vec_u64(&divisor);
// After the division, the degree of the quotient is N - M
assert_eq!(quotient.len() - 1, N - M);
let (quotient_to_u64, remainder_to_u64) = div_euclid::<Q>(&dividend_to_u64, &divisor_to_u64);
// Furthermore, the degree of the remainder must be strictly less than the degree of the divisor which is M
assert!(remainder.len() - 1 < M as usize);
// After the division, the degree of the quotient should be equal to DEG_DVD - DEG_DVS
assert_eq!(quotient_to_u64.len() - 1, DEG_DVD - DEG_DVS);
// Pad the remainder with 0s to make its degree equal to M - 1
let mut remainder = remainder;
while remainder.len() - 1 < M - 1 {
remainder.push(0);
// Furthermore, the degree of the remainder must be strictly less than the degree of the divisor
assert!(remainder_to_u64.len() - 1 < DEG_DVS);
// Pad the remainder with 0s to make its degree equal to DEG_DVS - 1
let mut remainder_to_u64 = remainder_to_u64;
while remainder_to_u64.len() - 1 < DEG_DVS - 1 {
remainder_to_u64.push(0);
}
// Now remainder must be of degree M - 1
assert_eq!(remainder.len() - 1, M - 1);
// Now remainder must be of degree DEG_DVS - 1
assert_eq!(remainder_to_u64.len() - 1, DEG_DVS - 1);
// Assign the quotient and remainder to the circuit
let mut quot_assigned = vec![];
let mut remainder_assigned = vec![];
let mut quotient = vec![];
let mut remainder = vec![];
for i in 0..N - M + 1 {
let val = F::from(quotient[i]);
for i in 0..DEG_DVD - DEG_DVS + 1 {
let val = F::from(quotient_to_u64[i]);
let assigned_val = ctx.load_witness(val);
quot_assigned.push(assigned_val);
quotient.push(assigned_val);
}
for i in 0..M {
let val = F::from(remainder[i]);
for i in 0..DEG_DVS {
let val = F::from(remainder_to_u64[i]);
let assigned_val = ctx.load_witness(val);
remainder_assigned.push(assigned_val);
remainder.push(assigned_val);
}
// Quotient is obtained by dividing the coefficients of the dividend by the highest degree coefficient of divisor
@@ -253,136 +238,91 @@ pub fn poly_divide_by_cyclo<const N: u64, const M: u64, const Q: u64, F: ScalarF
// The leading coefficient of divisor is 1 by assumption.
// Therefore, the coefficients of quotient have to be in the range [0, Q - 1]
// Since the quotient is computed outside the circuit, we need to enforce this constraint
for i in 0..(N - M + 1) as usize {
range.check_less_than_safe(ctx, quot_assigned[i], Q);
for i in 0..(DEG_DVD - DEG_DVS + 1) {
range.check_less_than_safe(ctx, quotient[i], Q);
}
// Remainder is equal to dividend - (quotient * divisor).
// The coefficients of dividend are in the range [0, Q - 1] by assumption.
// The coefficients of quotient are in the range [0, Q - 1] by constraint set above.
// The coefficients of divisior are either 0, 1 by assumption.
// The coefficients of divisior are either 0, 1 by assumption of the cyclotomic polynomial.
// It follows that the coefficients of quotient * divisor are in the range [0, Q - 1]
// The remainder might have coefficients that are negative. In that case we add Q to them to make them positive.
// Therefore, the coefficients of remainder are in the range [0, Q - 1]
// Since the remainder is computed outside the circuit, we need to enforce this constraint
for i in 0..M as usize {
range.check_less_than_safe(ctx, remainder_assigned[i], Q);
for i in 0..DEG_DVS {
range.check_less_than_safe(ctx, remainder[i], Q);
}
// ---- constraint check -----
// check that quotient * divisor + rem = dividend
// check that quotient * divisor + remainder = dividend
// DEGREE ANALYSIS
// Quotient is of degree N - M
// Divisor is of degree M
// Quotient * divisor is of degree N
// Rem is of degree M - 1
// Quotient * divisor + rem is of degree N
// Dividend is of degree N
// Quotient is of degree DEG_DVD - DEG_DVS
// Divisor is of degree DEG_DVS
// Quotient * divisor is of degree DEG_DVD
// Remainder is of degree DEG_DVS - 1
// Quotient * divisor + rem is of degree DEG_DVD
// Dividend is of degree DEG_DVD
// Perform the multiplication between quotient and divisor
// We use a polynomial multiplication algorithm that does not require the input polynomials to be of the same degree
// Initialize the resulting polynomial prod
// The degree of prod is N
let mut prod = vec![];
// degree of quotient
let a = (N - M) as usize;
// degree of divisor
let b = M as usize;
// No risk of overflowing the prime here when performing multiplication across coefficients
// No risk of overflowing the circuit prime here when performing multiplication across coefficients
// The coefficients of quotient are in the range [0, Q - 1] by constraint set above.
// The coefficients of divisor are either 0, 1 by assumption.
// Calculate the coefficients of the product polynomial
for i in 0..=a + b {
// Init the accumulator for the i-th coefficient of c
let mut coefficient_accumulator = vec![];
// We use a polynomial multiplication algorithm that does not require the input polynomials to be of the same degree
for j in 0..=i {
if j <= a && (i - j) <= b {
let quotient_coef = quot_assigned[j];
let divisor_coef = divisor_assigned[i - j];
// Update the accumulator
coefficient_accumulator.push(range.gate().mul(ctx, quotient_coef, divisor_coef));
}
}
let val = coefficient_accumulator
.iter()
.fold(ctx.load_witness(F::zero()), |acc, x| range.gate().add(ctx, acc, *x));
// Assign the accumulated value to the i-th coefficient of c
prod.push(val);
}
// assert that the degree of prod is N
assert_eq!(prod.len() - 1, N as usize);
let prod = poly_mul_diff_deg(ctx, quotient, divisor, range.gate());
// The coefficients of Divisor are either 0 or 1 given that this is a cyclotomic polynomial
// The coefficients of Quotient are in the range [0, Q - 1] as per constraint set above.
// Therefore, the coefficients of prod are in the range [0, Q - 1]
// Perform the addition between prod and remainder_assigned
// The degree of prod is N
// The degree of remainder_assigned is M - 1
// The degree of prod is DEG_DVD
assert_eq!(prod.len() - 1, DEG_DVD);
// We can pad the remainder_assigned with 0s to make its degree equal to N
let mut remainder_assigned = remainder_assigned;
while remainder_assigned.len() - 1 < N as usize {
remainder_assigned.push(ctx.load_witness(F::zero()));
// Perform the addition between prod and remainder
// The degree of prod is DEG_DVD
// The degree of remainder is DEG_DVS - 1
// We can pad the remainder with 0s to make its degree equal to DEG_DVD
let mut remainder = remainder;
while remainder.len() - 1 < DEG_DVD {
remainder.push(ctx.load_witness(F::zero()));
}
// Now remainder_assigned is of degree N
assert_eq!(remainder_assigned.len() - 1, N as usize);
// Now remainder is of degree DEG_DVD
assert_eq!(remainder.len() - 1, DEG_DVD);
// prod + rem_assigned
// prod + remainder
// No risk of overflowing the prime here when performing addition across coefficients
// No risk of overflowing the circuit prime here when performing addition across coefficients
// The coefficients of prod are in the range [0, Q - 1] by the operation above.
// The coefficients of rem_assigned are in the range [0, Q - 1] by constraint set above.
// The coefficients of remainder are in the range [0, Q - 1] by constraint set above.
let sum: Vec<AssignedValue<F>> = prod
.iter()
.zip(remainder_assigned.iter())
.take(2 * N as usize - 1)
.map(|(&a, &b)| range.gate().add(ctx, a, b))
.collect();
let sum = poly_add::<DEG_DVD, F>(ctx, prod, remainder.clone(), range.gate());
// assert that the degree of sum is N
assert_eq!(sum.len() - 1, N as usize);
// assert that the degree of sum is DEG_DVD
assert_eq!(sum.len() - 1, DEG_DVD);
// The coefficients of prod are in the range [0, Q - 1]
// The coefficients of rem_assigned are in the range [0, Q - 1] by constraint set above.
// The coefficients of remainder are in the range [0, Q - 1] by constraint set above.
// Therefore, the coefficients of sum are in the range [0, 2 * (Q - 1)].
// We can reduce the coefficients of sum modulo Q to make them in the range [0, Q - 1]
// Reduce the values of sum modulo Q
let mut sum_mod = vec![];
// get the number of bits needed to represent the value of 2 * (Q - 1)
let binary_representation = format!("{:b}", (2 * (Q - 1))); // Convert to binary (base-2)
let num_bits = binary_representation.len();
// `div_mod` assumes that sum[i] lives in q_bits which is true given the above analysis
for i in 0..=N as usize {
let assigned_val = range.div_mod(ctx, sum[i], Q, num_bits).1;
sum_mod.push(assigned_val);
let sum_mod = poly_reduce::<Q, DEG_DVD, F>(ctx, sum, range, num_bits);
// assert that the degree of sum_mod is DEG_DVD
assert_eq!(sum_mod.len() - 1, DEG_DVD);
// Enforce that sum_mod = dividend
for i in 0..=DEG_DVD {
let bool = range.gate().is_equal(ctx, sum_mod[i], dividend[i]);
range.gate().assert_is_const(ctx, &bool, &F::from(1))
}
// assert that the degree of sum_mod is N
assert_eq!(sum_mod.len() - 1, N as usize);
// Enforce that sum_mod = dividend_assigned
for i in 0..=N as usize{
range.gate().is_equal(ctx, sum_mod[i], dividend_assigned[i]);
}
// ---- constraint check -----
remainder_assigned
remainder
}

View File

@@ -1,18 +1,17 @@
use halo2_base::{utils::ScalarField, AssignedValue};
pub fn div_euclid(dividend: &Vec<u64>, divisor: &Vec<u64>) -> (Vec<u64>, Vec<u64>) {
/// Performs long polynomial division on two polynomials
/// Returns the quotient and remainder
/// Input polynomials are parsed as a vector of assigned coefficients [a_DEG, a_DEG-1, ..., a_1, a_0] where a_0 is the constant term
/// Q is the modulus of the Ring. All the coefficients will be in the range [0, Q-1]
pub fn div_euclid<const Q: u64>(dividend: &Vec<u64>, divisor: &Vec<u64>) -> (Vec<u64>, Vec<u64>) {
if divisor.is_empty() || divisor.iter().all(|&x| x == 0) {
panic!("Cannot divide by a zero polynomial!");
}
// transform the dividend and divisor into a vector of i64
let mut dividend = dividend.iter().map(|&x| x as i64).collect::<Vec<i64>>();
let mut divisor = divisor.iter().map(|&x| x as i64).collect::<Vec<i64>>();
// reverse the dividend and divisor
dividend.reverse();
divisor.reverse();
let divisor = divisor.iter().map(|&x| x as i64).collect::<Vec<i64>>();
let mut quotient = Vec::new();
let mut remainder = Vec::new();
@@ -50,14 +49,28 @@ pub fn div_euclid(dividend: &Vec<u64>, divisor: &Vec<u64>) -> (Vec<u64>, Vec<u64
}
// Convert remainder back to u64
let mut remainder = remainder.iter().map(|&x| x as u64).collect::<Vec<u64>>();
let remainder = remainder.iter().map(|&x| x as u64).collect::<Vec<u64>>();
// Convert quotient back to u64
let mut quotient = quotient.iter().map(|&x| x as u64).collect::<Vec<u64>>();
// reverse the quotient and remainder
quotient.reverse();
remainder.reverse();
let quotient = quotient.iter().map(|&x| x as u64).collect::<Vec<u64>>();
(quotient, remainder)
}
}
/// Convert a vector of AssignedValue to a vector of u64
/// Assumes that each element of AssignedValue can be represented in 8 bytes
pub fn vec_assigned_to_vec_u64<F: ScalarField>(vec: &Vec<AssignedValue<F>>) -> Vec<u64> {
let mut vec_u64 = Vec::new();
for i in 0..vec.len() {
let value_bytes_le = vec[i].value().to_bytes_le();
// slice value_to_bytes_le the first 8 bytes
let value_8_bytes_le = &value_bytes_le[..8];
let mut array_value_8_bytes_le = [0u8; 8];
array_value_8_bytes_le.copy_from_slice(value_8_bytes_le);
let num = u64::from_le_bytes(array_value_8_bytes_le);
vec_u64.push(num);
}
vec_u64
}