Barrett reduction!

src/cryptonum/barrett.rs

@@ -1,12 +1,19 @@
 use cryptonum::{U192,   U256,   U384,   U512,   U576,
                 U1024,  U2048,  U3072,  U4096,  U8192,
                 U15360};
+use cryptonum::addition::raw_addition;
+use cryptonum::comparison::{bignum_cmp,bignum_ge};
 use cryptonum::division::divmod;
+use cryptonum::multiplication::raw_multiplication;
+use cryptonum::subtraction::raw_subtraction;
+use std::cmp::{Ordering,min};
 use std::fmt;
 
 macro_rules! generate_barrett_implementations {
     ($bname: ident, $name: ident, $size: expr) => {
         pub struct $bname {
+            pub(crate) k:  usize,
+            pub(crate) m:  [u64; $size/32],
             pub(crate) mu: [u64; $size/32]
         }
 
@@ -32,25 +39,86 @@ macro_rules! generate_barrett_implementations {
             }
         }
 
-
         impl $bname {
-            fn new(m: &$name) -> $bname {
+            pub fn new(m: &$name) -> $bname {
                 let mut b      = [0; ($size/32) + 1];
                 let mut widerm = [0; ($size/32) + 1];
                 let mut quot   = [0; ($size/32) + 1];
                 let mut remndr = [0; ($size/32) + 1];
-                let mut result = $bname{ mu: [0; $size/32] };
+                let mut result = $bname{ k: 0,
+                                         m: [0; $size/32],
+                                         mu: [0; $size/32] };
 
-                b[$size/32] = 1;
-                for (idx, val) in m.values.iter().enumerate() { widerm[idx] = *val; }
+                for (idx, val) in m.values.iter().enumerate() {
+                    let x = *val;
+                    widerm[idx] = x;
+                    result.m[idx] = x;
+                    if x != 0 { result.k = idx; }
+                }
+                result.k += 1;
+                b[result.k*2] = 1;
                 divmod(&b, &widerm, &mut quot, &mut remndr);
                 for (idx, val) in result.mu.iter_mut().enumerate() { *val = quot[idx]; }
                 result
             }
+
+            pub fn reduce(&self, x: &mut $name) {
+                printvar("x", &x.values);
+                printvar("m", &self.m);
+                printvar("u", &self.mu);
+                // 1. q1←⌊x/bk−1⌋, q2←q1 · μ, q3←⌊q2/bk+1⌋.
+                let mut q1 = [0; $size/32];
+                shiftr(&x.values, self.k - 1, &mut q1);
+                let mut q2 = [0; $size/16];
+                raw_multiplication(&q1, &self.mu, &mut q2);
+                let mut q3 = [0; $size/16];
+                shiftr(&q2, self.k + 1,&mut q3);
+                // 2. r1←x mod bk+1, r2←q3 · m mod bk+1, r←r1 − r2.
+                let mut r = [0; $size/16];
+                let copylen = min(self.k+1, x.values.len());
+                for i in 0..copylen { r[i] = x.values[i]; }
+                let mut r2big = [0; $size/8];
+                let mut mwider = [0; $size/16];
+                for i in 0..$size/32 { mwider[i] = self.m[i]; }
+                raw_multiplication(&q3, &mwider, &mut r2big);
+                let mut r2 = [0; $size/16];
+                for i in 0..self.k+1 { r2[i] = r2big[i]; }
+                let went_negative = !bignum_ge(&r, &r2);
+                raw_subtraction(&mut r, &r2);
+                // 3. If r<0 then r←r+bk+1.
+                if went_negative {
+                    let mut bk1 = [0; $size/32];
+                    bk1[self.k+1] = 1;
+                    raw_addition(&mut r, &bk1);
+                }
+                // 4. While r≥m do: r←r−m.
+                while bignum_cmp(&r, &mwider) == Ordering::Greater {
+                    raw_subtraction(&mut r, &mwider);
+                }
+                // Copy it over.
+                for (idx, val) in x.values.iter_mut().enumerate() {
+                    *val = r[idx];
+                }
+            }
         }
     };
 }
 
+fn printvar(name: &'static str, val: &[u64]) {
+    print!("{}: 0x", name);
+    for x in val.iter().rev() {
+        print!("{:016X}", *x);
+    }
+    println!("");
+}
+
+fn shiftr(x: &[u64], amt: usize, dest: &mut [u64])
+{
+    for i in amt..x.len() {
+        dest[i-amt] = x[i];
+    }
+}
+
 
 generate_barrett_implementations!(BarrettU192,   U192,  192);
 generate_barrett_implementations!(BarrettU256,   U256,  256);
@@ -78,13 +146,19 @@ macro_rules! generate_tests {
                 fn $name() {
                     let fname = format!("tests/math/barrett_gen{}.test",
                                         stringify!($name));
-                    run_test(fname.to_string(), 2, |case| {
+                    run_test(fname.to_string(), 3, |case| {
                         let (neg0, mbytes) = case.get("m").unwrap();
                         let (neg1, ubytes) = case.get("u").unwrap();
+                        let (neg2, kbytes) = case.get("k").unwrap();
 
-                        assert!(!neg0 && !neg1);
+                        assert!(!neg0 && !neg1 && !neg2);
                         let m = $name::from_bytes(mbytes);
-                        let mut u = $bname{ mu: [0; $size/32]};
+                        let mut kbig = [0; 1];
+                        raw_decoder(&kbytes, &mut kbig);
+                        let mut u = $bname{ k: kbig[0] as usize,
+                                            m: [0; $size/32],
+                                            mu: [0; $size/32]};
+                        raw_decoder(&mbytes, &mut u.m);
                         raw_decoder(&ubytes, &mut u.mu);
                         let r = $bname::new(&m);
                         assert_eq!(u,r);
@@ -92,6 +166,42 @@ macro_rules! generate_tests {
                 }
             )*
         }
+
+        #[cfg(test)]
+        mod reduction {
+            use cryptonum::encoding::{Decoder,raw_decoder};
+            use super::*;
+            use testing::run_test;
+
+            $(
+                #[test]
+                #[allow(non_snake_case)]
+                fn $name() {
+                    let fname = format!("tests/math/barrett_reduce{}.test",
+                                        stringify!($name));
+                    run_test(fname.to_string(), 5, |case| {
+                        let (neg0, mbytes) = case.get("m").unwrap();
+                        let (neg1, ubytes) = case.get("u").unwrap();
+                        let (neg2, kbytes) = case.get("k").unwrap();
+                        let (neg3, xbytes) = case.get("x").unwrap();
+                        let (neg4, rbytes) = case.get("r").unwrap();
+
+                        assert!(!neg0 && !neg1 && !neg2 && !neg3 && !neg4);
+                        let mut kbig = [0; 1];
+                        raw_decoder(&kbytes, &mut kbig);
+                        let mut u = $bname{ k: kbig[0] as usize,
+                                            m: [0; $size/32],
+                                            mu: [0; $size/32]};
+                        raw_decoder(&mbytes, &mut u.m);
+                        raw_decoder(&ubytes, &mut u.mu);
+                        let mut x = $name::from_bytes(&xbytes);
+                        let r = $name::from_bytes(&rbytes);
+                        u.reduce(&mut x);
+                        assert_eq!(x, r);
+                    });
+                }
+            )*
+        }
     }
 }
 

tests/math/GenerateTests.hs

@@ -2,6 +2,7 @@ import Control.Monad
 import Data.Bits(shiftL,(.&.))
 import Data.Map.Strict(Map)
 import qualified Data.Map.Strict as Map
+import GHC.Integer.GMP.Internals(powModInteger)
 import Numeric(showHex)
 import Prelude hiding (log)
 import System.IO(hFlush,stdout,IOMode(..),withFile,Handle,hClose,hPutStrLn)
@@ -16,7 +17,10 @@ testTypes = [("addition", addTest),
              ("modmul", modmulTest),
              ("squaring", squareTest),
              ("modsq", modsqTest),
-             ("division", divTest)
+             ("modexp", modexpTest),
+             ("division", divTest),
+             ("barrett_gen", barrettGenTest),
+             ("barrett_reduce", barrettReduceTest)
             ]
 
 bitSizes :: [Int]
@@ -131,6 +135,19 @@ modsqTest bitsize gen0 = (res, gen1)
                             ("m", showHex m' ""),
                             ("r", showHex r  "")]
 
+modexpTest :: Int -> StdGen -> (Map String String, StdGen)
+modexpTest bitsize gen0 = (res, gen2)
+ where
+  (b, gen1)  = random gen0
+  (e, gen2)  = random gen1
+  (m, gen3)  = random gen2
+  [b',e',m'] = splitMod bitsize [b,e,m]
+  r          = powModInteger b' e' m'
+  res        = Map.fromList [("b", showHex b' ""),
+                             ("e", showHex e' ""),
+                             ("m", showHex m' ""),
+                             ("r", showHex r  "")]
+
 
 divTest :: Int -> StdGen -> (Map String String, StdGen)
 divTest bitsize gen0 = (res, gen2)
@@ -145,6 +162,45 @@ divTest bitsize gen0 = (res, gen2)
                             ("q", showHex q  ""),
                             ("r", showHex r  "")]
 
+barrettGenTest :: Int -> StdGen -> (Map String String, StdGen)
+barrettGenTest bitsize gen0 = (res, gen1)
+ where
+  (m, gen1) = random gen0
+  m'        = m .&. mask bitsize
+  k         = computeK m'
+  u         = barrett bitsize m'
+  res       = Map.fromList [("m", showHex m' ""),
+                            ("k", showHex k  ""),
+                            ("u", showHex u  "")]
+
+barrettReduceTest :: Int -> StdGen -> (Map String String, StdGen)
+barrettReduceTest bitsize gen0 = (res, gen2)
+ where
+  (m, gen1) = random gen0
+  (x, gen2) = random gen1
+  m'        = m .&. mask bitsize
+  x'        = x .&. mask (min bitsize (2 * k * 64))
+  k         = computeK m'
+  u         = barrett bitsize m'
+  r         = x' `mod` m'
+  res       = Map.fromList [("m", showHex m' ""),
+                            ("x", showHex x' ""),
+                            ("k", showHex k  ""),
+                            ("u", showHex u  ""),
+                            ("r", showHex r  "")]
+
+barrett :: Int -> Integer -> Integer
+barrett bitsize m = (b ^ (2 * k)) `div` m
+ where
+  b = 2 ^ 64
+  k = computeK m
+
+computeK :: Integer -> Int
+computeK v = go 0 1
+ where
+  go k acc | v < acc   = k + 1
+           | otherwise = go (k + 1) (acc * (2 ^ 64))
+
 log :: String -> IO ()
 log str =
   do putStr str