Implement ABGLSV-Pornin multiplication

.gitattributes

@@ -0,0 +1,4 @@
+curve25519-dalek/src/backend/serial/u32/constants/affine_odd_multiples_of_b_shl_128.rs linguist-generated=true
+curve25519-dalek/src/backend/serial/u64/constants/affine_odd_multiples_of_b_shl_128.rs linguist-generated=true
+curve25519-dalek/src/backend/vector/avx2/constants/b_shl_128_odd_lookup_table.rs linguist-generated=true
+curve25519-dalek/src/backend/vector/ifma/constants/b_shl_128_odd_lookup_table.rs linguist-generated=true

.github/workflows/workspace.yml

@@ -90,6 +90,26 @@ jobs:
         components: clippy
     - run: cargo clippy --target x86_64-unknown-linux-gnu --all-features
 
+  generated:
+    name: Check generated tables
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        bits: ['32', '64']
+    steps:
+    - uses: actions/checkout@v3
+    - uses: dtolnay/rust-toolchain@stable
+      with:
+        components: rustfmt
+    - run: cargo test table_generators
+      env:
+        RUSTFLAGS: '--cfg curve25519_dalek_bits="${{ matrix.bits }}" --cfg curve25519_dalek_generate_tables'
+    - name: Verify working directory is not clean
+      run: if $(git diff --quiet); then false; else true; fi
+    - run: cargo fmt
+    - name: Verify working directory is clean
+      run: git diff --exit-code
+
   rustfmt:
     name: Check formatting
     runs-on: ubuntu-latest

curve25519-dalek/CHANGELOG.md

@@ -5,6 +5,11 @@ major series.
 
 ## 4.x series
 
+### Unreleased
+
+* Add `EdwardsPoint::vartime_check_double_scalar_mul_basepoint`.
+* Add `RistrettoPoint::vartime_check_double_scalar_mul_basepoint`.
+
 ### 4.1.2
 
 * Fix nightly SIMD build

curve25519-dalek/benches/dalek_benchmarks.rs

@@ -56,6 +56,24 @@ mod edwards_benches {
         });
     }
 
+    fn vartime_check_double_scalar_mul_basepoint<M: Measurement>(c: &mut BenchmarkGroup<M>) {
+        c.bench_function(
+            "Variable-time 8(aA+bB)=8C, A&C variable, B fixed",
+            |bench| {
+                let mut rng = thread_rng();
+                let A = &Scalar::random(&mut rng) * constants::ED25519_BASEPOINT_TABLE;
+                let C = &Scalar::random(&mut rng) * constants::ED25519_BASEPOINT_TABLE;
+                bench.iter_batched(
+                    || (Scalar::random(&mut rng), Scalar::random(&mut rng)),
+                    |(a, b)| {
+                        EdwardsPoint::vartime_check_double_scalar_mul_basepoint(&a, &A, &b, &C)
+                    },
+                    BatchSize::SmallInput,
+                );
+            },
+        );
+    }
+
     pub(crate) fn edwards_benches() {
         let mut c = Criterion::default();
         let mut g = c.benchmark_group("edwards benches");
@@ -65,6 +83,7 @@ mod edwards_benches {
         consttime_fixed_base_scalar_mul(&mut g);
         consttime_variable_base_scalar_mul(&mut g);
         vartime_double_base_scalar_mul(&mut g);
+        vartime_check_double_scalar_mul_basepoint(&mut g);
     }
 }
 

curve25519-dalek/src/backend/mod.rs

@@ -249,3 +249,29 @@ pub fn vartime_double_base_mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar) -> Edwa
         BackendKind::Serial => serial::scalar_mul::vartime_double_base::mul(a, A, b),
     }
 }
+
+/// Computes \\([δa]A + [δb]B - [δ]C\\) in variable time.
+///
+/// - \\(B\\) is the Ed25519 basepoint.
+/// - \\(δ\\) is a value invertible \\( \mod \ell \\), which is selected internally to
+///   this function.
+///
+/// This corresponds to the signature verification optimisation presented in
+/// [Antipa et al 2005](http://cacr.uwaterloo.ca/techreports/2005/cacr2005-28.pdf).
+#[allow(non_snake_case)]
+pub fn scalar_mul_abglsv_pornin(
+    a: &Scalar,
+    A: &EdwardsPoint,
+    b: &Scalar,
+    C: &EdwardsPoint,
+) -> EdwardsPoint {
+    match get_selected_backend() {
+        #[cfg(curve25519_dalek_backend = "simd")]
+        BackendKind::Avx2 => vector::scalar_mul::abglsv_pornin::spec_avx2::mul(a, A, b, C),
+        #[cfg(all(curve25519_dalek_backend = "simd", nightly))]
+        BackendKind::Avx512 => {
+            vector::scalar_mul::abglsv_pornin::spec_avx512ifma_avx512vl::mul(a, A, b, C)
+        }
+        BackendKind::Serial => serial::scalar_mul::abglsv_pornin::mul(a, A, b, C),
+    }
+}

curve25519-dalek/src/backend/serial/fiat_u32/mod.rs

@@ -22,5 +22,5 @@ pub mod scalar;
 
 pub mod field;
 
-#[path = "../u32/constants.rs"]
+#[path = "../u32/constants/mod.rs"]
 pub mod constants;

curve25519-dalek/src/backend/serial/fiat_u64/mod.rs

@@ -24,5 +24,5 @@ pub mod scalar;
 
 pub mod field;
 
-#[path = "../u64/constants.rs"]
+#[path = "../u64/constants/mod.rs"]
 pub mod constants;

curve25519-dalek/src/backend/serial/scalar_mul/abglsv_pornin.rs

@@ -0,0 +1,192 @@
+// -*- mode: rust; -*-
+//
+// This file is part of curve25519-dalek.
+// Copyright (c) 2020-2024 Jack Grigg
+// See LICENSE for licensing information.
+//
+// Author:
+// Jack Grigg <thestr4d@gmail.com>
+#![allow(non_snake_case)]
+
+use core::cmp::Ordering;
+
+use crate::{
+    backend::serial::curve_models::{ProjectiveNielsPoint, ProjectivePoint},
+    constants,
+    edwards::EdwardsPoint,
+    scalar::{lattice_reduction::find_short_vector, Scalar},
+    traits::Identity,
+    window::NafLookupTable5,
+};
+
+/// Computes \\([δa]A + [δb]B - [δ]C\\) in variable time.
+///
+/// - \\(B\\) is the Ed25519 basepoint.
+/// - \\(δ\\) is a value invertible \\( \mod \ell \\), which is selected internally to
+///   this function.
+///
+/// This corresponds to the signature verification optimisation presented in
+/// [Antipa et al 2005](http://cacr.uwaterloo.ca/techreports/2005/cacr2005-28.pdf).
+pub fn mul(a: &Scalar, A: &EdwardsPoint, b: &Scalar, C: &EdwardsPoint) -> EdwardsPoint {
+    // Starting with the target equation:
+    //
+    //     [(δa mod l)]A + [(δb mod l)]B - [δ]C
+    //
+    // We can split δb mod l into two halves e_0 (128 bits) and e_1 (125 bits), and
+    // rewrite the equation as:
+    //
+    //     [(δa mod l)]A + [e_0]B + [e_1 2^128]B - [δ]C
+    //
+    // B and [2^128]B are precomputed, and their resulting scalar multiplications each
+    // have half as many doublings. We therefore want to find a pair of signed integers
+    //
+    //     (d_0, d_1) = (δa mod l, δ)
+    //
+    // that both have as few bits as possible, similarly reducing the number of doublings
+    // in the scalar multiplications [d_0]A and [d_1]C. This is equivalent to finding a
+    // short vector in a lattice of dimension 2.
+
+    // Find a short vector.
+    let (d_0, d_1) = find_short_vector(a);
+
+    // Move the signs of d_0 and d_1 into their corresponding bases and scalars.
+    let A = if d_0.is_negative() { -A } else { *A };
+    let (b, negC) = if d_1.is_negative() {
+        (-b, *C)
+    } else {
+        (*b, -C)
+    };
+    let d_0 = Scalar::from(d_0.unsigned_abs());
+    let d_1 = Scalar::from(d_1.unsigned_abs());
+
+    // Calculate the remaining scalars.
+    let (e_0, e_1) = {
+        let db = b * d_1;
+        let mut e_0 = [0; 32];
+        let mut e_1 = [0; 32];
+        e_0[..16].copy_from_slice(&db.as_bytes()[..16]);
+        e_1[..16].copy_from_slice(&db.as_bytes()[16..]);
+        (Scalar { bytes: e_0 }, Scalar { bytes: e_1 })
+    };
+
+    // Now we can compute the following using Straus's method:
+    //     [d_0]A + [e_0]B + [e_1][2^128]B + [d_1][-C]
+    //
+    // We inline it here so we can use precomputed multiples of [2^128]B.
+
+    let d_0_naf = d_0.non_adjacent_form(5);
+
+    #[cfg(feature = "precomputed-tables")]
+    let e_0_naf = e_0.non_adjacent_form(8);
+    #[cfg(not(feature = "precomputed-tables"))]
+    let e_0_naf = e_0.non_adjacent_form(5);
+
+    #[cfg(feature = "precomputed-tables")]
+    let e_1_naf = e_1.non_adjacent_form(8);
+    #[cfg(not(feature = "precomputed-tables"))]
+    let e_1_naf = e_1.non_adjacent_form(5);
+
+    let d_1_naf = d_1.non_adjacent_form(5);
+
+    // Find starting index
+    let mut i: usize = 255;
+    for j in (0..256).rev() {
+        i = j;
+        if d_0_naf[i] != 0 || e_0_naf[i] != 0 || e_1_naf[i] != 0 || d_1_naf[i] != 0 {
+            break;
+        }
+    }
+
+    let table_A = NafLookupTable5::<ProjectiveNielsPoint>::from(&A);
+
+    #[cfg(feature = "precomputed-tables")]
+    let table_B = &constants::AFFINE_ODD_MULTIPLES_OF_BASEPOINT;
+    #[cfg(not(feature = "precomputed-tables"))]
+    let table_B =
+        &NafLookupTable5::<ProjectiveNielsPoint>::from(&constants::ED25519_BASEPOINT_POINT);
+
+    #[cfg(feature = "precomputed-tables")]
+    let table_B_SHL_128 = &constants::AFFINE_ODD_MULTIPLES_OF_B_SHL_128;
+    #[cfg(not(feature = "precomputed-tables"))]
+    let table_B_SHL_128 =
+        &NafLookupTable5::<ProjectiveNielsPoint>::from(&constants::ED25519_BASEPOINT_SHL_128);
+
+    let table_negC = NafLookupTable5::<ProjectiveNielsPoint>::from(&negC);
+
+    let mut r = ProjectivePoint::identity();
+    loop {
+        let mut t = r.double();
+
+        match d_0_naf[i].cmp(&0) {
+            Ordering::Greater => t = &t.as_extended() + &table_A.select(d_0_naf[i] as usize),
+            Ordering::Less => t = &t.as_extended() - &table_A.select(-d_0_naf[i] as usize),
+            Ordering::Equal => {}
+        }
+
+        match e_0_naf[i].cmp(&0) {
+            Ordering::Greater => t = &t.as_extended() + &table_B.select(e_0_naf[i] as usize),
+            Ordering::Less => t = &t.as_extended() - &table_B.select(-e_0_naf[i] as usize),
+            Ordering::Equal => {}
+        }
+
+        match e_1_naf[i].cmp(&0) {
+            Ordering::Greater => {
+                t = &t.as_extended() + &table_B_SHL_128.select(e_1_naf[i] as usize);
+            }
+            Ordering::Less => t = &t.as_extended() - &table_B_SHL_128.select(-e_1_naf[i] as usize),
+            Ordering::Equal => {}
+        }
+
+        match d_1_naf[i].cmp(&0) {
+            Ordering::Greater => t = &t.as_extended() + &table_negC.select(d_1_naf[i] as usize),
+            Ordering::Less => t = &t.as_extended() - &table_negC.select(-d_1_naf[i] as usize),
+            Ordering::Equal => {}
+        }
+
+        r = t.as_projective();
+
+        if i == 0 {
+            break;
+        }
+        i -= 1;
+    }
+
+    r.as_extended()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::mul;
+    use crate::{constants::ED25519_BASEPOINT_POINT, scalar::Scalar, traits::IsIdentity};
+
+    #[test]
+    fn test_mul() {
+        let a = Scalar::from(2u8);
+        let A = ED25519_BASEPOINT_POINT.double(); // [2]B
+        let b = Scalar::from(4u8);
+        let C = A.double().double(); // [8]B
+
+        // The equation evaluates to the identity, so will be unaffected by δ.
+        assert_eq!(
+            mul(&a, &A, &b, &C),
+            (a * A) + (b * ED25519_BASEPOINT_POINT) - C
+        );
+
+        // Now test some random values.
+        let mut rng = rand::thread_rng();
+
+        for _ in 0..100 {
+            let a = Scalar::random(&mut rng);
+            let A = &Scalar::random(&mut rng) * ED25519_BASEPOINT_POINT;
+            let b = Scalar::random(&mut rng);
+
+            // With a correctly-constructed C, we get the identity.
+            let C = (a * A) + (b * ED25519_BASEPOINT_POINT);
+            assert!(mul(&a, &A, &b, &C).is_identity());
+
+            // With a random C, with high probability we do not get the identity.
+            let C = &Scalar::random(&mut rng) * ED25519_BASEPOINT_POINT;
+            assert!(!mul(&a, &A, &b, &C).is_identity());
+        }
+    }
+}

curve25519-dalek/src/backend/serial/scalar_mul/mod.rs

@@ -17,6 +17,8 @@
 //! scalar multiplication implementations, since it only uses one
 //! curve model.
 
+pub mod abglsv_pornin;
+
 #[allow(missing_docs)]
 pub mod variable_base;