diff --git a/algebra-core/Cargo.toml b/algebra-core/Cargo.toml index 953593715..d17b113e6 100644 --- a/algebra-core/Cargo.toml +++ b/algebra-core/Cargo.toml @@ -30,7 +30,6 @@ rand = { version = "0.7", default-features = false } rayon = { version = "1", optional = true } unroll = { version = "=0.1.4" } itertools = { version = "0.9.0", default-features = false } -voracious_radix_sort = { version = "0.1.0", optional = true } either = { version = "1.6.0", default-features = false } thread-id = { version = "3.3.0", optional = true } backtrace = { version = "0.3", optional = true } @@ -46,7 +45,7 @@ rand_xorshift = "0.2" [features] default = [ "std", "rand/default" ] -std = [ "voracious_radix_sort" ] +std = [] parallel = [ "std", "rayon", "rand/default" ] derive = [ "algebra-core-derive" ] prefetch = [ "std" ] diff --git a/algebra-core/src/curves/batch_arith.rs b/algebra-core/src/curves/batch_arith.rs index 4bfa80ab9..07c4cf630 100644 --- a/algebra-core/src/curves/batch_arith.rs +++ b/algebra-core/src/curves/batch_arith.rs @@ -7,6 +7,9 @@ use num_traits::Zero; /// inversion close to zero while not straining the CPU cache by generating and /// fetching from large w-NAF tables and slices [G] pub const BATCH_SIZE: usize = 4096; + +/// We code this in the second operand for the `batch_add_in_place_read_only` +/// method utilised in the `batch_scalar_mul_in_place` method. /// 0 == Identity; 1 == Neg; 2 == GLV; 3 == GLV + Neg pub const ENDO_CODING_BITS: usize = 2; @@ -29,8 +32,9 @@ where // Refer to e.g. Improved Techniques for Fast Exponentiation, Section 4 // Bodo M¨oller 2002. https://www.bmoeller.de/pdf/fastexp-icisc2002.pdf - /// Computes [[p, 3 * p, ..., (2^w - 1) * p], ..., [q, 3* q, ..., ]] - /// We need to manipulate the offsets when using the table + /// Computes [[p_1, 3 * p_1, ..., (2^w - 1) * p_1], ..., [p_n, 3*p_n, ..., + /// (2^w - 1) p_n]] We need to manipulate the offsets when using the + /// table fn batch_wnaf_tables(bases: &[Self], w: usize) -> Vec { let half_size = 1 << (w - 1); let batch_size = bases.len(); @@ -156,9 +160,9 @@ where /// or simply writes them to new_elems, using scratch space to store /// intermediate values. Scratch space is always cleared after use. - /// No-ops, or copies of the elem in the slice `lookup` in the position of the index - /// of the first operand to the new_elems vector, are encoded as !0u32 in the index - /// for the second operand + /// No-ops, or copies of the elem in the slice `lookup` in the position of + /// the index of the first operand to the new_elems vector, are encoded + /// as !0u32 in the index for the second operand fn batch_add_write( lookup: &[Self], index: &[(u32, u32)], @@ -169,9 +173,9 @@ where /// Similar to batch_add_write, only that the lookup for the first operand /// is performed in new_elems rather than lookup - /// No-ops, or copies of the elem in the slice `lookup` in the position of the index - /// of the first operand to the new_elems vector, are encoded as !0u32 in the index - /// for the second operand + /// No-ops, or copies of the elem in the slice `lookup` in the position of + /// the index of the first operand to the new_elems vector, are encoded + /// as !0u32 in the index for the second operand fn batch_add_write_read_self( lookup: &[Self], index: &[(u32, u32)], diff --git a/algebra-core/src/curves/bucketed_add.rs b/algebra-core/src/curves/bucketed_add.rs index 171a294cc..43b82c66d 100644 --- a/algebra-core/src/curves/bucketed_add.rs +++ b/algebra-core/src/curves/bucketed_add.rs @@ -3,42 +3,32 @@ use crate::{ AffineCurve, Vec, }; -#[cfg(feature = "std")] -use {core::cmp::Ordering, voracious_radix_sort::*}; - -#[cfg(not(feature = "std"))] -use crate::log2; - #[derive(Copy, Clone, Debug)] pub struct BucketPosition { pub bucket: u32, pub position: u32, } -#[cfg(feature = "std")] -impl PartialOrd for BucketPosition { - fn partial_cmp(&self, other: &Self) -> Option { - self.bucket.partial_cmp(&other.bucket) - } -} - -#[cfg(feature = "std")] -impl Radixable for BucketPosition { - type Key = u32; - #[inline] - fn key(&self) -> Self::Key { - self.bucket - } -} - -impl PartialEq for BucketPosition { - fn eq(&self, other: &Self) -> bool { - self.bucket == other.bucket - } -} - +/// The objective of this function is to identify an addition tree of +/// independent elliptic curve group additions for each bucket, and to batch the +/// independent additions using the batch affine inversion method. + +/// The strategy taken is to sort a list of bucket assignments of all the +/// elements (which we can for most intents and purposes, think of as being +/// uniformly random) by bucket, so that indices corresponding to elements that +/// must be added together are physically collocated in memory. Then, in the +/// first round, we proceed to perform independent additions producing +/// intermediate results at the greatest depth for each addition tree (each +/// corresponding to a bucket), and write the result to a new vector. We do so +/// to improve cache locality for future rounds, and take advantage of the +/// CPU-intensive nature of elliptic curve operations along with prfetching to +/// hide the latency of reading from essentially random locations in memory. + +/// Subsequently, we perform the additions in place, and the second operands +/// become junk data. Finally, when we only have the buckets left (no more +/// additions left to perform), we copy the result into a destination `res` +/// slice. #[inline] -#[cfg(feature = "std")] pub fn batch_bucketed_add( buckets: usize, elems: &[C], @@ -48,8 +38,11 @@ pub fn batch_bucketed_add( assert!(elems.len() > 0); let _now = timer!(); - dlsd_radixsort(bucket_positions, 8); - timer_println!(_now, "radixsort"); + // We sort the bucket positions so that indices of elements assigned + // to the same bucket are continguous. This way, we can easily identify + // how to construct the addition tree for that bucket. + bucket_positions.sort_unstable_by_key(|x| x.bucket); + timer_println!(_now, "sort"); let mut len = bucket_positions.len(); let mut all_ones = true; @@ -68,29 +61,45 @@ pub fn batch_bucketed_add( // Subsequently, we perform all the operations in place while glob < len { let current_bucket = bucket_positions[glob].bucket; + // We are iterating over elements using a global `glob` counter, and counting + // how many in a row are being assigned to the same bucket, using the `loc` + // counter. while glob + 1 < len && bucket_positions[glob + 1].bucket == current_bucket { glob += 1; loc += 1; } + // If the current bucket exceeds buckets, it encodes a noop if current_bucket >= buckets as u32 { loc = 1; } else if loc > 1 { // all ones is false if next len is not 1 + + // in other words, we have not reached the terminating + // condition that after the current round of addition + // there is only one element left in each addition tree + + // This would be the case, if each addition tree had at + // most 2 elements in the current round. if loc > 2 { all_ones = false; } let is_odd = loc % 2 == 1; let half = loc / 2; + // We encode instructions to add adjacent elements for i in 0..half { instr.push(( bucket_positions[glob - (loc - 1) + 2 * i].position, bucket_positions[glob - (loc - 1) + 2 * i + 1].position, )); + // Compactification of buckets bucket_positions[new_len + i] = BucketPosition { bucket: current_bucket, position: (new_len + i) as u32, }; } + // If there are an odd number of elements, the lone element + // without a partner will be copied over to the `new_elems` + // vector, a noop which is encoded as !0u32 if is_odd { instr.push((bucket_positions[glob].position, !0u32)); bucket_positions[new_len + half] = BucketPosition { @@ -99,6 +108,13 @@ pub fn batch_bucketed_add( }; } // Reset the local_counter and update state + + // We compactify the `bucket_positions` data by shifing left + // `new_len` is the len of the current compactified vector. + + // We also update the `batch` counter to decide when it is + // optimal to invoke the batch inversion, i.e. when we have + // accumulated enough independent additions. new_len += half + (loc % 2); batch += half; loc = 1; @@ -131,6 +147,9 @@ pub fn batch_bucketed_add( len = new_len; new_len = 0; + // We repeat the above procedure, except, since we are performing the addition + // trees in place, we do not need to encode noops to force a copy to a new + // vector. while !all_ones { all_ones = true; while glob < len { @@ -197,79 +216,3 @@ pub fn batch_bucketed_add( timer_println!(_now, "reassign"); res } - -#[cfg(not(feature = "std"))] -pub fn batch_bucketed_add( - buckets: usize, - elems: &[C], - bucket_assign: &[BucketPosition], -) -> Vec { - let mut elems = elems.to_vec(); - let num_split = 2i32.pow(log2(buckets) / 2 + 2) as usize; - let split_size = (buckets - 1) / num_split + 1; - let ratio = elems.len() / buckets * 2; - // Get the inverted index for the positions assigning to each bucket - let mut bucket_split = vec![vec![]; num_split]; - let mut index = vec![Vec::with_capacity(ratio); buckets]; - - for bucket_pos in bucket_assign.iter() { - let (bucket, position) = (bucket_pos.bucket as usize, bucket_pos.position as usize); - // Check the bucket assignment is valid - if bucket < buckets { - // index[bucket].push(position); - bucket_split[bucket / split_size].push((bucket, position)); - } - } - - for split in bucket_split { - for (bucket, position) in split { - index[bucket].push(position as u32); - } - } - - // Instructions for indexes for the in place addition tree - let mut instr: Vec> = vec![]; - // Find the maximum depth of the addition tree - let max_depth = index.iter() - // log_2 - .map(|x| log2(x.len())) - .max().unwrap(); - - // Generate in-place addition instructions that implement the addition tree - // for each bucket from the leaves to the root - for i in 0..max_depth { - let mut instr_row = Vec::<(u32, u32)>::with_capacity(buckets); - for to_add in index.iter_mut() { - if to_add.len() > 1 << (max_depth - i - 1) { - let mut new_to_add = vec![]; - for j in 0..(to_add.len() / 2) { - new_to_add.push(to_add[2 * j]); - instr_row.push((to_add[2 * j], to_add[2 * j + 1])); - } - if to_add.len() % 2 == 1 { - new_to_add.push(*to_add.last().unwrap()); - } - *to_add = new_to_add; - } - } - instr.push(instr_row); - } - - for instr_row in instr.iter() { - for instr in C::get_chunked_instr::<(u32, u32)>(&instr_row[..], BATCH_SIZE).iter() { - elems[..].batch_add_in_place_same_slice(&instr[..]); - } - } - - let zero = C::zero(); - let mut res = vec![zero; buckets]; - - for (i, to_add) in index.iter().enumerate() { - if to_add.len() == 1 { - res[i] = elems[to_add[0] as usize]; - } else if to_add.len() > 1 { - debug_assert!(false, "Did not successfully reduce to_add"); - } - } - res -} diff --git a/algebra-core/src/curves/glv.rs b/algebra-core/src/curves/glv.rs index c3993f357..eb4af4a35 100644 --- a/algebra-core/src/curves/glv.rs +++ b/algebra-core/src/curves/glv.rs @@ -3,7 +3,7 @@ use core::ops::Neg; /// The GLV parameters here require the following conditions to be satisfied: /// 1. MODULUS_BITS < NUM_LIMBS * 64 - 1. So 2 * n < 1 << (64 * NUM_LIMBS) -/// We also assume that |b1| * |b2| < 2 * n +/// We also assume that (|b1| + 2) * (|b2| + 2) < 2 * n /// We also know that either B1 is neg or B2 is. pub trait GLVParameters: Send + Sync + 'static + ModelParameters { type WideBigInt: BigInteger; @@ -53,9 +53,12 @@ pub trait GLVParameters: Send + Sync + 'static + ModelParameters { let c2 = &c2_wide.as_ref()[..limbs]; // We first assume that the final 2 bits of the representation for the modulus - // is not set, so that 2 * n < R = 1 << (64 * NUM_LIMBS). Then, since we - // know that |b_i| < \sqrt{2n}, wlog k|b1|/n * |b2| < 2 * k < 2 * n < - // R. + // is not set, so that 2 * n < R = 1 << (64 * NUM_LIMBS). + + // wlog c1 = round(k * round(|b_1|R / n) / R) < ceil(k * ceil(|b_1|* R / n) / R) + // < k * (b_1 * R / n + 1) / R + 1 < b_1 * k / n + 2 < b_1 + 2, so a + // bound like (|b1| + 2) * (|b2| + 2) < 2 * n is good enough for wlog d1 + // < 2 * n let mut d1 = ::BigInt::mul_no_reduce_lo(&c1, Self::B1.as_ref()); if d1 > modulus { @@ -66,7 +69,8 @@ pub trait GLVParameters: Send + Sync + 'static + ModelParameters { if d2 > modulus { d2.sub_noborrow(&modulus); } - // We compute k_2 = -(c1.b1 + c1.b1) = sign(b1)*(c2|b2| - c1|b1|) = sign(b1)(d2 - d1) + // We compute k_2 = -(c1.b1 + c1.b1) = sign(b1)*(c2|b2| - c1|b1|) = sign(b1)(d2 + // - d1) let k2_field = if !Self::B1_IS_NEG { Self::ScalarField::from(d2) - &Self::ScalarField::from(d1) } else { @@ -91,3 +95,28 @@ pub trait GLVParameters: Send + Sync + 'static + ModelParameters { ((neg1, k1), (neg2, k2)) } } + +#[macro_export] +macro_rules! impl_glv_for_sw { + () => { + #[inline(always)] + fn has_glv() -> bool { + true + } + + #[inline(always)] + fn glv_endomorphism_in_place(elem: &mut Self::BaseField) { + *elem *= &::OMEGA; + } + + #[inline] + fn glv_scalar_decomposition( + k: ::BigInt, + ) -> ( + (bool, ::BigInt), + (bool, ::BigInt), + ) { + ::glv_scalar_decomposition_inner(k) + } + }; +} diff --git a/algebra-core/src/curves/mod.rs b/algebra-core/src/curves/mod.rs index 553fd5345..1ba08682d 100644 --- a/algebra-core/src/curves/mod.rs +++ b/algebra-core/src/curves/mod.rs @@ -20,6 +20,7 @@ pub use self::batch_arith::*; pub mod bucketed_add; pub use self::bucketed_add::*; +#[macro_use] pub mod glv; pub use self::glv::*; @@ -202,6 +203,8 @@ pub trait ProjectiveCurve: self = res; self } + + fn get_x(&mut self) -> &mut Self::BaseField; } /// Affine representation of an elliptic curve point guaranteed to be diff --git a/algebra-core/src/curves/models/short_weierstrass_jacobian.rs b/algebra-core/src/curves/models/short_weierstrass_jacobian.rs index db1a1d001..3b06ff835 100644 --- a/algebra-core/src/curves/models/short_weierstrass_jacobian.rs +++ b/algebra-core/src/curves/models/short_weierstrass_jacobian.rs @@ -367,6 +367,10 @@ impl ProjectiveCurve for GroupProjective

{ self } } + + fn get_x(&mut self) -> &mut Self::BaseField { + &mut self.x + } } impl Neg for GroupProjective

{ diff --git a/algebra-core/src/curves/models/short_weierstrass_projective.rs b/algebra-core/src/curves/models/short_weierstrass_projective.rs index 4ed54f0a7..854268ee8 100644 --- a/algebra-core/src/curves/models/short_weierstrass_projective.rs +++ b/algebra-core/src/curves/models/short_weierstrass_projective.rs @@ -141,6 +141,10 @@ impl ProjectiveCurve for GroupProjective

{ type ScalarField = P::ScalarField; type Affine = GroupAffine

; + fn get_x(&mut self) -> &mut Self::BaseField { + &mut self.x + } + #[inline] fn prime_subgroup_generator() -> Self { GroupAffine::prime_subgroup_generator().into() diff --git a/algebra-core/src/curves/models/twisted_edwards_extended.rs b/algebra-core/src/curves/models/twisted_edwards_extended.rs index 1fe50a957..772c5c714 100644 --- a/algebra-core/src/curves/models/twisted_edwards_extended.rs +++ b/algebra-core/src/curves/models/twisted_edwards_extended.rs @@ -703,6 +703,10 @@ impl ProjectiveCurve for GroupProjective

{ // Z3 = F*G self.z = f * &g; } + + fn get_x(&mut self) -> &mut Self::BaseField { + &mut self.x + } } impl Neg for GroupProjective

{ diff --git a/algebra-core/src/fields/macros.rs b/algebra-core/src/fields/macros.rs index b4fbc8232..f6815591a 100644 --- a/algebra-core/src/fields/macros.rs +++ b/algebra-core/src/fields/macros.rs @@ -67,6 +67,10 @@ macro_rules! impl_Fp { impl Field for $Fp

{ type BasePrimeField = Self; + fn extension_degree() -> usize { + 1 + } + #[inline] fn double(&self) -> Self { let mut temp = *self; diff --git a/algebra-core/src/fields/mod.rs b/algebra-core/src/fields/mod.rs index 3ce3c30fd..a317a94ae 100644 --- a/algebra-core/src/fields/mod.rs +++ b/algebra-core/src/fields/mod.rs @@ -103,6 +103,8 @@ pub trait Field: { type BasePrimeField: PrimeField; + fn extension_degree() -> usize; + /// Returns the characteristic of the field. fn characteristic<'a>() -> &'a [u64] { Self::BasePrimeField::characteristic() diff --git a/algebra-core/src/fields/models/cubic_extension.rs b/algebra-core/src/fields/models/cubic_extension.rs index 85f4a6585..b7a78e0ab 100644 --- a/algebra-core/src/fields/models/cubic_extension.rs +++ b/algebra-core/src/fields/models/cubic_extension.rs @@ -139,6 +139,10 @@ impl One for CubicExtField

{ impl Field for CubicExtField

{ type BasePrimeField = P::BasePrimeField; + fn extension_degree() -> usize { + 3 * P::BaseField::extension_degree() + } + fn double(&self) -> Self { let mut result = self.clone(); result.double_in_place(); diff --git a/algebra-core/src/fields/models/quadratic_extension.rs b/algebra-core/src/fields/models/quadratic_extension.rs index 9f9c52def..845e98fb8 100644 --- a/algebra-core/src/fields/models/quadratic_extension.rs +++ b/algebra-core/src/fields/models/quadratic_extension.rs @@ -165,6 +165,10 @@ impl One for QuadExtField

{ impl Field for QuadExtField

{ type BasePrimeField = P::BasePrimeField; + fn extension_degree() -> usize { + 2 * P::BaseField::extension_degree() + } + fn double(&self) -> Self { let mut result = self.clone(); result.double_in_place(); diff --git a/algebra/src/bls12_377/curves/g1.rs b/algebra/src/bls12_377/curves/g1.rs index 801b3b49b..3c318afda 100644 --- a/algebra/src/bls12_377/curves/g1.rs +++ b/algebra/src/bls12_377/curves/g1.rs @@ -1,7 +1,10 @@ use algebra_core::{ - biginteger::{BigInteger256, BigInteger384}, - curves::models::{ModelParameters, SWModelParameters}, - field_new, Zero, + biginteger::{BigInteger256, BigInteger384, BigInteger512}, + curves::{ + models::{ModelParameters, SWModelParameters}, + GLVParameters, + }, + field_new, impl_glv_for_sw, PrimeField, Zero, }; use crate::bls12_377::{Fq, Fr}; @@ -14,6 +17,40 @@ impl ModelParameters for Parameters { type ScalarField = Fr; } +impl GLVParameters for Parameters { + type WideBigInt = BigInteger512; + const OMEGA: Self::BaseField = field_new!( + Fq, + BigInteger384([ + 15766275933608376691, + 15635974902606112666, + 1934946774703877852, + 18129354943882397960, + 15437979634065614942, + 101285514078273488 + ]) + ); + const LAMBDA: Self::ScalarField = field_new!( + Fr, + BigInteger256([ + 12574070832645531618, + 10005695704657941814, + 1564543351912391449, + 657300228442948690 + ]) + ); + /// |round(B1 * R / n)| + const Q2: ::BigInt = + BigInteger256([9183663392111466540, 12968021215939883360, 3, 0]); + const B1: ::BigInt = + BigInteger256([725501752471715841, 4981570305181876225, 0, 0]); + const B1_IS_NEG: bool = false; + /// |round(B2 * R / n)| + const Q1: ::BigInt = BigInteger256([13, 0, 0, 0]); + const B2: ::BigInt = BigInteger256([1, 0, 0, 0]); + const R_BITS: u32 = 256; +} + impl SWModelParameters for Parameters { /// COEFF_A = 0 const COEFF_A: Fq = field_new!(Fq, BigInteger384([0x0, 0x0, 0x0, 0x0, 0x0, 0x0])); @@ -50,6 +87,8 @@ impl SWModelParameters for Parameters { fn mul_by_a(_: &Self::BaseField) -> Self::BaseField { Self::BaseField::zero() } + + impl_glv_for_sw!(); } /// G1_GENERATOR_X = diff --git a/algebra/src/bls12_377/curves/g2.rs b/algebra/src/bls12_377/curves/g2.rs index 98b5040ea..efab698de 100644 --- a/algebra/src/bls12_377/curves/g2.rs +++ b/algebra/src/bls12_377/curves/g2.rs @@ -1,11 +1,13 @@ +use crate::bls12_377::{g1, Fq, Fq2, Fr}; use algebra_core::{ - biginteger::{BigInteger256, BigInteger384}, - curves::models::{ModelParameters, SWModelParameters}, - field_new, Zero, + biginteger::{BigInteger256, BigInteger384, BigInteger512}, + curves::{ + models::{ModelParameters, SWModelParameters}, + GLVParameters, + }, + field_new, impl_glv_for_sw, PrimeField, Zero, }; -use crate::bls12_377::{g1, Fq, Fq2, Fr}; - #[derive(Clone, Default, PartialEq, Eq)] pub struct Parameters; @@ -14,6 +16,44 @@ impl ModelParameters for Parameters { type ScalarField = Fr; } +impl GLVParameters for Parameters { + type WideBigInt = BigInteger512; + const OMEGA: Self::BaseField = field_new!( + Fq2, + field_new!( + Fq, + BigInteger384([ + 3203870859294639911, + 276961138506029237, + 9479726329337356593, + 13645541738420943632, + 7584832609311778094, + 101110569012358506 + ]) + ), + field_new!(Fq, BigInteger384([0, 0, 0, 0, 0, 0])) + ); + const LAMBDA: Self::ScalarField = field_new!( + Fr, + BigInteger256([ + 12574070832645531618, + 10005695704657941814, + 1564543351912391449, + 657300228442948690 + ]) + ); + /// |round(B1 * R / n)| + const Q2: ::BigInt = + BigInteger256([9183663392111466540, 12968021215939883360, 3, 0]); + const B1: ::BigInt = + BigInteger256([725501752471715841, 4981570305181876225, 0, 0]); + const B1_IS_NEG: bool = false; + /// |round(B2 * R / n)| + const Q1: ::BigInt = BigInteger256([13, 0, 0, 0]); + const B2: ::BigInt = BigInteger256([1, 0, 0, 0]); + const R_BITS: u32 = 256; +} + impl SWModelParameters for Parameters { /// COEFF_A = [0, 0] #[rustfmt::skip] @@ -73,6 +113,8 @@ impl SWModelParameters for Parameters { fn mul_by_a(_: &Self::BaseField) -> Self::BaseField { Self::BaseField::zero() } + + impl_glv_for_sw!(); } #[rustfmt::skip] diff --git a/algebra/src/bls12_381/curves/g1.rs b/algebra/src/bls12_381/curves/g1.rs index 65e17283f..b7508f27f 100644 --- a/algebra/src/bls12_381/curves/g1.rs +++ b/algebra/src/bls12_381/curves/g1.rs @@ -1,12 +1,13 @@ use crate::{ - biginteger::{BigInteger256, BigInteger384}, + biginteger::{BigInteger256, BigInteger384, BigInteger512}, bls12_381, bls12_381::*, curves::{ bls12, models::{ModelParameters, SWModelParameters}, + GLVParameters, }, - field_new, Zero, + field_new, impl_glv_for_sw, PrimeField, Zero, }; pub type G1Affine = bls12::G1Affine; @@ -20,6 +21,39 @@ impl ModelParameters for Parameters { type ScalarField = Fr; } +impl GLVParameters for Parameters { + type WideBigInt = BigInteger512; + const OMEGA: Self::BaseField = field_new!( + Fq, + BigInteger384([ + 3526659474838938856, + 17562030475567847978, + 1632777218702014455, + 14009062335050482331, + 3906511377122991214, + 368068849512964448, + ]) + ); + const LAMBDA: Self::ScalarField = field_new!( + Fr, + BigInteger256([ + 7865245318337523249, + 18346590209729131401, + 15545362854776399464, + 6505881510324251116, + ]) + ); + /// |round(B1 * R / n)| + const Q2: ::BigInt = + BigInteger256([7203196592358157870, 8965520006802549469, 1, 0]); + const B1: ::BigInt = + BigInteger256([4294967295, 12413508272118670338, 0, 0]); + const B1_IS_NEG: bool = true; + /// |round(B2 * R / n)| + const Q1: ::BigInt = BigInteger256([2, 0, 0, 0]); + const B2: ::BigInt = BigInteger256([1, 0, 0, 0]); + const R_BITS: u32 = 256; +} impl SWModelParameters for Parameters { /// COEFF_A = 0 const COEFF_A: Fq = field_new!(Fq, BigInteger384([0x0, 0x0, 0x0, 0x0, 0x0, 0x0])); @@ -56,6 +90,8 @@ impl SWModelParameters for Parameters { fn mul_by_a(_: &Self::BaseField) -> Self::BaseField { Self::BaseField::zero() } + + impl_glv_for_sw!(); } /// G1_GENERATOR_X = diff --git a/algebra/src/bls12_381/curves/g2.rs b/algebra/src/bls12_381/curves/g2.rs index 65ba55d67..a851d53e0 100644 --- a/algebra/src/bls12_381/curves/g2.rs +++ b/algebra/src/bls12_381/curves/g2.rs @@ -1,12 +1,13 @@ use crate::{ - biginteger::{BigInteger256, BigInteger384}, + biginteger::{BigInteger256, BigInteger384, BigInteger512}, bls12_381, bls12_381::*, curves::{ bls12, models::{ModelParameters, SWModelParameters}, + GLVParameters, }, - field_new, Zero, + field_new, impl_glv_for_sw, PrimeField, Zero, }; pub type G2Affine = bls12::G2Affine; @@ -20,6 +21,44 @@ impl ModelParameters for Parameters { type ScalarField = Fr; } +impl GLVParameters for Parameters { + type WideBigInt = BigInteger512; + const OMEGA: Self::BaseField = field_new!( + Fq2, + field_new!( + Fq, + BigInteger384([ + 14772873186050699377, + 6749526151121446354, + 6372666795664677781, + 10283423008382700446, + 286397964926079186, + 1796971870900422465 + ]) + ), + field_new!(Fq, BigInteger384([0, 0, 0, 0, 0, 0])) + ); + const LAMBDA: Self::ScalarField = field_new!( + Fr, + BigInteger256([ + 7865245318337523249, + 18346590209729131401, + 15545362854776399464, + 6505881510324251116 + ]) + ); + /// |round(B1 * R / n)| + const Q2: ::BigInt = + BigInteger256([7203196592358157870, 8965520006802549469, 1, 0]); + const B1: ::BigInt = + BigInteger256([4294967295, 12413508272118670338, 0, 0]); + const B1_IS_NEG: bool = true; + /// |round(B2 * R / n)| + const Q1: ::BigInt = BigInteger256([2, 0, 0, 0]); + const B2: ::BigInt = BigInteger256([1, 0, 0, 0]); + const R_BITS: u32 = 256; +} + impl SWModelParameters for Parameters { /// COEFF_A = [0, 0] const COEFF_A: Fq2 = field_new!(Fq2, g1::Parameters::COEFF_A, g1::Parameters::COEFF_A,); @@ -60,6 +99,8 @@ impl SWModelParameters for Parameters { fn mul_by_a(_: &Self::BaseField) -> Self::BaseField { Self::BaseField::zero() } + + impl_glv_for_sw!(); } pub const G2_GENERATOR_X: Fq2 = field_new!(Fq2, G2_GENERATOR_X_C0, G2_GENERATOR_X_C1); diff --git a/algebra/src/bn254/curves/g1.rs b/algebra/src/bn254/curves/g1.rs index 8f0a81952..b9b59ce23 100644 --- a/algebra/src/bn254/curves/g1.rs +++ b/algebra/src/bn254/curves/g1.rs @@ -1,7 +1,7 @@ use algebra_core::{ - biginteger::BigInteger256, + biginteger::{BigInteger256, BigInteger512}, curves::models::{ModelParameters, SWModelParameters}, - field_new, Zero, + field_new, impl_glv_for_sw, GLVParameters, PrimeField, Zero, }; use crate::bn254::{Fq, Fr}; @@ -14,6 +14,41 @@ impl ModelParameters for Parameters { type ScalarField = Fr; } +impl GLVParameters for Parameters { + type WideBigInt = BigInteger512; + const OMEGA: Self::BaseField = field_new!( + Fq, + BigInteger256([ + 3697675806616062876, + 9065277094688085689, + 6918009208039626314, + 2775033306905974752 + ]) + ); + const LAMBDA: Self::ScalarField = field_new!( + Fr, + BigInteger256([ + 244305545194690131, + 8351807910065594880, + 14266533074055306532, + 404339206190769364 + ]) + ); + /// |round(B1 * R / n)| + /// |round(B1 * R / n)| + const Q2: ::BigInt = + BigInteger256([6023842690951505253, 5534624963584316114, 2, 0]); + const B1: ::BigInt = + BigInteger256([857057580304901387, 8020209761171036669, 0, 0]); + const B1_IS_NEG: bool = false; + /// |round(B2 * R / n)| + const Q1: ::BigInt = + BigInteger256([15644699364383830999, 2, 0, 0]); + const B2: ::BigInt = + BigInteger256([9931322734385697763, 0, 0, 0]); + const R_BITS: u32 = 256; +} + impl SWModelParameters for Parameters { /// COEFF_A = 0 const COEFF_A: Fq = field_new!(Fq, BigInteger256([0x0, 0x0, 0x0, 0x0])); @@ -47,6 +82,8 @@ impl SWModelParameters for Parameters { fn mul_by_a(_: &Self::BaseField) -> Self::BaseField { Self::BaseField::zero() } + + impl_glv_for_sw!(); } /// G1_GENERATOR_X = diff --git a/algebra/src/bn254/curves/g2.rs b/algebra/src/bn254/curves/g2.rs index e85a89968..d4c51e6f5 100644 --- a/algebra/src/bn254/curves/g2.rs +++ b/algebra/src/bn254/curves/g2.rs @@ -1,7 +1,7 @@ use algebra_core::{ - biginteger::BigInteger256, + biginteger::{BigInteger256, BigInteger512}, curves::models::{ModelParameters, SWModelParameters}, - field_new, Zero, + field_new, impl_glv_for_sw, GLVParameters, PrimeField, Zero, }; use crate::bn254::{g1, Fq, Fq2, Fr}; @@ -14,6 +14,44 @@ impl ModelParameters for Parameters { type ScalarField = Fr; } +impl GLVParameters for Parameters { + type WideBigInt = BigInteger512; + const OMEGA: Self::BaseField = field_new!( + Fq2, + field_new!( + Fq, + BigInteger256([ + 8183898218631979349, + 12014359695528440611, + 12263358156045030468, + 3187210487005268291 + ]) + ), + field_new!(Fq, BigInteger256([0, 0, 0, 0])), + ); + const LAMBDA: Self::ScalarField = field_new!( + Fr, + BigInteger256([ + 244305545194690131, + 8351807910065594880, + 14266533074055306532, + 404339206190769364 + ]) + ); + /// |round(B1 * R / n)| + const Q2: ::BigInt = + BigInteger256([6023842690951505253, 5534624963584316114, 2, 0]); + const B1: ::BigInt = + BigInteger256([857057580304901387, 8020209761171036669, 0, 0]); + const B1_IS_NEG: bool = false; + /// |round(B2 * R / n)| + const Q1: ::BigInt = + BigInteger256([15644699364383830999, 2, 0, 0]); + const B2: ::BigInt = + BigInteger256([9931322734385697763, 0, 0, 0]); + const R_BITS: u32 = 256; +} + impl SWModelParameters for Parameters { /// COEFF_A = [0, 0] #[rustfmt::skip] @@ -68,6 +106,8 @@ impl SWModelParameters for Parameters { fn mul_by_a(_: &Self::BaseField) -> Self::BaseField { Self::BaseField::zero() } + + impl_glv_for_sw!(); } #[rustfmt::skip] diff --git a/algebra/src/bw6_761/curves/g1.rs b/algebra/src/bw6_761/curves/g1.rs index 778171755..a6512199e 100644 --- a/algebra/src/bw6_761/curves/g1.rs +++ b/algebra/src/bw6_761/curves/g1.rs @@ -8,6 +8,7 @@ use crate::{ }, field_new, fields::PrimeField, + impl_glv_for_sw, }; pub type G1Affine = GroupAffine; @@ -68,34 +69,34 @@ impl GLVParameters for Parameters { ); /// |round(B1 * R / n)| const Q2: ::BigInt = BigInteger384([ - 14430678704534329733, - 14479735877321354361, - 6958676793196883088, - 21, + 11941976086484053770, + 4826578625773784813, + 2319558931065627696, + 7, 0, 0, ]); const B1: ::BigInt = BigInteger384([ - 9586122913090633729, - 9963140610363752448, - 2588746559005780992, + 6390748608727089153, + 3321046870121250816, + 862915519668593664, 0, 0, 0, ]); - const B1_IS_NEG: bool = false; + const B1_IS_NEG: bool = true; /// |round(B2 * R / n)| const Q1: ::BigInt = BigInteger384([ - 11941976086484053770, - 4826578625773784813, + 8993470605275773807, + 4826578625773784734, 2319558931065627696, 7, 0, 0, ]); const B2: ::BigInt = BigInteger384([ - 6390748608727089153, - 3321046870121250816, + 15251369769346007039, + 3321046870121250815, 862915519668593664, 0, 0, @@ -160,25 +161,7 @@ impl SWModelParameters for Parameters { Self::BaseField::zero() } - #[inline(always)] - fn has_glv() -> bool { - true - } - - #[inline(always)] - fn glv_endomorphism_in_place(elem: &mut Self::BaseField) { - *elem *= &::OMEGA; - } - - #[inline] - fn glv_scalar_decomposition( - k: ::BigInt, - ) -> ( - (bool, ::BigInt), - (bool, ::BigInt), - ) { - ::glv_scalar_decomposition_inner(k) - } + impl_glv_for_sw!(); } /// G1_GENERATOR_X = diff --git a/algebra/src/bw6_761/curves/g2.rs b/algebra/src/bw6_761/curves/g2.rs index 88be2a2fc..a3d363067 100644 --- a/algebra/src/bw6_761/curves/g2.rs +++ b/algebra/src/bw6_761/curves/g2.rs @@ -8,6 +8,7 @@ use crate::{ }, field_new, fields::PrimeField, + impl_glv_for_sw, }; pub type G2Affine = GroupAffine; @@ -61,34 +62,34 @@ impl GLVParameters for Parameters { ); /// |round(B1 * R / n)| const Q2: ::BigInt = BigInteger384([ - 14430678704534329733, - 14479735877321354361, - 6958676793196883088, - 21, + 11941976086484053770, + 4826578625773784813, + 2319558931065627696, + 7, 0, 0, ]); const B1: ::BigInt = BigInteger384([ - 9586122913090633729, - 9963140610363752448, - 2588746559005780992, + 6390748608727089153, + 3321046870121250816, + 862915519668593664, 0, 0, 0, ]); - const B1_IS_NEG: bool = false; + const B1_IS_NEG: bool = true; /// |round(B2 * R / n)| const Q1: ::BigInt = BigInteger384([ - 11941976086484053770, - 4826578625773784813, + 8993470605275773807, + 4826578625773784734, 2319558931065627696, 7, 0, 0, ]); const B2: ::BigInt = BigInteger384([ - 6390748608727089153, - 3321046870121250816, + 15251369769346007039, + 3321046870121250815, 862915519668593664, 0, 0, @@ -153,25 +154,7 @@ impl SWModelParameters for Parameters { Self::BaseField::zero() } - #[inline(always)] - fn has_glv() -> bool { - true - } - - #[inline(always)] - fn glv_endomorphism_in_place(elem: &mut Self::BaseField) { - *elem *= &::OMEGA; - } - - #[inline] - fn glv_scalar_decomposition( - k: ::BigInt, - ) -> ( - (bool, ::BigInt), - (bool, ::BigInt), - ) { - ::glv_scalar_decomposition_inner(k) - } + impl_glv_for_sw!(); } /// G2_GENERATOR_X = diff --git a/algebra/src/tests/curves.rs b/algebra/src/tests/curves.rs index cd642f82b..f51cdc49f 100644 --- a/algebra/src/tests/curves.rs +++ b/algebra/src/tests/curves.rs @@ -136,6 +136,29 @@ fn random_multiplication_test() { } } +fn glv_multiplication_test() { + use algebra_core::fields::BitIteratorBE; + let mut rng = XorShiftRng::seed_from_u64(1231275789u64); + + for _ in 0..ITERATIONS { + let mut a = G::rand(&mut rng); + let b = a; + + let s = G::ScalarField::rand(&mut rng); + a.mul_assign(s); + + let mut res = G::zero(); + for bit in BitIteratorBE::without_leading_zeros(s.into()) { + res.double_in_place(); + if bit { + res += b; + } + } + + assert_eq!(a, res); + } +} + fn random_doubling_test() { let mut rng = XorShiftRng::seed_from_u64(1231275789u64); @@ -655,6 +678,7 @@ pub fn curve_tests() { random_doubling_test::(); random_negation_test::(); random_transformation_test::(); + glv_multiplication_test::(); } pub fn batch_affine_test() { diff --git a/scripts/glv_lattice_basis/examples/main.rs b/scripts/glv_lattice_basis/examples/main.rs index a4ba58203..536e67f0f 100644 --- a/scripts/glv_lattice_basis/examples/main.rs +++ b/scripts/glv_lattice_basis/examples/main.rs @@ -1,107 +1,10 @@ extern crate algebra; extern crate algebra_core; -extern crate num_traits; -use algebra::bw6_761::{Fq, Fr}; -use algebra_core::{ - biginteger::{BigInteger, BigInteger384, BigInteger768}, - fields::PrimeField, -}; +use algebra::bw6_761::G1Projective as GroupProjective; +use algebra_core::{BigInteger768 as BaseFieldBigInt, BigInteger768 as FrWideBigInt}; use glv_lattice_basis::*; -use num_traits::Zero; -use std::ops::Neg; fn main() { - let _omega_g1 = BigInteger768([ - 0x962140000000002a, - 0xc547ba8a4000002f, - 0xb6290012d96f8819, - 0xf2f082d4dcb5e37c, - 0xc65759fc45183151, - 0x8e0a235a0a398300, - 0xab5e57926fa70184, - 0xee4a737f73b6f952, - 0x2d17be416c5e4426, - 0x6c1f31e53bd9603c, - 0xaa846c61024e4cca, - 0x531dc16c6ecd27, - ]); - let _omega_g2 = BigInteger768([ - 0x5e7bc00000000060, - 0x214983de30000053, - 0x5fe3f89c11811c1e, - 0xa5b093ed79b1c57b, - 0xab8579e02ed3cddc, - 0xf87fa59308c07a8f, - 0x5870636cb60d217f, - 0x823132b971cdefc6, - 0x256ab7ae14297a1a, - 0x4d06e68545f7e64c, - 0x27035cdf02acb274, - 0xcfca638f1500e3, - ]); - println!( - "const OMEGA: Self::BaseField = {:?};", - Fq::from_repr(_omega_g2).unwrap() - ); - let n = BigInteger384([ - 0x8508c00000000001, - 0x170b5d4430000000, - 0x1ef3622fba094800, - 0x1a22d9f300f5138f, - 0xc63b05c06ca1493b, - 0x1ae3a4617c510ea, - ]); - let lambda = BigInteger384([ - 0x8508c00000000001, - 0x452217cc90000000, - 0xc5ed1347970dec00, - 0x619aaf7d34594aab, - 0x9b3af05dd14f6ec, - 0x0, - ]); - println!( - "const LAMBDA: Self::ScalarField = {:?};", - Fr::from_repr(lambda).unwrap() - ); - - let vecs = get_lattice_basis::(n, lambda); - - for (i, vec) in [vecs.0, vecs.1].iter().enumerate() { - // println!("vec: {:?}", vec); - let (s1, (flag, t1)) = vec; - - let mut t1_big = BigInteger768::from_slice(t1.as_ref()); - let n_big = BigInteger768::from_slice(n.as_ref()); - t1_big.muln(BigInteger384::NUM_LIMBS as u32 * 64); - let (g1_big, _) = div_with_remainder::(t1_big, n_big); - let g1 = BigInteger384::from_slice(g1_big.as_ref()); - - println!("/// |round(B{} * R / n)|", i + 1); - println!( - "const Q{}: ::BigInt = {:?};", - ((i + 1) % 2) + 1, - g1 - ); - println!( - "const B{}: ::BigInt = {:?};", - i + 1, - t1 - ); - println!("const B{}_IS_NEG: bool = {:?};", i + 1, flag); - - debug_assert_eq!( - recompose_integer( - Fr::from_repr(*s1).unwrap(), - if !flag { - Fr::from_repr(*t1).unwrap() - } else { - Fr::from_repr(*t1).unwrap().neg() - }, - Fr::from_repr(lambda).unwrap() - ), - Fr::zero() - ); - } - println!("const R_BITS: u32 = {:?};", BigInteger384::NUM_LIMBS * 64); + print_glv_params::(); } diff --git a/scripts/glv_lattice_basis/src/lib.rs b/scripts/glv_lattice_basis/src/lib.rs index 2cf1a291c..172bc230e 100644 --- a/scripts/glv_lattice_basis/src/lib.rs +++ b/scripts/glv_lattice_basis/src/lib.rs @@ -4,9 +4,166 @@ extern crate num_traits; mod arithmetic; -use algebra::BigInteger; -use algebra_core::fields::PrimeField; +use algebra_core::{BigInteger, Field, PrimeField, ProjectiveCurve}; pub use arithmetic::*; +use num_traits::Zero; +use std::ops::Neg; + +/// Takes data from two endomorphisms and sorts out which corresponds to which +fn which_endo( + base_roots: (G::BaseField, G::BaseField), + scalar_roots: (G::ScalarField, G::ScalarField), +) -> ( + (G::BaseField, G::ScalarField), + (G::BaseField, G::ScalarField), +) { + // println!("{:?}, {:?}", base_roots, scalar_roots); + let g = G::prime_subgroup_generator(); + + let mut g_endo = g; + *g_endo.get_x() *= &base_roots.0; + + let d1 = if g.mul(scalar_roots.0) == g_endo { + (base_roots.0, scalar_roots.0) + } else { + let mut g_endo = g; + *g_endo.get_x() *= &base_roots.1; + assert!(g.mul(scalar_roots.0) == g_endo); + + (base_roots.1, scalar_roots.0) + }; + + let d2 = if g.mul(scalar_roots.1) == g_endo { + (base_roots.0, scalar_roots.1) + } else { + let mut g_endo = g; + *g_endo.get_x() *= &base_roots.1; + assert!(g.mul(scalar_roots.1) == g_endo); + + (base_roots.1, scalar_roots.1) + }; + + (d1, d2) +} + +fn cube_root_unity() -> (F, F) { + let char = B::from_slice(F::characteristic()); + let deg = F::extension_degree(); + let mut modulus = char; + for _ in 1..deg { + modulus = B::mul_no_reduce_lo(&modulus.as_ref(), &char.as_ref()); + } + + modulus.sub_noborrow(&B::from(1)); + let (q, r) = div_with_remainder(modulus, B::from(3)); + assert!(r == B::from(0)); + + let mut g = 2u32; + let mut root1 = F::one(); + loop { + if root1 != F::one() { + break; + } + let x = F::from(g); + root1 = x.pow(q); + g += 1; + } + let root2 = root1 * root1; + assert!(root1.pow(&[3]) == F::one()); + assert!(root2.pow(&[3]) == F::one()); + assert!(root1 != root2); + + (root1, root2) +} + +fn get_endo_data() -> (G::BaseField, G::ScalarField) { + // println!("{:?}", which_endo::( + // cube_root_unity::(), + // cube_root_unity::(), + // )); + which_endo::( + cube_root_unity::(), + cube_root_unity::::BigInt>(), + ) + .1 +} + +pub fn print_glv_params() { + let (omega, lambda) = get_endo_data::(); + let g = G::prime_subgroup_generator(); + let mut g_endo = g; + *g_endo.get_x() *= ω + assert!(g.mul(lambda) == g_endo); + + println!("const OMEGA: Self::BaseField = {:?};", omega); + let n = ::modulus(); + println!("const LAMBDA: Self::ScalarField = {:?};", lambda); + + let vecs = get_lattice_basis::(n, lambda.into_repr()); + + // We check that `(|B1| + 2) * (|B2| + 2) < 2n` + // and `B_i^2 < 2n` e.g. `|B_i| < \sqrt{2n}$ + // We use this to prove some bounds later + let wide_modulus = WideBigInt::from_slice(&n.as_ref()[..]); + let two_modulus = WideBigInt::mul_no_reduce_lo( + &wide_modulus.as_ref()[..], + &WideBigInt::from(2).as_ref()[..], + ); + + let mut b1 = ((vecs.0).1).1; + let mut b2 = ((vecs.1).1).1; + let two = ::BigInt::from(2); + let b1b1 = WideBigInt::mul_no_reduce(&b1.as_ref()[..], &b1.as_ref()[..]); + let b2b2 = WideBigInt::mul_no_reduce(&b2.as_ref()[..], &b2.as_ref()[..]); + + b1.add_nocarry(&two); + b2.add_nocarry(&two); + let b1b2 = WideBigInt::mul_no_reduce(&b1.as_ref()[..], &b2.as_ref()[..]); + + assert!(b1b1 < two_modulus); + assert!(b2b2 < two_modulus); + assert!(b1b2 < two_modulus); + + for (i, vec) in [vecs.0, vecs.1].iter().enumerate() { + let (s1, (flag, t1)) = vec; + + let mut t1_big = WideBigInt::from_slice(t1.as_ref()); + let n_big = WideBigInt::from_slice(n.as_ref()); + t1_big.muln(::BigInt::NUM_LIMBS as u32 * 64); + let (g1_big, _) = div_with_remainder::(t1_big, n_big); + let g1 = ::BigInt::from_slice(g1_big.as_ref()); + + println!("/// |round(B{} * R / n)|", i + 1); + println!( + "const Q{}: ::BigInt = {:?};", + ((i + 1) % 2) + 1, + g1 + ); + println!( + "const B{}: ::BigInt = {:?};", + i + 1, + t1 + ); + println!("const B{}_IS_NEG: bool = {:?};", i + 1, flag); + + debug_assert_eq!( + recompose_integer( + G::ScalarField::from_repr(*s1).unwrap(), + if !flag { + G::ScalarField::from_repr(*t1).unwrap() + } else { + G::ScalarField::from_repr(*t1).unwrap().neg() + }, + lambda + ), + G::ScalarField::zero() + ); + } + println!( + "const R_BITS: u32 = {:?};", + ::BigInt::NUM_LIMBS * 64 + ); +} // We work on arrays of size 3 // We assume that |E(F_q)| < R = 2^{ceil(limbs/2) * 64} @@ -23,16 +180,19 @@ pub fn get_lattice_basis( let mut t: [F; 3] = [zero, one, zero]; let max_num_bits_lattice = (F::BigInt::from_slice(F::characteristic()).num_bits() - 1) / 2 + 1; + // We can use an approximation as we are merely using a heuristic. We should + // check that the parameters obtained from this heuristic satisfies the + // required conditions separately. let sqrt_n = as_f64(n.as_ref()).sqrt(); - println!("Log sqrtn: {}", sqrt_n.log2()); + // println!("Log sqrtn: {}", sqrt_n.log2()); let mut i = 0; // While r_i >= sqrt(n), we perform the extended euclidean algorithm so that // si*n + ti*lambda = ri then return the vectors (r_i, (sign(t_i), |t_i|)), // (r_i+1, (sign(t_i+1), |t_i+1|)) Notice this makes ri + (-ti)*lambda = 0 // mod n, which is what we desire for our short lattice basis - while as_f64(r[i % 3].as_ref()) >= sqrt_n { + while as_f64(r[(i + 1) % 3].as_ref()) >= sqrt_n { // while i < 20 { let (q, rem): (F::BigInt, F::BigInt) = div_with_remainder::(r[i % 3], r[(i + 1) % 3]);