{-# LANGUAGE BangPatterns #-} -- | -- Module : Number.Basic -- License : BSD-style -- Maintainer : Vincent Hanquez -- Stability : experimental -- Portability : Good module Number.Basic ( sqrti , gcde , gcde_binary , areEven ) where import Data.Bits -- | sqrti returns two integer (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 :: Integer -> Integer -> (Integer, Integer, Integer) gcde 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)) -- | get the extended GCD of two integer using the extended binary algorithm (HAC 14.61) -- get (x,y,d) where d = gcd(a,b) and x,y satisfying ax + by = d gcde_binary :: Integer -> Integer -> (Integer, Integer, Integer) gcde_binary a' b' | b' == 0 = (1,0,a') | a' >= b' = compute a' b' | otherwise = (\(x,y,d) -> (y,x,d)) \$ compute b' a' where getEvenMultiplier !g !x !y | areEven [x,y] = getEvenMultiplier (g `shiftL` 1) (x `shiftR` 1) (y `shiftR` 1) | otherwise = (x,y,g) halfLoop !x !y !u !i !j | areEven [u,i,j] = halfLoop x y (u `shiftR` 1) (i `shiftR` 1) (j `shiftR` 1) | even u = halfLoop x y (u `shiftR` 1) ((i + y) `shiftR` 1) ((j - x) `shiftR` 1) | otherwise = (u, i, j) compute a b = let (x,y,g) = getEvenMultiplier 1 a b in loop g x y x y 1 0 0 1 loop g _ _ 0 !v _ _ !c !d = (c, d, g * v) loop g x y !u !v !a !b !c !d = let (u2,a2,b2) = halfLoop x y u a b in let (v2,c2,d2) = halfLoop x y v c d in if u2 >= v2 then loop g x y (u2 - v2) v2 (a2 - c2) (b2 - d2) c2 d2 else loop g x y u2 (v2 - u2) a2 b2 (c2 - a2) (d2 - b2) -- | check if a list of integer are all even areEven :: [Integer] -> Bool areEven = and . map even