101 lines
3.5 KiB
Haskell
101 lines
3.5 KiB
Haskell
-- |
|
|
-- 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
|
|
) where
|
|
|
|
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)
|