Optimise class group impl

fastcrypto/benches/groups.rs

@@ -185,22 +185,29 @@ mod group_benches {
             self.0 * scalar
         }
     }
-
     fn class_group_ops(c: &mut Criterion) {
-        let mut group: BenchmarkGroup<_> = c.benchmark_group("Class Group Operation");
+        let mut group: BenchmarkGroup<_> = c.benchmark_group("Class Group");
         let d = Discriminant::try_from(BigInt::from_str_radix("-9458193260787340859710210783898414376413627187338129653105774703043377776905956484932486183722303201135571583745806165441941755833466966188398807387661571", 10).unwrap()).unwrap();
         let x = QuadraticForm::generator(&d).mul(&BigInt::from(1234));
         let y = QuadraticForm::generator(&d).mul(&BigInt::from(4321));
+        let z = y.clone();
         group.bench_function("Compose (512 bit discriminant)", move |b| {
             b.iter(|| x.compose(&y))
         });
+        group.bench_function("Double (512 bit discriminant)", move |b| {
+            b.iter(|| z.double())
+        });
 
         let d = Discriminant::try_from(BigInt::from_str_radix("-173197108158285529655099692042166386683260486655764503111574151459397279244340625070436917386670107433539464870917173822190635872887684166173874718269704667936351650895772937202272326332043347073303124000059154982400685660701006453457007094026343973435157790533480400962985543272080923974737725172126369794019", 10).unwrap()).unwrap();
         let x = QuadraticForm::generator(&d).mul(&BigInt::from(1234));
         let y = QuadraticForm::generator(&d).mul(&BigInt::from(4321));
+        let z = y.clone();
         group.bench_function("Compose (1024 bit discriminant)", move |b| {
             b.iter(|| x.compose(&y))
         });
+        group.bench_function("Double (1024 bit discriminant)", move |b| {
+            b.iter(|| z.double())
+        });
     }
 
     criterion_group! {

fastcrypto/src/groups/class_group/bigint_utils.rs

@@ -0,0 +1,111 @@
+// Copyright (c) 2022, Mysten Labs, Inc.
+// SPDX-License-Identifier: Apache-2.0
+
+use num_bigint::BigInt;
+use num_integer::Integer;
+use num_traits::{One, Signed, Zero};
+use std::mem;
+use std::ops::Neg;
+
+pub struct EuclideanAlgorithmOutput {
+    pub gcd: BigInt,
+    pub x: BigInt,
+    pub y: BigInt,
+    pub a_divided_by_gcd: BigInt,
+    pub b_divided_by_gcd: BigInt,
+}
+
+impl EuclideanAlgorithmOutput {
+    fn flip(self) -> Self {
+        Self {
+            gcd: self.gcd,
+            x: self.y,
+            y: self.x,
+            a_divided_by_gcd: self.b_divided_by_gcd,
+            b_divided_by_gcd: self.a_divided_by_gcd,
+        }
+    }
+}
+
+/// Compute the greatest common divisor gcd of a and b. The output also returns the Bezout coefficients
+/// x and y such that ax + by = gcd and also the quotients a / gcd and b / gcd.
+pub fn extended_euclidean_algorithm(a: &BigInt, b: &BigInt) -> EuclideanAlgorithmOutput {
+    if b < a {
+        return extended_euclidean_algorithm(b, a).flip();
+    }
+
+    let mut s = (BigInt::zero(), BigInt::one());
+    let mut t = (BigInt::one(), BigInt::zero());
+    let mut r = (a.clone(), b.clone());
+
+    while !r.0.is_zero() {
+        let (q, r_prime) = r.1.div_rem(&r.0);
+        r.1 = r.0;
+        r.0 = r_prime;
+
+        let f = |mut x: (BigInt, BigInt)| {
+            mem::swap(&mut x.0, &mut x.1);
+            x.0 -= &q * &x.1;
+            x
+        };
+        s = f(s);
+        t = f(t);
+    }
+
+    // The last coefficients are equal to +/- a / gcd(a,b) and b / gcd(a,b) respectively.
+    let a_divided_by_gcd = if a.sign() != s.0.sign() {
+        s.0.neg()
+    } else {
+        s.0
+    };
+    let b_divided_by_gcd = if b.sign() != t.0.sign() {
+        t.0.neg()
+    } else {
+        t.0
+    };
+
+    if !r.1.is_negative() {
+        EuclideanAlgorithmOutput {
+            gcd: r.1,
+            x: t.1,
+            y: s.1,
+            a_divided_by_gcd,
+            b_divided_by_gcd,
+        }
+    } else {
+        EuclideanAlgorithmOutput {
+            gcd: r.1.neg(),
+            x: t.1.neg(),
+            y: s.1.neg(),
+            a_divided_by_gcd,
+            b_divided_by_gcd,
+        }
+    }
+}
+
+#[test]
+fn test_xgcd() {
+    let a = BigInt::from(240);
+    let b = BigInt::from(46);
+    let output = extended_euclidean_algorithm(&a, &b);
+    assert_eq!(output.gcd, a.gcd(&b));
+    assert_eq!(&output.x * &a + &output.y * &b, output.gcd);
+    assert_eq!(output.a_divided_by_gcd, &a / &output.gcd);
+    assert_eq!(output.b_divided_by_gcd, &b / &output.gcd);
+
+    let a = BigInt::from(240);
+    let b = BigInt::from(-46);
+    let output = extended_euclidean_algorithm(&a, &b);
+    assert_eq!(output.gcd, a.gcd(&b));
+    assert_eq!(&output.x * &a + &output.y * &b, output.gcd);
+    assert_eq!(output.a_divided_by_gcd, &a / &output.gcd);
+    assert_eq!(output.b_divided_by_gcd, &b / &output.gcd);
+
+    let a = BigInt::from(-240);
+    let b = BigInt::from(-46);
+    let output = extended_euclidean_algorithm(&a, &b);
+    assert_eq!(output.gcd, a.gcd(&b));
+    assert_eq!(&output.x * &a + &output.y * &b, output.gcd);
+    assert_eq!(output.a_divided_by_gcd, &a / &output.gcd);
+    assert_eq!(output.b_divided_by_gcd, &b / &output.gcd);
+}