Initial commit

bundled/Crypto/Number/Basic.hs

@@ -0,0 +1,116 @@
+-- |
+-- Module      : Crypto.Number.Basic
+-- License     : BSD-style
+-- Maintainer  : Vincent Hanquez <vincent@snarc.org>
+-- Stability   : experimental
+-- Portability : Good
+
+{-# LANGUAGE BangPatterns #-}
+module Crypto.Number.Basic
+    ( sqrti
+    , gcde
+    , areEven
+    , log2
+    , numBits
+    , numBytes
+    , asPowerOf2AndOdd
+    ) where
+
+import Data.Bits
+
+import Crypto.Number.Compat
+
+-- | @sqrti@ returns two integers @(l,b)@ so that @l <= sqrt i <= b@.
+-- The implementation is quite naive, use an approximation for the first number
+-- and use a dichotomy algorithm to compute the bound relatively efficiently.
+sqrti :: Integer -> (Integer, Integer)
+sqrti i
+    | i < 0     = error "cannot compute negative square root"
+    | i == 0    = (0,0)
+    | i == 1    = (1,1)
+    | i == 2    = (1,2)
+    | otherwise = loop x0
+        where
+            nbdigits = length $ show i
+            x0n = (if even nbdigits then nbdigits - 2 else nbdigits - 1) `div` 2
+            x0  = if even nbdigits then 2 * 10 ^ x0n else 6 * 10 ^ x0n
+            loop x = case compare (sq x) i of
+                LT -> iterUp x
+                EQ -> (x, x)
+                GT -> iterDown x
+            iterUp lb = if sq ub >= i then iter lb ub else iterUp ub
+                where ub = lb * 2
+            iterDown ub = if sq lb >= i then iterDown lb else iter lb ub
+                where lb = ub `div` 2
+            iter lb ub
+                | lb == ub   = (lb, ub)
+                | lb+1 == ub = (lb, ub)
+                | otherwise  =
+                    let d = (ub - lb) `div` 2 in
+                    if sq (lb + d) >= i
+                        then iter lb (ub-d)
+                        else iter (lb+d) ub
+            sq a = a * a
+
+-- | Get the extended GCD of two integer using integer divMod
+--
+-- gcde 'a' 'b' find (x,y,gcd(a,b)) where ax + by = d
+--
+gcde :: Integer -> Integer -> (Integer, Integer, Integer)
+gcde a b = onGmpUnsupported (gmpGcde a b) $
+    if d < 0 then (-x,-y,-d) else (x,y,d)
+  where
+    (d, x, y)                     = f (a,1,0) (b,0,1)
+    f t              (0, _, _)    = t
+    f (a', sa, ta) t@(b', sb, tb) =
+        let (q, r) = a' `divMod` b' in
+        f t (r, sa - (q * sb), ta - (q * tb))
+
+-- | Check if a list of integer are all even
+areEven :: [Integer] -> Bool
+areEven = and . map even
+
+-- | Compute the binary logarithm of a integer
+log2 :: Integer -> Int
+log2 n = onGmpUnsupported (gmpLog2 n) $ imLog 2 n
+  where
+    -- http://www.haskell.org/pipermail/haskell-cafe/2008-February/039465.html
+    imLog b x = if x < b then 0 else (x `div` b^l) `doDiv` l
+      where
+        l = 2 * imLog (b * b) x
+        doDiv x' l' = if x' < b then l' else (x' `div` b) `doDiv` (l' + 1)
+{-# INLINE log2 #-}
+
+-- | Compute the number of bits for an integer
+numBits :: Integer -> Int
+numBits n = gmpSizeInBits n `onGmpUnsupported` (if n == 0 then 1 else computeBits 0 n)
+  where computeBits !acc i
+            | q == 0 =
+                if r >= 0x80 then acc+8
+                else if r >= 0x40 then acc+7
+                else if r >= 0x20 then acc+6
+                else if r >= 0x10 then acc+5
+                else if r >= 0x08 then acc+4
+                else if r >= 0x04 then acc+3
+                else if r >= 0x02 then acc+2
+                else if r >= 0x01 then acc+1
+                else acc -- should be catch by previous loop
+            | otherwise = computeBits (acc+8) q
+          where (q,r) = i `divMod` 256
+
+-- | Compute the number of bytes for an integer
+numBytes :: Integer -> Int
+numBytes n = gmpSizeInBytes n `onGmpUnsupported` ((numBits n + 7) `div` 8)
+
+-- | Express an integer as an odd number and a power of 2
+asPowerOf2AndOdd :: Integer -> (Int, Integer)
+asPowerOf2AndOdd a
+    | a == 0       = (0, 0)
+    | odd a        = (0, a)
+    | a < 0        = let (e, a1) = asPowerOf2AndOdd $ abs a in (e, -a1)
+    | isPowerOf2 a = (log2 a, 1)
+    | otherwise    = loop a 0
+        where      
+          isPowerOf2 n = (n /= 0) && ((n .&. (n - 1)) == 0)
+          loop n pw = if n `mod` 2 == 0 then loop (n `div` 2) (pw + 1)
+                      else (pw, n)