fn acc_column>()

in prover/src/constraints/evaluation_table.rs [281:341]

• 41 lines of code
• 3 McCabe index (conditional complexity)

fn acc_column<B: StarkField, E: FieldElement<BaseField = B>>(
    column: Vec<E>,
    divisor: &ConstraintDivisor<B>,
    domain_offset: B,
    result: &mut [E],
) {
    let numerator = divisor.numerator();
    assert_eq!(numerator.len(), 1, "complex divisors are not yet supported");
    assert!(
        divisor.exclude().len() <= 1,
        "multiple exclusion points are not yet supported"
    );

    // compute inverse evaluations of the divisor's numerator, which has the form (x^a - b)
    let domain_size = column.len();
    let z = get_inv_evaluation(divisor, domain_size, domain_offset);

    // divide column values by the divisor; for boundary constraints this computed simply as
    // multiplication of column value by the inverse of divisor numerator; for transition
    // constraints, it is computed similarly, but the result is also multiplied by the divisor's
    // denominator (exclusion point).
    if divisor.exclude().is_empty() {
        // the column represents merged evaluations of boundary constraints, and divisor has the
        // form of (x^a - b); thus to divide the column by the divisor, we compute: value * z,
        // where z = 1 / (x^a - 1) and has already been computed above.
        iter_mut!(result, 1024)
            .zip(column)
            .enumerate()
            .for_each(|(i, (acc_value, value))| {
                // determine which value of z corresponds to the current domain point
                let z = E::from(z[i % z.len()]);
                // compute value * z and add it to the result
                *acc_value += value * z;
            });
    } else {
        // the column represents merged evaluations of transition constraints, and divisor has the
        // form of (x^a - 1) / (x - b); thus, to divide the column by the divisor, we compute:
        // value * (x - b) * z, where z = 1 / (x^a - 1) and has already been computed above.

        // set up variables for computing x at every point in the domain
        let g = B::get_root_of_unity(domain_size.trailing_zeros());
        let b = divisor.exclude()[0];

        batch_iter_mut!(
            result,
            128, // min batch size
            |batch: &mut [E], batch_offset: usize| {
                let mut x = domain_offset * g.exp((batch_offset as u64).into());
                for (i, acc_value) in batch.iter_mut().enumerate() {
                    // compute value of (x - b) and compute next value of x
                    let e = x - b;
                    x *= g;
                    // determine which value of z corresponds to the current domain point
                    let z = z[i % z.len()];
                    // compute value * (x - b) * z and add it to the result
                    *acc_value += column[batch_offset + i] * E::from(z * e);
                }
            }
        );
    }
}