Basic mathematics for unsigned integers ported over, tests pass.

generate.hs

@@ -2,6 +2,7 @@ import Control.Monad(foldM_,forM_,when)
 import Data.Bits(shiftL,shiftR)
 import Data.List(sort)
 import qualified Data.Map.Strict as Map
+import GHC.Integer.GMP.Internals(powModInteger)
 import Numeric(showHex)
 import Prelude hiding (log)
 import System.Directory(createDirectoryIfMissing)
@@ -19,6 +20,7 @@ data Operation = Add
                | ModSq
                | Mul
                | Shifts
+               | Square
                | Sub
                | Convert Int
  deriving (Eq, Ord, Show)
@@ -33,7 +35,8 @@ needs = [ Need ModExp   (\ size -> [Req size ModMul
                                    ,Req size ModSq
                                    ,Req size Barretts])
         , Need ModSq    (\ size -> [Req (size * 2) Div
-                                   ,Req size Barretts])
+                                   ,Req size Barretts
+                                   ,Req size Square])
         , Need ModMul   (\ size -> [Req size Mul
                                    ,Req size Barretts
                                    ,Req size (Convert (size * 2))
@@ -127,9 +130,10 @@ generateInvocs =
            Div        -> hPutStrLn hndl ("div_impls!(U" ++ show size ++ ", U" ++ show (size * 2) ++ ");")
            ModExp     -> hPutStrLn hndl ("modexp_impls!(U" ++ show size ++ ");")
            ModMul     -> hPutStrLn hndl ("modmul_impls!(U" ++ show size ++ ", U" ++ show (size * 2) ++ ");")
-           ModSq      -> hPutStrLn hndl ("modsq_impls!(U" ++ show size ++ ");")
+           ModSq      -> hPutStrLn hndl ("modsq_impls!(U" ++ show size ++ ", U" ++ show (size * 2) ++ ");")
            Mul        -> hPutStrLn hndl ("multiply_impls!(U" ++ show size ++ ", U" ++ show (size * 2) ++ ");")
            Shifts     -> hPutStrLn hndl ("shift_impls!(U" ++ show size ++ ", " ++ show (size `div` 64) ++ ");")
+           Square     -> hPutStrLn hndl ("square_impls!(U" ++ show size ++ ", U" ++ show (size * 2) ++ ", " ++ show size ++ ");")
            Sub        -> hPutStrLn hndl ("subtraction_impls!(U" ++ show size ++ ", " ++ show (size `div` 64) ++ ");")
            Convert to -> hPutStrLn hndl ("conversion_impls!(U" ++ show size ++ ", U" ++ show to ++ ");")
        hPutStrLn hndl ""
@@ -147,6 +151,10 @@ generateInvocs =
        generateTestBlock hndl "div"         Div      True  []
        generateTestBlock hndl "barrett_gen" Barretts True  [(+ 64)]
        generateTestBlock hndl "barrett_red" Barretts True  [(+ 64), (* 2)]
+       generateTestBlock hndl "modsq"       ModSq    True  []
+       generateTestBlock hndl "modmul"      ModMul   True  []
+       generateTestBlock hndl "modexp"      ModExp   True  []
+       generateTestBlock hndl "square"      Square   True  [(* 2)]
        hPutStrLn hndl "}"
 
 log :: String -> IO ()
@@ -243,7 +251,7 @@ generateAllTheTests =
            res          = Map.fromList [("m", showX m), ("k", showX k),
                                         ("u", showX u)]
        in (res, u, memory1)
-     let (db3, gen3) = emptyDatabase gen1
+     let (db3, gen3) = emptyDatabase gen2
      generateTests Barretts "barrett_reduce" db3 $ \ size memory0 ->
        let (m, memory1) = generateNum memory0 "m" size
            (x, memory2) = generateNum memory1 "x" (min size (2 * k * 64))
@@ -254,7 +262,7 @@ generateAllTheTests =
                                         ("k", showX k), ("u", showX u),
                                         ("r", showX r)]
        in (res, r, memory2)
-     let (db4, gen4) = emptyDatabase gen2
+     let (db4, gen4) = emptyDatabase gen3
      generateTests BaseOps "base" db4 $ \ size memory0 ->
        let (x, memory1) = generateNum memory0 "x" size
            (m, memory2) = generateNum memory1 "m" size
@@ -264,7 +272,7 @@ generateAllTheTests =
                                         ("e", showB (even x)), ("o", showB (odd x)),
                                         ("m", showX m'),       ("r", showX r)]
        in (res, x, memory2)
-     let (db5, gen5) = emptyDatabase gen3
+     let (db5, gen5) = emptyDatabase gen4
      generateTests BaseOps "cmp" db5 $ \ size memory0 ->
        let (a, memory1) = generateNum memory0 "a" size
            (b, memory2) = generateNum memory1 "b" size
@@ -272,7 +280,7 @@ generateAllTheTests =
                                         ("g", showB (a > b)), ("l", showB (a < b)),
                                         ("e", showB (a == b))]
        in (res, a, memory2)
-     let (db6, gen6) = emptyDatabase gen4
+     let (db6, gen6) = emptyDatabase gen5
      generateTests Div "div" db6 $ \ size memory0 ->
        let (a, memory1) = generateNum memory0 "a" size
            (b, memory2) = generateNum memory1 "b" size
@@ -281,7 +289,7 @@ generateAllTheTests =
            res          = Map.fromList [("a", showX a), ("b", showX b),
                                         ("q", showX q), ("r", showX r)]
        in (res, q, memory2)
-     let (db7, gen7) = emptyDatabase gen5
+     let (db7, gen7) = emptyDatabase gen6
      generateTests Mul "mul" db7 $ \ size memory0 ->
        let (a, memory1) = generateNum memory0 "a" size
            (b, memory2) = generateNum memory1 "b" size
@@ -289,7 +297,7 @@ generateAllTheTests =
            res          = Map.fromList [("a", showX a), ("b", showX b),
                                         ("c", showX c)]
        in (res, c, memory2)
-     let (db8, gen8) = emptyDatabase gen6
+     let (db8, gen8) = emptyDatabase gen7
      generateTests Shifts "shiftl" db8 $ \ size memory0 ->
        let (a, memory1) = generateNum memory0 "a" size
            (l, memory2) = generateNum memory1 "l" size
@@ -297,7 +305,7 @@ generateAllTheTests =
            r            = modulate (a `shiftL` fromIntegral l') size
            res          = Map.fromList [("a", showX a), ("l", showX l'), ("r", showX r)]
        in (res, r, memory2)
-     let (db9, gen9) = emptyDatabase gen6
+     let (db9, gen9) = emptyDatabase gen8
      generateTests Shifts "shiftr" db9 $ \ size memory0 ->
        let (a, memory1) = generateNum memory0 "a" size
            (l, memory2) = generateNum memory1 "l" size
@@ -305,7 +313,7 @@ generateAllTheTests =
            r            = modulate (a `shiftR` fromIntegral l') size
            res          = Map.fromList [("a", showX a), ("l", showX l'), ("r", showX r)]
        in (res, l, memory2)
-     let (dbA, genA) = emptyDatabase gen7
+     let (dbA, genA) = emptyDatabase gen9
      generateTests Sub "sub" dbA $ \ size memory0 ->
        let (a, memory1) = generateNum memory0 "a" size
            (b, memory2) = generateNum memory1 "b" size
@@ -313,6 +321,38 @@ generateAllTheTests =
            res          = Map.fromList [("a", showX a), ("b", showX b),
                                         ("c", showX c)]
        in (res, c, memory2)
+     let (dbB, genB) = emptyDatabase genA
+     generateTests ModSq "modsq" dbB $ \ size memory0 ->
+       let (a, memory1) = generateNum memory0 "a" size
+           (m, memory2) = generateNum memory1 "m" size
+           c            = (a * a) `mod` m
+           res          = Map.fromList [("a", showX a), ("m", showX m),
+                                        ("c", showX c)]
+       in (res, c, memory2)
+     let (dbC, genC) = emptyDatabase genB
+     generateTests ModMul "modmul" dbC $ \ size memory0 ->
+       let (a, memory1) = generateNum memory0 "a" size
+           (b, memory2) = generateNum memory1 "b" size
+           (m, memory3) = generateNum memory2 "m" size
+           c            = (a * b) `mod` m
+           res          = Map.fromList [("a", showX a), ("b", showX b),
+                                        ("m", showX m), ("c", showX c)]
+       in (res, c, memory3)
+     let (dbD, genD) = emptyDatabase genC
+     generateTests ModExp "modexp" dbD $ \ size memory0 ->
+       let (b, memory1) = generateNum memory0 "b" size
+           (e, memory2) = generateNum memory1 "e" size
+           (m, memory3) = generateNum memory2 "m" size
+           r            = powModInteger b e m
+           res          = Map.fromList [("b", showX b), ("e", showX e),
+                                        ("m", showX m), ("r", showX r)]
+       in (res, r, memory3)
+     let (dbE, genE) = emptyDatabase genC
+     generateTests Square "square" dbE $ \ size memory0 ->
+       let (a, memory1) = generateNum memory0 "a" size
+           r            = modulate    (a * a)     (2 * size)
+           res          = Map.fromList [("a", showX a), ("r", showX r)]
+       in (res, r, memory1)
 
 main :: IO ()
 main =

src/unsigned/invoc.rs

@@ -2,8 +2,9 @@ base_impls!(U192, 3);
 barrett_impl!(BarrettU192, U192, U256, U384, U448);
 modexp_impls!(U192);
 modmul_impls!(U192, U384);
-modsq_impls!(U192);
+modsq_impls!(U192, U384);
 multiply_impls!(U192, U384);
+square_impls!(U192, U384, 192);
 conversion_impls!(U192, U256);
 conversion_impls!(U192, U384);
 conversion_impls!(U192, U448);
@@ -11,8 +12,9 @@ base_impls!(U256, 4);
 barrett_impl!(BarrettU256, U256, U320, U512, U576);
 modexp_impls!(U256);
 modmul_impls!(U256, U512);
-modsq_impls!(U256);
+modsq_impls!(U256, U512);
 multiply_impls!(U256, U512);
+square_impls!(U256, U512, 256);
 conversion_impls!(U256, U320);
 conversion_impls!(U256, U384);
 conversion_impls!(U256, U448);
@@ -27,9 +29,10 @@ barrett_impl!(BarrettU384, U384, U448, U768, U832);
 div_impls!(U384, U768);
 modexp_impls!(U384);
 modmul_impls!(U384, U768);
-modsq_impls!(U384);
+modsq_impls!(U384, U768);
 multiply_impls!(U384, U768);
 shift_impls!(U384, 6);
+square_impls!(U384, U768, 384);
 conversion_impls!(U384, U448);
 conversion_impls!(U384, U768);
 conversion_impls!(U384, U832);
@@ -48,9 +51,10 @@ barrett_impl!(BarrettU512, U512, U576, U1024, U1088);
 div_impls!(U512, U1024);
 modexp_impls!(U512);
 modmul_impls!(U512, U1024);
-modsq_impls!(U512);
+modsq_impls!(U512, U1024);
 multiply_impls!(U512, U1024);
 shift_impls!(U512, 8);
+square_impls!(U512, U1024, 512);
 conversion_impls!(U512, U576);
 conversion_impls!(U512, U1024);
 conversion_impls!(U512, U1088);
@@ -60,9 +64,10 @@ barrett_impl!(BarrettU576, U576, U640, U1152, U1216);
 div_impls!(U576, U1152);
 modexp_impls!(U576);
 modmul_impls!(U576, U1152);
-modsq_impls!(U576);
+modsq_impls!(U576, U1152);
 multiply_impls!(U576, U1152);
 shift_impls!(U576, 9);
+square_impls!(U576, U1152, 576);
 subtraction_impls!(U576, 9);
 conversion_impls!(U576, U640);
 conversion_impls!(U576, U1024);
@@ -100,9 +105,10 @@ barrett_impl!(BarrettU1024, U1024, U1088, U2048, U2112);
 div_impls!(U1024, U2048);
 modexp_impls!(U1024);
 modmul_impls!(U1024, U2048);
-modsq_impls!(U1024);
+modsq_impls!(U1024, U2048);
 multiply_impls!(U1024, U2048);
 shift_impls!(U1024, 16);
+square_impls!(U1024, U2048, 1024);
 subtraction_impls!(U1024, 16);
 conversion_impls!(U1024, U1088);
 conversion_impls!(U1024, U2048);
@@ -150,9 +156,10 @@ barrett_impl!(BarrettU2048, U2048, U2112, U4096, U4160);
 div_impls!(U2048, U4096);
 modexp_impls!(U2048);
 modmul_impls!(U2048, U4096);
-modsq_impls!(U2048);
+modsq_impls!(U2048, U4096);
 multiply_impls!(U2048, U4096);
 shift_impls!(U2048, 32);
+square_impls!(U2048, U4096, 2048);
 subtraction_impls!(U2048, 32);
 conversion_impls!(U2048, U2112);
 conversion_impls!(U2048, U4096);
@@ -186,8 +193,9 @@ base_impls!(U3072, 48);
 barrett_impl!(BarrettU3072, U3072, U3136, U6144, U6208);
 modexp_impls!(U3072);
 modmul_impls!(U3072, U6144);
-modsq_impls!(U3072);
+modsq_impls!(U3072, U6144);
 multiply_impls!(U3072, U6144);
+square_impls!(U3072, U6144, 3072);
 conversion_impls!(U3072, U3136);
 conversion_impls!(U3072, U6144);
 conversion_impls!(U3072, U6208);
@@ -201,9 +209,10 @@ barrett_impl!(BarrettU4096, U4096, U4160, U8192, U8256);
 div_impls!(U4096, U8192);
 modexp_impls!(U4096);
 modmul_impls!(U4096, U8192);
-modsq_impls!(U4096);
+modsq_impls!(U4096, U8192);
 multiply_impls!(U4096, U8192);
 shift_impls!(U4096, 64);
+square_impls!(U4096, U8192, 4096);
 subtraction_impls!(U4096, 64);
 conversion_impls!(U4096, U4160);
 conversion_impls!(U4096, U8192);
@@ -242,8 +251,9 @@ base_impls!(U7680, 120);
 barrett_impl!(BarrettU7680, U7680, U7744, U15360, U15424);
 modexp_impls!(U7680);
 modmul_impls!(U7680, U15360);
-modsq_impls!(U7680);
+modsq_impls!(U7680, U15360);
 multiply_impls!(U7680, U15360);
+square_impls!(U7680, U15360, 7680);
 conversion_impls!(U7680, U7744);
 conversion_impls!(U7680, U15360);
 conversion_impls!(U7680, U15424);
@@ -257,9 +267,10 @@ barrett_impl!(BarrettU8192, U8192, U8256, U16384, U16448);
 div_impls!(U8192, U16384);
 modexp_impls!(U8192);
 modmul_impls!(U8192, U16384);
-modsq_impls!(U8192);
+modsq_impls!(U8192, U16384);
 multiply_impls!(U8192, U16384);
 shift_impls!(U8192, 128);
+square_impls!(U8192, U16384, 8192);
 subtraction_impls!(U8192, 128);
 conversion_impls!(U8192, U8256);
 conversion_impls!(U8192, U16384);
@@ -294,9 +305,10 @@ barrett_impl!(BarrettU15360, U15360, U15424, U30720, U30784);
 div_impls!(U15360, U30720);
 modexp_impls!(U15360);
 modmul_impls!(U15360, U30720);
-modsq_impls!(U15360);
+modsq_impls!(U15360, U30720);
 multiply_impls!(U15360, U30720);
 shift_impls!(U15360, 240);
+square_impls!(U15360, U30720, 15360);
 conversion_impls!(U15360, U15424);
 conversion_impls!(U15360, U30720);
 conversion_impls!(U15360, U30784);
@@ -899,4 +911,72 @@ mod tests {
     generate_barrett_red_tests!(U8192, u8192, U8256, U16384);
     generate_barrett_red_tests!(U15360, u15360, U15424, U30720);
   }
+  mod modsq {
+    use super::super::*;
+    use testing::run_test;
+
+    generate_modsq_tests!(U192, u192);
+    generate_modsq_tests!(U256, u256);
+    generate_modsq_tests!(U384, u384);
+    generate_modsq_tests!(U512, u512);
+    generate_modsq_tests!(U576, u576);
+    generate_modsq_tests!(U1024, u1024);
+    generate_modsq_tests!(U2048, u2048);
+    generate_modsq_tests!(U3072, u3072);
+    generate_modsq_tests!(U4096, u4096);
+    generate_modsq_tests!(U7680, u7680);
+    generate_modsq_tests!(U8192, u8192);
+    generate_modsq_tests!(U15360, u15360);
+  }
+  mod modmul {
+    use super::super::*;
+    use testing::run_test;
+
+    generate_modmul_tests!(U192, u192);
+    generate_modmul_tests!(U256, u256);
+    generate_modmul_tests!(U384, u384);
+    generate_modmul_tests!(U512, u512);
+    generate_modmul_tests!(U576, u576);
+    generate_modmul_tests!(U1024, u1024);
+    generate_modmul_tests!(U2048, u2048);
+    generate_modmul_tests!(U3072, u3072);
+    generate_modmul_tests!(U4096, u4096);
+    generate_modmul_tests!(U7680, u7680);
+    generate_modmul_tests!(U8192, u8192);
+    generate_modmul_tests!(U15360, u15360);
+  }
+  mod modexp {
+    use super::super::*;
+    use testing::run_test;
+
+    generate_modexp_tests!(U192, u192);
+    generate_modexp_tests!(U256, u256);
+    generate_modexp_tests!(U384, u384);
+    generate_modexp_tests!(U512, u512);
+    generate_modexp_tests!(U576, u576);
+    generate_modexp_tests!(U1024, u1024);
+    generate_modexp_tests!(U2048, u2048);
+    generate_modexp_tests!(U3072, u3072);
+    generate_modexp_tests!(U4096, u4096);
+    generate_modexp_tests!(U7680, u7680);
+    generate_modexp_tests!(U8192, u8192);
+    generate_modexp_tests!(U15360, u15360);
+  }
+  mod square {
+    use super::super::*;
+    use testing::run_test;
+
+    generate_square_tests!(U192, u192, U384);
+    generate_square_tests!(U256, u256, U512);
+    generate_square_tests!(U384, u384, U768);
+    generate_square_tests!(U512, u512, U1024);
+    generate_square_tests!(U576, u576, U1152);
+    generate_square_tests!(U1024, u1024, U2048);
+    generate_square_tests!(U2048, u2048, U4096);
+    generate_square_tests!(U3072, u3072, U6144);
+    generate_square_tests!(U4096, u4096, U8192);
+    generate_square_tests!(U7680, u7680, U15360);
+    generate_square_tests!(U8192, u8192, U16384);
+    generate_square_tests!(U15360, u15360, U30720);
+  }
 }

src/unsigned/mod.rs

@@ -19,10 +19,14 @@ mod modexp;
 #[macro_use]
 mod modmul;
 #[macro_use]
+mod modsq;
+#[macro_use]
 mod mul;
 #[macro_use]
 mod shifts;
 #[macro_use]
+mod square;
+#[macro_use]
 mod sub;
 
 use self::add::{addition,unsafe_addition};
@@ -31,9 +35,12 @@ use self::cmp::compare;
 use self::codec::{Encoder,Decoder,raw_decoder};
 use self::div::{DivMod,get_number_size};
 use self::formatter::tochar;
+use self::modexp::ModExp;
 use self::modmul::ModMul;
+use self::modsq::ModSquare;
 use self::mul::multiply;
 use self::shifts::{shiftl,shiftr};
+use self::square::Square;
 use self::sub::subtract;
 use std::cmp::{Ordering,min};
 use std::fmt;
@@ -57,10 +64,4 @@ macro_rules! base_impls
     }
 }
 
-macro_rules! modsq_impls
-{
-    ($name: ident) => {
-    }
-}
-
 include!("invoc.rs");

src/unsigned/modexp.rs

@@ -4,36 +4,60 @@ pub trait ModExp<T> {
 
 macro_rules! modexp_impls {
     ($name: ident) => {
-//        impl ModExp<$name> for $name {
-//            fn modexp(&self, e: &$name, m: &$name) -> $name {
-//                // S <- g
-//                let mut s = self.clone();
-//                // A <- 1
-//                let mut a = $name::from(1u64);
-//                // We do a quick skim through and find the highest index that
-//                // actually has a value in it.
-//                let mut e = ine.clone();
-//                // While e != 0 do the following:
-//                while e.values.iter().any(|x| *x != 0) {
-//                    // If e is odd then A <- A * S
-//                    if e.values[0] & 1 != 0 {
-//                        a.modmul(&s, m);
-//                    }
-//                    // e <- floor(e / 2)
-//                    let mut carry = 0;
-//                    e.values.iter_mut().rev().for_each(|x| {
-//                        let new_carry = *x & 1;
-//                        *x = (*x >> 1) | (carry << 63);
-//                        carry = new_carry;
-//                    });
-//                    // If e != 0 then S <- S * S
-//                    s.modsq(m);
-//                }
-//                // Return A
-//                a
-//            }
-//        }
+        impl ModExp<$name> for $name {
+            fn modexp(&self, ine: &$name, m: &$name) -> $name {
+                // S <- g
+                let mut s = self.clone();
+                // A <- 1
+                let mut a = $name::from(1u64);
+                // We do a quick skim through and find the highest index that
+                // actually has a value in it.
+                let mut e = ine.clone();
+                // While e != 0 do the following:
+                while e.value.iter().any(|x| *x != 0) {
+                    // If e is odd then A <- A * S
+                    if e.value[0] & 1 != 0 {
+                        a = a.modmul(&s, m);
+                        println!("Updating a to {:X}", a);
+                    }
+                    // e <- floor(e / 2)
+                    let mut carry = 0;
+                    e.value.iter_mut().rev().for_each(|x| {
+                        let new_carry = *x & 1;
+                        *x = (*x >> 1) | (carry << 63);
+                        carry = new_carry;
+                    });
+                    // If e != 0 then S <- S * S
+                    s = s.modsq(m);
+                }
+                // Return A
+                a
+            }
+        }
     };
     ($name: ident, $barrett: ident) => {
     };
 }
+
+#[cfg(test)]
+macro_rules! generate_modexp_tests {
+    ($name: ident, $lname: ident) => {
+        #[test]
+        fn $lname() {
+            let fname = format!("testdata/modexp/{}.tests", stringify!($name));
+            run_test(fname.to_string(), 4, |case| {
+                let (neg0, bbytes) = case.get("b").unwrap();
+                let (neg1, ebytes) = case.get("e").unwrap();
+                let (neg2, mbytes) = case.get("m").unwrap();
+                let (neg3, rbytes) = case.get("r").unwrap();
+                assert!(!neg0 && !neg1 && !neg2 && !neg3);
+
+                let b = $name::from_bytes(bbytes);
+                let e = $name::from_bytes(ebytes);
+                let m = $name::from_bytes(mbytes);
+                let r = $name::from_bytes(rbytes);
+                assert_eq!(r, b.modexp(&e, &m));
+            });
+        }
+    };
+}

src/unsigned/modmul.rs

@@ -14,3 +14,26 @@ macro_rules! modmul_impls {
         }
     };
 }
+
+#[cfg(test)]
+macro_rules! generate_modmul_tests {
+    ($name: ident, $lname: ident) => {
+        #[test]
+        fn $lname() {
+            let fname = format!("testdata/modmul/{}.tests", stringify!($name));
+            run_test(fname.to_string(), 4, |case| {
+                let (neg0, abytes) = case.get("a").unwrap();
+                let (neg1, bbytes) = case.get("b").unwrap();
+                let (neg2, mbytes) = case.get("m").unwrap();
+                let (neg3, cbytes) = case.get("c").unwrap();
+                assert!(!neg0 && !neg1 && !neg2 && !neg3);
+
+                let a = $name::from_bytes(abytes);
+                let b = $name::from_bytes(bbytes);
+                let m = $name::from_bytes(mbytes);
+                let c = $name::from_bytes(cbytes);
+                assert_eq!(c, a.modmul(&b, &m));
+            });
+        }
+    };
+}

src/unsigned/modsq.rs

@@ -0,0 +1,37 @@
+pub trait ModSquare<T> {
+    fn modsq(&self, m: &T) -> Self;
+}
+
+macro_rules! modsq_impls {
+    ($name: ident, $dbl: ident) => {
+        impl ModSquare<$name> for $name {
+            fn modsq(&self, m: &$name) -> $name {
+                let bigsquare = self.square();
+                let bigm      = $dbl::from(m);
+                let (_, res)  = bigsquare.divmod(&bigm);
+                $name::from(res)
+            }
+        }
+    };
+}
+
+#[cfg(test)]
+macro_rules! generate_modsq_tests {
+    ($name: ident, $lname: ident) => {
+        #[test]
+        fn $lname() {
+            let fname = format!("testdata/modsq/{}.tests", stringify!($name));
+            run_test(fname.to_string(), 3, |case| {
+                let (neg0, abytes) = case.get("a").unwrap();
+                let (neg1, mbytes) = case.get("m").unwrap();
+                let (neg2, cbytes) = case.get("c").unwrap();
+                assert!(!neg0 && !neg1 && !neg2);
+
+                let a = $name::from_bytes(abytes);
+                let m = $name::from_bytes(mbytes);
+                let c = $name::from_bytes(cbytes);
+                assert_eq!(c, a.modsq(&m));
+            });
+        }
+    };
+}

src/unsigned/square.rs

@@ -0,0 +1,72 @@
+pub trait Square<Output> {
+    fn square(&self) -> Output;
+}
+
+macro_rules! square_impls {
+    ($name: ident, $bigger: ident, $size: expr) => {
+        impl Square<$bigger> for $name {
+            fn square(&self) -> $bigger {
+               let mut w = [0; $size/32];
+               let     t = $size / 64;
+
+               for i in 0..t {
+                   let x128 = self.value[i] as u128;
+                   let mut uvb = (w[2*i] as u128) + (x128 * x128);
+                   w[2*i] = uvb & 0xFFFFFFFFFFFFFFFF;
+                   let mut c = uvb >> 64;
+                   for j in (i+1)..t {
+                       let xj128  = self.value[j] as u128;
+                       let xi128  = self.value[i] as u128;
+                       // this first product is safely 128 bits or less,
+                       // because the input arguments are both 64 bits.
+                       let xij128 = xj128 * xi128;
+                       // this next bit may overflow, but will do so by exactly
+                       // one bit.
+                       let twoxij128 = xij128 << 1;
+                       let carried_shl = (xij128 & (1 << 127)) != 0;
+                       // this next bit may *also* overflow, but should also do
+                       // so by no more than one bit.
+                       let (new,carry1) = twoxij128.overflowing_add(c);
+                       // ditto ...
+                       let wij = w[i+j];
+                       let (uvb2,carry2) = new.overflowing_add(wij as u128);
+                       // for the value we're going to save for this digit, we
+                       // only care about the low bits, so we can forget about
+                       // the carry stuff.
+                       w[i+j] = uvb2 & 0xFFFFFFFFFFFFFFFF;
+                       // for c, though, we do care about the carries, above.
+                       // Fortunately, they were both by only one bit, so we
+                       // should be able to just back-fix them.
+                       c = uvb2 >> 64;
+                       if carried_shl  { c += 1 << 64; }
+                       if carry1       { c += 1 << 64; }
+                       if carry2       { c += 1 << 64; }
+                   }
+                   w[i+t] = c;
+               }
+               let mut res = $bigger::zero();
+               for i in 0..w.len() { res.value[i] = w[i] as u64; }
+               res
+            }
+        }
+    };
+}
+
+#[cfg(test)]
+macro_rules! generate_square_tests {
+    ($name: ident, $lname: ident, $dbl: ident) => {
+        #[test]
+        fn $lname() {
+            let fname = format!("testdata/square/{}.tests", stringify!($name));
+            run_test(fname.to_string(), 2, |case| {
+                let (neg0, abytes) = case.get("a").unwrap();
+                let (neg1, rbytes) = case.get("r").unwrap();
+                assert!(!neg0 && !neg1);
+
+                let a = $name::from_bytes(abytes);
+                let r = $dbl::from_bytes(rbytes);
+                assert_eq!(r, a.square());
+            });
+        }
+    };
+}