Add montgomery multiplication

curve/examples/schnorr.rs

@@ -43,5 +43,5 @@ fn main() {
     );
     assert_eq!(lhs, rhs, "points did not match");
     println!("Schnorr signature verified");
-    println!("Took {} milliseconds", time.elapsed().unwrap().as_millis());
+    println!("Took {} microseconds", time.elapsed().unwrap().as_micros());
 }

curve/src/lib.rs

@@ -40,39 +40,43 @@ pub trait Curve<F: Field>: fmt::Debug + Clone + Copy + PartialEq + Eq {
 pub struct Secp256k1Curve;
 
 impl Secp256k1Curve {
-    pub const A: FieldElement<Secp256k1FieldOrder> =
-        FieldElement::from_limbs_unchecked([0, 0, 0, 0]);
-    pub const B: FieldElement<Secp256k1FieldOrder> =
-        FieldElement::from_limbs_unchecked([7, 0, 0, 0]);
+    pub fn a() -> FieldElement<Secp256k1FieldOrder> {
+        FieldElement::from_limbs_unchecked([0, 0, 0, 0])
+    }
 
-    pub const GENERATOR: Point<Secp256k1FieldOrder> = Point::from_jacobian_unchecked(
-        FieldElement::from_limbs_unchecked([
-            0x59F2815B16F81798,
-            0x029BFCDB2DCE28D9,
-            0x55A06295CE870B07,
-            0x79BE667EF9DCBBAC,
-        ]),
-        FieldElement::from_limbs_unchecked([
-            0x9C47D08FFB10D4B8,
-            0xFD17B448A6855419,
-            0x5DA4FBFC0E1108A8,
-            0x483ADA7726A3C465,
-        ]),
-        FieldElement::from_limbs_unchecked([1, 0, 0, 0]),
-    );
+    pub fn b() -> FieldElement<Secp256k1FieldOrder> {
+        FieldElement::from_limbs_unchecked([7, 0, 0, 0])
+    }
+
+    pub fn generator() -> Point<Secp256k1FieldOrder> {
+        Point::from_affine(
+            FieldElement::from_limbs_unchecked([
+                0x59F2815B16F81798,
+                0x029BFCDB2DCE28D9,
+                0x55A06295CE870B07,
+                0x79BE667EF9DCBBAC,
+            ]),
+            FieldElement::from_limbs_unchecked([
+                0x9C47D08FFB10D4B8,
+                0xFD17B448A6855419,
+                0x5DA4FBFC0E1108A8,
+                0x483ADA7726A3C465,
+            ]),
+        )
+    }
 
     pub fn point_from_scalar(scalar: Scalar<Secp256k1GroupOrder>) -> Point<Secp256k1FieldOrder> {
-        Self.multiply(scalar, Self::GENERATOR)
+        Self.multiply(scalar, Self::generator())
     }
 }
 
 impl Curve<Secp256k1FieldOrder> for Secp256k1Curve {
     fn a(&self) -> FieldElement<Secp256k1FieldOrder> {
-        Self::A
+        Self::a()
     }
 
     fn b(&self) -> FieldElement<Secp256k1FieldOrder> {
-        Self::B
+        Self::b()
     }
 }
 
@@ -448,38 +452,38 @@ mod tests {
 
     #[test]
     fn struct_generator_matches_local_generator() {
-        assert_eq!(Secp256k1Curve::GENERATOR, generator());
+        assert_eq!(Secp256k1Curve::generator(), generator());
     }
 
     #[test]
     fn struct_generator_mul_by_one_is_itself() {
         assert_eq!(
-            CURVE.multiply(Sc::from_u128(1), Secp256k1Curve::GENERATOR),
-            Secp256k1Curve::GENERATOR
+            CURVE.multiply(Sc::from_u128(1), Secp256k1Curve::generator()),
+            Secp256k1Curve::generator()
         );
     }
 
     #[test]
     fn struct_generator_mul_by_two_matches_2g() {
         assert_eq!(
-            CURVE.multiply(Sc::from_u128(2), Secp256k1Curve::GENERATOR),
+            CURVE.multiply(Sc::from_u128(2), Secp256k1Curve::generator()),
             two_g()
         );
     }
 
     #[test]
     fn struct_generator_mul_by_three_matches_3g() {
         assert_eq!(
-            CURVE.multiply(Sc::from_u128(3), Secp256k1Curve::GENERATOR),
+            CURVE.multiply(Sc::from_u128(3), Secp256k1Curve::generator()),
             three_g()
         );
     }
 
     #[test]
     fn struct_generator_mul_scalar_additivity() {
-        let a = CURVE.multiply(Sc::from_u128(7), Secp256k1Curve::GENERATOR);
-        let b = CURVE.multiply(Sc::from_u128(11), Secp256k1Curve::GENERATOR);
-        let sum = CURVE.multiply(Sc::from_u128(18), Secp256k1Curve::GENERATOR);
+        let a = CURVE.multiply(Sc::from_u128(7), Secp256k1Curve::generator());
+        let b = CURVE.multiply(Sc::from_u128(11), Secp256k1Curve::generator());
+        let sum = CURVE.multiply(Sc::from_u128(18), Secp256k1Curve::generator());
         assert_eq!(CURVE.add(a, b), sum);
     }
 
@@ -490,7 +494,7 @@ mod tests {
             0xAE, 0xBA, 0xFE, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
             0xFF, 0xFF, 0xFF, 0xFF,
         ]);
-        assert!(CURVE.multiply(n, Secp256k1Curve::GENERATOR).is_infinity());
+        assert!(CURVE.multiply(n, Secp256k1Curve::generator()).is_infinity());
     }
 
     #[test]

field/scripts/mont_constants.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python3
+"""Print Montgomery constants for a prime modulus.
+
+Usage:
+    python3 mont_constants.py <modulus>
+    python3 mont_constants.py 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F
+
+Outputs, for R = 2^(64 * n) where n is the number of u64 limbs needed to hold p:
+  * R^2 mod p as little-endian u64 limbs
+  * -p^-1 mod 2^64 (the CIOS/SOS reduction multiplier)
+"""
+
+import sys
+
+
+def parse_modulus(s: str) -> int:
+    s = s.strip().replace("_", "")
+    if s.lower().startswith("0x"):
+        return int(s, 16)
+    return int(s)
+
+
+def limbs_le(x: int, n: int) -> list[int]:
+    mask = (1 << 64) - 1
+    return [(x >> (64 * i)) & mask for i in range(n)]
+
+
+def main() -> int:
+    if len(sys.argv) != 2:
+        print(f"usage: {sys.argv[0]} <modulus>", file=sys.stderr)
+        return 1
+    p = parse_modulus(sys.argv[1])
+    if p <= 1 or p % 2 == 0:
+        print("modulus must be an odd integer > 1", file=sys.stderr)
+        return 1
+
+    n = (p.bit_length() + 63) // 64
+    r = 1 << (64 * n)
+    r_squared = (r * r) % p
+    p_inv_pos = pow(p, -1, 1 << 64)
+    p_inv_neg = (-p_inv_pos) & ((1 << 64) - 1)
+
+    print(f"limbs (n): {n}")
+    print(f"R^2 mod p:")
+    for i, limb in enumerate(limbs_le(r_squared, n)):
+        print(f"  [{i}] 0x{limb:016X}")
+    print(f"-p^-1 mod 2^64: 0x{p_inv_neg:016X}")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())