Skip to content
This repository has been archived by the owner on Dec 9, 2021. It is now read-only.

Commit

Permalink
Issue 68: Refactor existing modular arithmetic code (#80)
Browse files Browse the repository at this point in the history 
* Add test case generation script.

* Change test case names.

* Move gentests to subdirectory; ignore generated files.

* Document gentests.py file output format.

* Add BLS12-377 base field to 'gentests.py'.

* First more-or-less full version of test suite for Tweedledum.

* Ensure input is normalised.

* Support unary operators in test suite.

* Use zero instead of modulus for the 'None' value.

* Finish unaryop/binaryop refactor; make division work.

* Remove dependency on Tweedledum curve.

* Move test utils out of TweedledumBase (on their way to their own file eventually).

* Generate tests with a macro; instantiate in each base field.

* Use Field::BYTES when reading in test cases.

* Generate test inputs in a function.

* Completed test suite.

* Remove unused files.

* Document some functions.

* Addressed PR comments.

* Used const generics to refactor bigint arithmetic.

* Fix bad merge.

* Remove `mul` use clause and fix `cmp` visibility.
  • Loading branch information
@unzvfu
unzvfu authored Sep 27, 2020
1 parent 1bd0f5c commit a676a71
Show file tree
Hide file tree
Showing 8 changed files with 139 additions and 367 deletions.
12 changes: 4 additions & 8 deletions benches/bigint_arithmetic.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,26 +2,22 @@ use criterion::{black_box, Criterion};
use criterion::criterion_group;
use criterion::criterion_main;

use plonky::{cmp_6_6, mul_6_6};
use plonky::cmp;

fn criterion_benchmark(c: &mut Criterion) {
let x = [11111111, 22222222, 33333333, 44444444, 55555555, 66666666];
let y = [44444444, 55555555, 66666666, 77777777, 88888888, 99999999];

c.bench_function("[u64; 6] widening multiplication", move |b| b.iter(|| {
mul_6_6(black_box(x), black_box(y))
}));

c.bench_function("[u64; 6] comparison (lhs == rhs)", move |b| b.iter(|| {
cmp_6_6(black_box(x), black_box(x))
cmp(black_box(x), black_box(x))
}));

c.bench_function("[u64; 6] comparison (lhs < rhs)", move |b| b.iter(|| {
cmp_6_6(black_box(x), black_box(y))
cmp(black_box(x), black_box(y))
}));

c.bench_function("[u64; 6] comparison (lhs > rhs)", move |b| b.iter(|| {
cmp_6_6(black_box(y), black_box(x))
cmp(black_box(y), black_box(x))
}));
}

Expand Down
253 changes: 37 additions & 216 deletions src/bigint/bigint_arithmetic.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,22 +8,8 @@ use unroll::unroll_for_loops;
/// This module provides functions for big integer arithmetic using little-endian encoded u64
/// arrays.
#[unroll_for_loops]
pub(crate) fn cmp_4_4(a: [u64; 4], b: [u64; 4]) -> Ordering {
for i in (0..4).rev() {
if a[i] < b[i] {
return Less;
}
if a[i] > b[i] {
return Greater;
}
}
Equal
}

/// Public only for use with Criterion benchmarks.
#[unroll_for_loops]
pub fn cmp_6_6(a: [u64; 6], b: [u64; 6]) -> Ordering {
for i in (0..6).rev() {
pub fn cmp<const N: usize>(a: [u64; N], b: [u64; N]) -> Ordering {
for i in (0..N).rev() {
if a[i] < b[i] {
return Less;
}
Expand All @@ -36,27 +22,10 @@ pub fn cmp_6_6(a: [u64; 6], b: [u64; 6]) -> Ordering {

/// Computes `a + b`. Assumes that there is no overflow; this is verified only in debug builds.
#[unroll_for_loops]
pub(crate) fn add_4_4_no_overflow(a: [u64; 4], b: [u64; 4]) -> [u64; 4] {
pub(crate) fn add_no_overflow<const N: usize>(a: [u64; N], b: [u64; N]) -> [u64; N] {
let mut carry = false;
let mut sum = [0; 4];
for i in 0..4 {
// First add a[i] + b[i], then add in the carry.
let result1 = a[i].overflowing_add(b[i]);
let result2 = result1.0.overflowing_add(carry as u64);
sum[i] = result2.0;
// If either sum overflowed, set the carry bit.
carry = result1.1 | result2.1;
}
debug_assert_eq!(carry, false, "Addition overflowed");
sum
}

/// Computes `a + b`. Assumes that there is no overflow; this is verified only in debug builds.
#[unroll_for_loops]
pub(crate) fn add_6_6_no_overflow(a: [u64; 6], b: [u64; 6]) -> [u64; 6] {
let mut carry = false;
let mut sum = [0; 6];
for i in 0..6 {
let mut sum = [0; N];
for i in 0..N {
// First add a[i] + b[i], then add in the carry.
let result1 = a[i].overflowing_add(b[i]);
let result2 = result1.0.overflowing_add(carry as u64);
Expand All @@ -70,10 +39,10 @@ pub(crate) fn add_6_6_no_overflow(a: [u64; 6], b: [u64; 6]) -> [u64; 6] {

/// Computes `a - b`. Assumes `a >= b`, otherwise the behavior is undefined.
#[unroll_for_loops]
pub(crate) fn sub_4_4(a: [u64; 4], b: [u64; 4]) -> [u64; 4] {
pub(crate) fn sub<const N: usize>(a: [u64; N], b: [u64; N]) -> [u64; N] {
let mut borrow = false;
let mut difference = [0; 4];
for i in 0..4 {
let mut difference = [0; N];
for i in 0..N {
let result1 = a[i].overflowing_sub(b[i]);
let result2 = result1.0.overflowing_sub(borrow as u64);
difference[i] = result2.0;
Expand All @@ -85,121 +54,25 @@ pub(crate) fn sub_4_4(a: [u64; 4], b: [u64; 4]) -> [u64; 4] {
difference
}

/// Computes `a - b`. Assumes `a >= b`, otherwise the behavior is undefined.
#[unroll_for_loops]
pub(crate) fn sub_6_6(a: [u64; 6], b: [u64; 6]) -> [u64; 6] {
let mut borrow = false;
let mut difference = [0; 6];

for i in 0..6 {
let result1 = a[i].overflowing_sub(b[i]);
let result2 = result1.0.overflowing_sub(borrow as u64);
difference[i] = result2.0;
// If either difference underflowed, set the borrow bit.
borrow = result1.1 | result2.1;
}

debug_assert!(!borrow, "a < b: {:?} < {:?}", a, b);
difference
}

#[allow(dead_code)]
#[unroll_for_loops]
pub(crate) fn mul_4_4(a: [u64; 4], b: [u64; 4]) -> [u64; 8] {
// Grade school multiplication. To avoid carrying at each of O(n^2) steps, we first add each
// intermediate product to a 128-bit accumulator, then propagate carries at the end.
let mut acc128 = [0u128; 4 + 4];

for i in 0..4 {
for j in 0..4 {
let a_i_b_j = a[i] as u128 * b[j] as u128;
// Add the less significant chunk to the less significant accumulator.
acc128[i + j] += a_i_b_j as u64 as u128;
// Add the more significant chunk to the more significant accumulator.
acc128[i + j + 1] += a_i_b_j >> 64;
}
}

let mut acc = [0u64; 8];
acc[0] = acc128[0] as u64;
let mut carry = false;
for i in 1..8 {
let last_chunk_big = (acc128[i - 1] >> 64) as u64;
let curr_chunk_small = acc128[i] as u64;
// Note that last_chunk_big won't get anywhere near 2^64, since it's essentially a carry
// from some additions in the previous phase, so we can add the carry bit to it without
// fear of overflow.
let result = curr_chunk_small.overflowing_add(last_chunk_big + carry as u64);
acc[i] += result.0;
carry = result.1;
}
debug_assert!(!carry);
acc
}

/// Public only for use with Criterion benchmarks.
#[unroll_for_loops]
pub fn mul_6_6(a: [u64; 6], b: [u64; 6]) -> [u64; 6 + 6] {
// Grade school multiplication. To avoid carrying at each of O(n^2) steps, we first add each
// intermediate product to a 128-bit accumulator, then propagate carries at the end.
let mut acc128 = [0u128; 6 + 6];

for i in 0..6 {
for j in 0..6 {
let a_i_b_j = a[i] as u128 * b[j] as u128;
// Add the less significant chunk to the less significant accumulator.
acc128[i + j] += a_i_b_j as u64 as u128;
// Add the more significant chunk to the more significant accumulator.
acc128[i + j + 1] += a_i_b_j >> 64;
}
}

let mut acc = [0u64; 6 + 6];
acc[0] = acc128[0] as u64;
let mut carry = false;
for i in 1..(6 + 6) {
let last_chunk_big = (acc128[i - 1] >> 64) as u64;
let curr_chunk_small = acc128[i] as u64;
// Note that last_chunk_big won't get anywhere near 2^64, since it's essentially a carry
// from some additions in the previous phase, so we can add the carry bit to it without
// fear of overflow.
let result = curr_chunk_small.overflowing_add(last_chunk_big + carry as u64);
acc[i] += result.0;
carry = result.1;
}
debug_assert!(!carry);
acc
}

#[inline(always)]
pub(crate) fn is_even_4(x: [u64; 4]) -> bool {
x[0] & 1 == 0
}

#[inline(always)]
pub(crate) fn is_even_6(x: [u64; 6]) -> bool {
pub(crate) fn is_even<const N: usize>(x: [u64; N]) -> bool {
x[0] & 1 == 0
}

#[inline(always)]
pub(crate) fn is_odd_4(x: [u64; 4]) -> bool {
x[0] & 1 != 0
}

#[inline(always)]
pub(crate) fn is_odd_6(x: [u64; 6]) -> bool {
pub(crate) fn is_odd<const N: usize>(x: [u64; N]) -> bool {
x[0] & 1 != 0
}

/// Shift in the direction of increasing significance by one bit. Equivalent to integer
/// multiplication by two. Assumes that the input's most significant bit is clear.
#[unroll_for_loops]
pub(crate) fn mul2_6(x: [u64; 6]) -> [u64; 6] {
debug_assert_eq!(x[5] >> 63, 0, "Most significant bit should be clear");
pub(crate) fn mul2<const N: usize>(x: [u64; N]) -> [u64; N] {
debug_assert_eq!(x[N-1] >> 63, 0, "Most significant bit should be clear");

let mut result = [0; 6];
let mut result = [0; N];
result[0] = x[0] << 1;
for i in 1..6 {
for i in 1..N {
result[i] = x[i] << 1 | x[i - 1] >> 63;
}
result
Expand All @@ -208,114 +81,62 @@ pub(crate) fn mul2_6(x: [u64; 6]) -> [u64; 6] {
/// Shift in the direction of increasing significance by one bit. Equivalent to integer
/// division by two.
#[unroll_for_loops]
pub(crate) fn div2_4(x: [u64; 4]) -> [u64; 4] {
let mut result = [0; 4];
for i in 0..3 {
pub(crate) fn div2<const N: usize>(x: [u64; N]) -> [u64; N] {
let mut result = [0; N];
for i in 0..N-1 {
result[i] = x[i] >> 1 | x[i + 1] << 63;
}
result[3] = x[3] >> 1;
result[N-1] = x[N-1] >> 1;
result
}

/// Shift in the direction of increasing significance by one bit. Equivalent to integer
/// division by two.
#[unroll_for_loops]
pub(crate) fn div2_6(x: [u64; 6]) -> [u64; 6] {
let mut result = [0; 6];
for i in 0..5 {
result[i] = x[i] >> 1 | x[i + 1] << 63;
}
result[5] = x[5] >> 1;
result
}

pub(crate) fn rand_range_6(limit_exclusive: [u64; 6]) -> [u64; 6] {
rand_range_6_from_rng(limit_exclusive, &mut OsRng)
pub(crate) fn rand_range<const N: usize>(limit_exclusive: [u64; N]) -> [u64; N] {
rand_range_from_rng(limit_exclusive, &mut OsRng)
}

// Same as `rand_range_6` but specifying a RNG (useful when dealing with seeded RNGs).
pub(crate) fn rand_range_6_from_rng<R: Rng>(limit_exclusive: [u64; 6], rng: &mut R) -> [u64; 6] {
// Same as `rand_range` but specifying a RNG (useful when dealing with seeded RNGs).
pub(crate) fn rand_range_from_rng<R: Rng, const N: usize>(limit_exclusive: [u64; N], rng: &mut R) -> [u64; N] {
// Our approach is to repeatedly generate random u64 arrays until one of them happens to be
// within the limit. This could take a lot of attempts if the limit has many leading zero bits,
// though. It is more efficient to generate n-bit random numbers, where n is the number of bits
// in the limit.
let bits_to_strip = limit_exclusive[5].leading_zeros();
let bits_to_strip = limit_exclusive[N-1].leading_zeros();

let mut limbs = [0; 6];
let mut limbs = [0; N];

loop {
for limb_i in &mut limbs {
*limb_i = rng.next_u64();
}
limbs[5] >>= bits_to_strip;
limbs[N-1] >>= bits_to_strip;

if cmp_6_6(limbs, limit_exclusive) == Less {
if cmp(limbs, limit_exclusive) == Less {
return limbs;
}
}
}

pub(crate) fn rand_range_4(limit_exclusive: [u64; 4]) -> [u64; 4] {
rand_range_4_from_rng(limit_exclusive, &mut OsRng)
}

// Same as `rand_range_4` but specifying a RNG (useful when dealing with seeded RNGs).
pub(crate) fn rand_range_4_from_rng<R: Rng>(limit_exclusive: [u64; 4], rng: &mut R) -> [u64; 4] {
// Our approach is to repeatedly generate random u64 arrays until one of them happens to be
// within the limit. This could take a lot of attempts if the limit has many leading zero bits,
// though. It is more efficient to generate n-bit random numbers, where n is the number of bits
// in the limit.
let bits_to_strip = limit_exclusive[3].leading_zeros();

let mut limbs = [0; 4];

loop {
for limb_i in &mut limbs {
*limb_i = rng.next_u64();
}
limbs[3] >>= bits_to_strip;

if cmp_4_4(limbs, limit_exclusive) == Less {
return limbs;
#[macro_export]
macro_rules! one_array {
[$type:ty; $N:expr] => {
{
let mut tmp: [$type; $N] = [0 as $type; $N];
tmp[0] = 1;
tmp
}
}
};
}

#[cfg(test)]
mod tests {
use crate::conversions::u64_slice_to_biguint;
use crate::{div2_6, mul_6_6};

#[test]
fn test_mul_6_6() {
let a = [
11111111u64,
22222222,
33333333,
44444444,
55555555,
66666666,
];
let b = [
77777777u64,
88888888,
99999999,
11111111,
22222222,
33333333,
];
assert_eq!(
u64_slice_to_biguint(&mul_6_6(a, b)),
u64_slice_to_biguint(&a) * u64_slice_to_biguint(&b)
);
}
use crate::div2;

#[test]
fn test_div2() {
assert_eq!(div2_6([40, 0, 0, 0, 0, 0]), [20, 0, 0, 0, 0, 0]);
assert_eq!(div2([40, 0, 0, 0, 0, 0]), [20, 0, 0, 0, 0, 0]);

assert_eq!(
div2_6([
div2([
15668009436471190370,
3102040391300197453,
4166322749169705801,
Expand Down
Loading

0 comments on commit a676a71

Please sign in to comment.