implement REDC (no serious testing yet)

README.md

@@ -17,4 +17,4 @@ and [`staging-agda`](https://github.com/faulhornlabs/staging-agda/).
 - [ ] add a 64 bit version
 - [ ] add more Poseidon2 state widths than `t=3`
 - [ ] implement `circomlib`-compatible Poseidon
-
+- [ ] further optimizations

sanity.hs

@@ -22,6 +22,10 @@ big1   = 0x3142edd042d60bd1c80322cb76b2f85cdc3f0f00b0151d71e4db8a57191eb5c7
 big2   = 0x7d7cf5c0086f97e30ac1f5a6987ee9529c417d290b0c5fe203e4f3728cbc1b97
 big3   = 0xbbfc0dc5195a120a2317351be185427ea931be03406c5062f9fff49bb1efa2d1
 
+felt1  = 0x22ac1ee66024036b5ab6f194bf51bf7fd03cc3b9ca5c5b8a00d4796720dc4a9f
+felt2  = 0x2063f06a7b59b4fa596d982aa362c094703bcf35eb5f4c171f9c48a6e34d39b8
+felt3  = 0x07d89e99d3bbaebca574aefcf8f492c5188cc945728b9d4da1c33f8831bfd5d3
+
 mont1  = 0x095fad8ecdabd8686a74e6c0c03384162216c484c62c6b29f7cb6585c7a4b4fc
 mont2  = 0x01a6741e42e6d6bc5f830d505947fe1ef1226480a55c65600b62d2a19c974559
 mont3  = 0x2e493a4b8bf5eb71facdfcdc73d11ffa9c15ed08671492cefbb9d8d024a2de63

src/bigint.rs

@@ -27,7 +27,7 @@ pub type BigInt512 = BigInt<16>;
 
 //------------------------------------------------------------------------------
 
-fn boolToU32(c: bool) -> u32 {
+pub fn boolToU32(c: bool) -> u32 {
   if c { 1 } else { 0 }
 }
 
@@ -57,6 +57,13 @@ impl<const N: usize> BigInt<N> {
     BigInt { limbs: ls }
   }
 
+  pub fn truncate1(big : &BigInt<{N+1}>) -> BigInt<N> {
+    // let small: [u32; N] = &big.limbs[0..N];
+    let mut small: [u32; N] = [0; N];
+    for i in 0..N { small[i] = big.limbs[i]; }
+    BigInt { limbs: small }
+  }
+
   pub fn zero() -> BigInt<N> {
     BigInt { limbs: [0; N] }
   }

src/bin/testmain.rs

@@ -43,6 +43,7 @@ fn main() {
   Mont::print_internal("MONT1",&MONT1);
   Mont::print_internal("MONT2",&MONT2);
   Mont::print_internal("MONT3",&MONT3);
+  println!("-----");
   Mont::print_standard("MONT1",&MONT1);
   Mont::print_standard("MONT2",&MONT2);
   Mont::print_standard("MONT3",&MONT3);
@@ -92,8 +93,8 @@ fn main() {
   println!("-----");
 
   println!("F1*F2 = {}", Felt::mul( &FELT1, &FELT2) );  
-  println!("F1*F2 = {}", Felt::mul( &FELT1, &FELT2) );  
-  println!("F1*F2 = {}", Felt::mul( &FELT1, &FELT2) );  
+  println!("F2*F3 = {}", Felt::mul( &FELT2, &FELT3) );  
+  println!("F3*F1 = {}", Felt::mul( &FELT3, &FELT1) );  
 
 }
 

src/constant.rs

@@ -6,9 +6,10 @@ use crate::bigint::*;
 
 type Big = BigInt<8>;
 
-pub const FIELD_PRIME  : Big = BigInt::make( [ 0xf0000001 , 0x43e1f593 , 0x79b97091 , 0x2833e848 , 0x8181585d , 0xb85045b6 , 0xe131a029 , 0x30644e72 ] );
-pub const PRIME_PLUS_1 : Big = BigInt::make( [ 0xf0000002 , 0x43e1f593 , 0x79b97091 , 0x2833e848 , 0x8181585d , 0xb85045b6 , 0xe131a029 , 0x30644e72 ] );
-pub const HALFP_PLUS_1 : Big = BigInt::make( [ 0xf8000001 , 0xa1f0fac9 , 0x3cdcb848 , 0x9419f424 , 0x40c0ac2e , 0xdc2822db , 0x97098d01 , 0x01832273 ] ); 
+pub const PRIME_EXT    : BigInt<9> = BigInt::make( [ 0xf0000001 , 0x43e1f593 , 0x79b97091 , 0x2833e848 , 0x8181585d , 0xb85045b6 , 0xe131a029 , 0x30644e72 , 0x00000000 ] );
+pub const FIELD_PRIME  : Big       = BigInt::make( [ 0xf0000001 , 0x43e1f593 , 0x79b97091 , 0x2833e848 , 0x8181585d , 0xb85045b6 , 0xe131a029 , 0x30644e72 ] );
+pub const PRIME_PLUS_1 : Big       = BigInt::make( [ 0xf0000002 , 0x43e1f593 , 0x79b97091 , 0x2833e848 , 0x8181585d , 0xb85045b6 , 0xe131a029 , 0x30644e72 ] );
+pub const HALFP_PLUS_1 : Big       = BigInt::make( [ 0xf8000001 , 0xa1f0fac9 , 0x3cdcb848 , 0x9419f424 , 0x40c0ac2e , 0xdc2822db , 0x97098d01 , 0x01832273 ] ); 
 
 //------------------------------------------------------------------------------
 // montgomery constants

src/montgomery.rs

@@ -8,6 +8,7 @@
 
 use std::fmt;
 
+use crate::platform::*;
 use crate::bigint::*;
 use crate::constant::*;
 
@@ -85,11 +86,41 @@ impl Mont {
   }
 
   // the Montgomery reduction algorithm
+  // <https://en.wikipedia.org/wiki/Montgomery_modular_multiplication#Montgomery_arithmetic_on_multiprecision_integers>
   fn redc(input: BigInt<16>) -> Big {
-    let mut T: [u32; 17] = [0; 17];
-    for i in 0..16 { T[i] = input.limbs[i] }
 
-    BigInt::zero()
+    let mut T: [u32; 17] = [0; 17];
+    for i in 0..16 { T[i] = input.limbs[i]; }
+
+    for i in 0..8 {
+      let mut carry: u32 = 0;
+      let m: u32 = truncMul32( T[i] , MONT_Q );
+      for j in 0..8 {
+        let (lo,hi) = mulAddAdd32( m, FIELD_PRIME.limbs[j], carry, T[i+j] );
+        T[i+j] = lo;
+        carry  = hi;
+      }
+      for j in 8..(17-i) {
+        let (x,c) = addCarry32_( T[i+j] , carry );
+        T[i+j] = x;
+        carry  = boolToU32(c);
+      }
+    }
+
+    let mut S : [u32; 9] = [0; 9];
+    for i in 0..9 { S[i] = T[8+i]; }
+
+    let A     :  BigInt<9>       = BigInt { limbs: S };
+    let (B,c) : (BigInt<9>,bool) = BigInt::subBorrow( &A , &PRIME_EXT );
+
+    if c {
+      // `A - prime < 0` is equivalent to `A < prime` 
+      BigInt::truncate1(&A)
+    }
+    else {
+      // `A - prime >= 0` is equivalent to `A >= prime`
+      BigInt::truncate1(&B)
+    }
   }
 
   pub fn sqr(mont: &Mont) -> Mont {

src/platform.rs

@@ -21,8 +21,13 @@ pub fn subBorrow32(x: u32, y: u32, cin: bool) -> (u32,bool) {
   u32::borrowing_sub(x,y,cin)
 }
 
+pub fn truncMul32(x: u32, y: u32) -> u32 {
+  u32::wrapping_mul(x,y)
+}
+
 pub fn extMul32(x: u32, y: u32) -> (u32,u32) {
-  u32::carrying_mul(x,y,0)
+  // u32::carrying_mul(x,y,0)
+  u32::widening_mul(x,y)
 }
 
 pub fn mulAdd32(x: u32, y: u32, a: u32) -> (u32,u32) {