diff --git a/CHANGELOG.md b/CHANGELOG.md index e9bfa6806..720131102 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -29,6 +29,7 @@ - [\#691](https://github.com/arkworks-rs/algebra/pull/691) (`ark-poly`) Implement `Polynomial` for `SparseMultilinearExtension` and `DenseMultilinearExtension`. - [\#693](https://github.com/arkworks-rs/algebra/pull/693) (`ark-serialize`) Add `serialize_to_vec!` convenience macro. - [\#713](https://github.com/arkworks-rs/algebra/pull/713) (`ark-ff`) Add support for bitwise operations AND, OR, and XOR between `BigInteger`. +- [\#761](https://github.com/arkworks-rs/algebra/pull/761) (`ark-ff`) Implementation of modulus operator for `BigInteger`. ### Improvements diff --git a/ff/src/biginteger/mod.rs b/ff/src/biginteger/mod.rs index 46a236fcb..b253fc2fd 100644 --- a/ff/src/biginteger/mod.rs +++ b/ff/src/biginteger/mod.rs @@ -9,12 +9,13 @@ use ark_serialize::{ }; use ark_std::{ borrow::Borrow, + cmp::Ordering, convert::TryFrom, fmt::{Debug, Display, UpperHex}, io::{Read, Write}, ops::{ - BitAnd, BitAndAssign, BitOr, BitOrAssign, BitXor, BitXorAssign, Not, Shl, ShlAssign, Shr, - ShrAssign, + BitAnd, BitAndAssign, BitOr, BitOrAssign, BitXor, BitXorAssign, Not, Rem, RemAssign, Shl, + ShlAssign, Shr, ShrAssign, }, rand::{ distributions::{Distribution, Standard}, @@ -294,6 +295,7 @@ impl BigInt { impl BigInteger for BigInt { const NUM_LIMBS: usize = N; + const RADIX: u128 = 1 << 64; #[inline] fn add_with_carry(&mut self, other: &Self) -> bool { @@ -916,6 +918,328 @@ impl Not for BigInt { } } +/// Computes a division between to BigInt with remainder. +/// +/// This algorithm is taken from "The Art of Computer Programming" by Donald +/// Knuth, Third Edition, Section 4.3.2, Algorithm D, with some elements taken +/// from the repository https://github.com/rust-num/num-bigint. +fn div_with_rem>, const N: usize>( + a: BigInt, + d: B, +) -> (BigInt, BigInt) { + let d = d.borrow(); + + // Basic cases + if d.is_zero() { + panic!("trying to divide by zero"); + } + if d.is_one() { + return (a, BigInt::::zero()); + } + if a.is_zero() { + return (BigInt::::zero(), BigInt::::zero()); + } + + // Covers the case of N = 1, for which the quotient and remainder are + // computed as usual. + if N == 1 { + let q: u64 = a.0[0] / d.0[0]; + let r: u64 = a.0[0] % d.0[0]; + return (BigInt::::from(q), BigInt::::from(r)); + } + + // Covers immediate cases where a <= d. + match a.cmp(d) { + Ordering::Less => return (BigInt::::zero(), a), + Ordering::Equal => return (BigInt::::one(), BigInt::::zero()), + Ordering::Greater => {}, // This case is covered in what follows. + } + + // Here, we are in the case where a >= d, and it's not a trivial case. We + // will apply the TAOCP algorithm. + + // Considers the case where N is not 1, but the only non-zero limb of v is + // the least significant limb. We can apply unwrap here because we know that a + // is not zero, so there should be at least one limb with a non-zero value. + let index_non_zero_d = d.0.iter().rposition(|&e| e != 0).unwrap(); + if index_non_zero_d == 0 { + // Removes the most significant limbs with value zero from a. + let index_non_zero_a = a.0.iter().rposition(|&e| e != 0).unwrap(); + let u = a.0[..index_non_zero_a + 1].to_vec(); + + // Takes the only non-zero limb of v. + let v = d.0[index_non_zero_d]; + + return div_with_rem_one_digit(u, v); + } + + // Computes the normalization factor. + let shift = d.0[index_non_zero_d].leading_zeros() as u32; + + if shift == 0 { + div_with_rem_core(a.0.to_vec(), d.0.to_vec()) + } else { + // Here we construct the shifted version of a without chopping the overflow. + // Instead, we add another limb to store the overflow and putting it as + // the most significant limb. + + // Captures the most sifnificant limb of a before the shifting. Be careful: + // this limb can be zero. + let most_sig_a = *a.0.last().unwrap(); + + // Shift a + let shifted_a_bi = a << shift; // bi stands for BigInteger. + let mut a_shifted = shifted_a_bi.0.to_vec(); + + // Check if we need to add an overflow to a. If the most significant + // limb of a before the sifting is zero, we don't need to add an + // additional overflow. + if most_sig_a != 0 { + let overflow = most_sig_a >> (64 - shift); + a_shifted.push(overflow); + } else { + // If the most significant limb is zero we chop the most significant + // limbs that are zero after shifting a. + let index_non_zero_a_shifted = a_shifted.iter().rposition(|&e| e != 0).unwrap(); + a_shifted = a_shifted[0..index_non_zero_a_shifted + 1].to_vec(); + } + + // Apply the shift to d and chop the most significant limbs with value + // zero to d before applying the core algorithm. + let shifted_d_bi = *d << shift; + let d_shifted = shifted_d_bi.0[0..index_non_zero_d + 1].to_vec(); + + let (q, r) = div_with_rem_core(a_shifted, d_shifted); + return (q, r >> shift); + } +} + +/// Computes the division with remainder between u = [u_0, ..., u_(n + 1)] and v, where +/// v has just one non-zero limbs. +fn div_with_rem_one_digit(mut u: Vec, v: u64) -> (BigInt, BigInt) { + if v == 0 { + panic!("trying to divide by zero."); + } + + let mut rem = 0; + + for d in u.iter_mut().rev() { + let (q, r) = div_with_rem_wide(rem, *d, v); + *d = q; + rem = r; + } + + let mut q_result = [0; N]; + let mut r_result = [0; N]; + q_result.copy_from_slice(&u); + r_result[0] = rem; + + (BigInt::new(q_result), BigInt::new(r_result)) +} + +/// Divides a two-limb number by a one limb number returning the quotient and the remainder. +/// +/// More specifically, this method computes the divission between [lo, hi] / div. +#[inline] +fn div_with_rem_wide(hi: u64, lo: u64, div: u64) -> (u64, u64) { + let dividend = from_two_u64_to_u128(lo, hi); + let divisor = div as u128; + let q = (dividend / divisor) as u64; + let r = (dividend % divisor) as u64; + (q, r) +} + +/// Core of the algorithm for long division. +/// +/// Computes the long division between two BigInt. This function receives two +/// `Vec` given that it is more easy to manipulate the data this way. +/// In this algorithm, the notation is [less_sig, ..., most_sig] for +/// the radix-b representation. +fn div_with_rem_core(mut u: Vec, v: Vec) -> (BigInt, BigInt) { + // Length of the resulting quotient + let q_length = u.len() - v.len() + 1; + + // This vector will store the result of the quotient. + let mut q = vec![0; q_length]; + + let v0 = *v.last().unwrap(); // Most significant bit of v. + let v1 = v[v.len() - 2]; + + let mut u0 = 0; // Additional most significant limb as mentioned in Step D1. + + // We procceed as in the TAOCP algorithm, Steps D2-D7. + for j in (0..q_length).rev() { + let u1 = *u.last().unwrap(); // Most significant bit of u. + let u2 = u[u.len() - 2]; // We can do this given that N > 1. + + // Step D3. + let (mut q_hat, rem) = div_with_rem_wide(u0, u1, v0); + let mut r_hat = rem as u128; + loop { + if (q_hat as u128) == BigInt::::RADIX + || from_two_u64_to_u128(u2, r_hat as u64) < (q_hat as u128) * (v1 as u128) + { + q_hat -= 1; + r_hat += v0 as u128; // Promoted to u128 to avoid overflow. + } else { + break; + } + + // Test the previous condition again in the case that r_hat < u64::MAX, as + // specified in Step D3. + if r_hat >= BigInt::::RADIX { + break; + } + } + + // Step D4. Step D5 is not explicitly stated given that the result of the + // division has just one limb. + let mut borrow = sub_mul_with_borrow(&mut u[j..], &v, q_hat); + + // Step D6. + if borrow > u0 { + q_hat -= 1; + borrow -= add_different_len(&mut u[j..], &v); + } + + debug_assert!(borrow == u0); + + // Deletes the most significant bit of u. At the end, u will have length N. + u0 = u.pop().unwrap(); + + // Reconstruct BigInt + q[j] = q_hat; // It moves just one limb given that the quotient needs + // just one limb by definition. + } + + u.push(u0); + + // Collect the remainder and the quotient into arrays. + let mut remainder = [0; N]; + let mut quotient = [0; N]; + + // Append missing limbs to u and q + let miss_zeros_u = N - u.len(); + for _ in 0..miss_zeros_u { + u.push(0); + } + + let miss_zeros_q = N - q.len(); + for _ in 0..miss_zeros_q { + q.push(0); + } + + remainder.copy_from_slice(&u); + quotient.copy_from_slice(&q); + + (BigInt::new(quotient), BigInt::new(remainder)) +} + +/// Computes `u = u - q * v`, where `u = [u_0, ..., u_(n + 1)]`, `q` has just one limb and +/// `v = [v_0, ..., v_n]`, and returns the borrow. This source code is based in the one presented +/// in https://github.com/rust-num/num-bigint/blob/master/src/bigint/division.rs. +fn sub_mul_with_borrow(u: &mut [u64], v: &Vec, q: u64) -> u64 { + // We do this because the borrow is in the interval [-u64::MAX, 0]. Hence, to avoid overflow + // we compute the carry plus an offset. Then we retrieve the possitive value of the borrow + // by removing the offset. + let mut offset_carry = u64::MAX; + + for (x, y) in u.iter_mut().zip(v) { + // x - y * q + carry + let offset_sum = from_two_u64_to_u128(*x, u64::MAX) - (u64::MAX as u128) + + (offset_carry as u128) + - (*y as u128) * (q as u128); + + let (new_offset_carry, new_x) = get_u64_limbs_from_u128(offset_sum); + offset_carry = new_offset_carry; + + *x = new_x; + } + + u64::MAX - offset_carry +} + +/// Obtains the 64 most significant bits as an u68 value from a u128 value. +#[inline] +fn get_high_u64_from_u128(n: u128) -> u64 { + (n >> 64) as u64 +} + +/// Obtains the 64 least significant bits as an u68 value from a u128 value. +#[inline] +fn get_low_u64_from_u128(n: u128) -> u64 { + // This is equivalent to the 64 least significant bits set to 1 (0xFFFFFFFFFFFFFFFF) + const LOW_MASK: u128 = (1_u128 << 64) - 1; + (n & LOW_MASK) as u64 +} + +/// Obtains the two limbs of a u128 value as u64 values. +#[inline] +fn get_u64_limbs_from_u128(n: u128) -> (u64, u64) { + (get_high_u64_from_u128(n), get_low_u64_from_u128(n)) +} + +/// Computes the sum of `carry + a + b` storing the result in `output` and returning +/// the carry flag as `u8`. +#[inline] +fn add_with_carry(carry_in: u8, a: u64, b: u64, output: &mut u64) -> u8 { + let (out_a_b, carry_a_b) = a.overflowing_add(b); + let (out_ab_carry, carry_ab_carry) = out_a_b.overflowing_add(carry_in as u64); + let carry_sum = (carry_a_b as u8) + (carry_ab_carry as u8); + let (out_sum_carry, final_carry) = out_ab_carry.overflowing_add(carry_sum as u64); + *output = out_sum_carry; + final_carry as u8 +} + +/// Converts two u64 numbers `low` and `high` into an u128 representation of the number +/// `[low, high]`. +#[inline] +fn from_two_u64_to_u128(low: u64, high: u64) -> u128 { + (u128::from(high) << 64) | u128::from(low) +} + +/// Computes the addition two values in radix u64 and returns the carry. Those +/// two values have not necessarily the same number of limbs. This algorithm is +/// inspired from https://github.com/rust-num/num-bigint/blob/master/src/bigint/addition.rs. +fn add_different_len(u: &mut [u64], v: &Vec) -> u64 { + let mut carry = 0; + let (u_low, u_high) = u.split_at_mut(v.len()); + + for (u, v) in u_low.iter_mut().zip(v) { + carry = add_with_carry(carry, *u, *v, u); + } + + if carry != 0 { + for u in u_high { + carry = add_with_carry(carry, *u, 0, u); + if carry == 0 { + break; + } + } + } + + carry as u64 +} + +impl, const N: usize> Rem for BigInt { + type Output = Self; + + /// Computes the remainder between two `BigInt`. + fn rem(mut self, rhs: B) -> Self::Output { + self %= rhs; + self + } +} + +impl, const N: usize> RemAssign for BigInt { + /// Computes the remainder between `self` and `rhs`, assigning the result to + /// `self`. + fn rem_assign(&mut self, rhs: B) { + let (_, rem) = div_with_rem(*self, rhs); + (0..N).for_each(|i| self.0[i] = rem.0[i]); + } +} + /// Compute the signed modulo operation on a u64 representation, returning the result. /// If n % modulus > modulus / 2, return modulus - n /// # Example @@ -989,10 +1313,18 @@ pub trait BigInteger: + ShrAssign + Shl + ShlAssign + + Rem + + for<'a> Rem<&'a Self, Output = Self> + + RemAssign + + for<'a> RemAssign<&'a Self> { /// Number of 64-bit limbs representing `Self`. const NUM_LIMBS: usize; + /// Radix of the representation. This is the number of different "digits" that + /// can be represented by one limb. + const RADIX: u128; + /// Add another [`BigInteger`] to `self`. This method stores the result in `self`, /// and returns a carry bit. /// @@ -1226,6 +1558,19 @@ pub trait BigInteger: /// ``` fn is_zero(&self) -> bool; + /// Returns true iff this number is one. + /// # Example + /// + /// ``` + /// use ark_ff::{biginteger::BigInteger64 as B, BigInteger as _}; + /// + /// let mut one = B::from(1u64); + /// assert!(one.is_one()); + /// ``` + fn is_one(&self) -> bool { + *self == Self::from(1u64) + } + /// Compute the minimum number of bits needed to encode this number. /// # Example /// ``` diff --git a/ff/src/biginteger/tests.rs b/ff/src/biginteger/tests.rs index dad4730c1..22f9d5ce1 100644 --- a/ff/src/biginteger/tests.rs +++ b/ff/src/biginteger/tests.rs @@ -9,6 +9,10 @@ fn biginteger_arithmetic_test(a: B, b: B, zero: B, max: B) { // zero.is_zero() == true assert_eq!(zero.is_zero(), true); + // one.is_one() == true + let one = B::from(1u64); + assert!(one.is_one()); + // a == a assert_eq!(a, a); @@ -139,6 +143,38 @@ fn biginteger_shl() { assert_eq!(b.get_bit(2), false); } +fn biginteger_modulus() { + let mut rng = ark_std::test_rng(); + + // Basic tests + let a = B::rand(&mut rng); + assert_eq!(B::from(0u64) % a, B::from(0u64)); + assert_eq!(a % B::from(1u64), B::from(0u64)); + + // Check that the modulus computation agrees with parity test. + if a.is_even() { + assert_eq!(a % B::from(2u64), B::from(0u64)); + } else { + assert_eq!(a % B::from(2u64), B::from(1u64)); + } + + // Check modulus identity. + let a = B::rand(&mut rng); + let b = B::rand(&mut rng); + let a_mod_b = a % b; + let a_mod_b_mod_b = a_mod_b % b; + assert_eq!(a_mod_b, a_mod_b_mod_b); + + // Check trivial cases from comparisons. + let a = B::rand(&mut rng); + let b = B::rand(&mut rng); + match a.cmp(&b) { + core::cmp::Ordering::Less => assert_eq!(a % b, a), + core::cmp::Ordering::Equal => assert_eq!(a % b, B::from(0u64)), + core::cmp::Ordering::Greater => {}, + } +} + // Test for BigInt's bitwise operations fn biginteger_bitwise_ops_test() { let mut rng = ark_std::test_rng(); @@ -230,6 +266,160 @@ fn test_biginteger(max: B, zero: B) { biginteger_bitwise_ops_test::(); biginteger_shr::(); biginteger_shl::(); + biginteger_modulus::(); +} + +#[test] +fn biginteger_modulus_static_tests() { + use crate::biginteger::BigInt; + + let a64 = BigInt::new([3601073082870250418]); + let b64 = BigInt::new([2193000201251463648]); + let mod64 = BigInt::new([1408072881618786770]); + assert_eq!(a64 % b64, mod64); + + let a128 = BigInt::new([4174451077990245596, 3669629040712867839]); + let b128 = BigInt::new([8908506066818206150, 17331115722498298823]); + let mod128 = BigInt::new([4174451077990245596, 3669629040712867839]); + assert_eq!(a128 % b128, mod128); + + let a256 = BigInt::new([ + 14621467787333493418, + 12108485322752294880, + 11620129110562858835, + 6771525239094520794, + ]); + let b256 = BigInt::new([ + 9457583366609244046, + 3292446216252191244, + 15683190688725518430, + 12376509759692993049, + ]); + let mod256 = BigInt::new([ + 14621467787333493418, + 12108485322752294880, + 11620129110562858835, + 6771525239094520794, + ]); + assert_eq!(a256 % b256, mod256); + + let a512 = BigInt::new([ + 12510516664288822641, + 14513930181904105826, + 13449486877885219817, + 8005255619594419542, + 6228641990551128361, + 8770155517167400132, + 15617099880276166326, + 17143852322472171649, + ]); + let b512 = BigInt::new([ + 14002016314341253844, + 18398549445620213883, + 5645345381311304630, + 16592555506467753304, + 16984857277391876626, + 16427282640994978259, + 11465536577168429635, + 698543981412520453, + ]); + let mod512 = BigInt::new([ + 8503518446870659473, + 15670601256048211400, + 7088406242380769985, + 15612293085978475791, + 4420436954756224867, + 1896997681188505830, + 17145383133877129305, + 378796768571680762, + ]); + assert_eq!(a512 % b512, mod512); + + let a1024 = BigInt::new([ + 813765468505266982, + 9461391675911782696, + 9082012781690260231, + 8393500664090969669, + 4133856364806303958, + 11149080897905582441, + 14731346490539363185, + 989680701211629440, + 8526880266768469178, + 2316258228357047461, + 5239140938378412139, + 13856807587119809799, + 17095289661235426508, + 15245639872422554895, + 15979549614507724819, + 1709763898347752414, + ]); + let b1024 = BigInt::new([ + 14373198534555697120, + 1857962328286269016, + 15002453475225291270, + 12505120393686880765, + 17843441582302809739, + 8589716293423618750, + 684922508168804796, + 3583840145186657945, + 14969815686562072035, + 13969943195137168376, + 11597877104598378561, + 8138617899241207781, + 11747787798498493404, + 10447752210747488036, + 2262843315480796728, + 6590054619087534118, + ]); + let mod1024 = BigInt::new([ + 813765468505266982, + 9461391675911782696, + 9082012781690260231, + 8393500664090969669, + 4133856364806303958, + 11149080897905582441, + 14731346490539363185, + 989680701211629440, + 8526880266768469178, + 2316258228357047461, + 5239140938378412139, + 13856807587119809799, + 17095289661235426508, + 15245639872422554895, + 15979549614507724819, + 1709763898347752414, + ]); + assert_eq!(a1024 % b1024, mod1024); +} + +#[test] +fn test_from_two_u64_to_128() { + use crate::biginteger::from_two_u64_to_u128; + + let mut rng = ark_std::test_rng(); + let a = u64::rand(&mut rng); + let b = u64::rand(&mut rng); + + let r = from_two_u64_to_u128(a, b); + + let ref_val = (b as u128) * ((u64::MAX as u128) + 1) + (a as u128); + assert_eq!(ref_val, r); +} + +#[test] +fn test_conversion_u64_u128() { + use crate::biginteger::{from_two_u64_to_u128, get_u64_limbs_from_u128}; + + let mut rng = ark_std::test_rng(); + let a = u64::rand(&mut rng); + let b = u64::rand(&mut rng); + + let r = from_two_u64_to_u128(a, b); + + let (b_hat, a_hat) = get_u64_limbs_from_u128(r); + + assert_eq!(a, a_hat); + assert_eq!(b, b_hat); } #[test]