text:
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}

import Math.NumberTheory.Primes.Testing
import Math.NumberTheory.Moduli
import Options.Applicative
import System.Random
import System.IO
import System.Exit
import Data.Bits
import Data.Function
import Data.List
import Data.Maybe
import Test.QuickCheck
import Debug.Trace
import Control.Monad
import Statistics.Sample
import Statistics.Quantile
import Control.Exception (evaluate)
import Data.Bifunctor
import GHC.Stack
import Text.Read
import Text.Printf
import Data.List.Split
import qualified Data.Vector as V
import qualified Data.Sequence as S
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.Monoid

-- Some convenient names for some parameters
type BitLength = Int
type Width = Int

-- http://stackoverflow.com/a/13301611/946226
rndPrime :: BitLength -> Gen Integer
rndPrime 1 = error "Cannot create distinct p and q for 1 bit"
rndPrime bits = do
    x <- fmap (.|. 1) $ choose (2^(bits - 1), 2^bits - 1)
    if isPrime x then return x else rndPrime bits

data RSAParams = RSAParams { bits :: Int, p, q, n, e, dp, dq, kp, kq :: Integer } deriving Show

instance Arbitrary RSAParams where
    arbitrary = sized (\s -> genRSAParams (min 1024 (max s 2)))

genRSAParams :: BitLength -> Gen RSAParams
genRSAParams bits = do
    p <- rndPrime bits
    q <- rndPrime bits
    let n = p*q
        e = 65537
        mdp = invertMod e (p-1)
        mdq = invertMod e (q-1)
    case (mdp, mdq) of
        (Just dp, Just dq) -> do
            let kp = (e*dp-1) `div` (p-1)
                kq = (e*dq-1) `div` (q-1)
            return (RSAParams {..})
        _ -> genRSAParams bits -- retry

data Step = S | M deriving (Eq, Ord, Show) -- square or multiply
type Leak = V.Vector Step

leftToRight :: BitLength -> Width -> Integer -> Leak
leftToRight bits w s | w <= 0 = error "width must be positive"
leftToRight bits w s = V.fromList $ go (bits-1)
   where
    go i | i < 0         = []
    go i | s `testBit` i = replicate (slideLength-1) S ++ M : go (i - slideLength)
      where
        Just slideLength = find (\l -> s `testBit` (i-l+1)) $ [w,w-1..1]
    go i | otherwise     = S : go (i-1)

-- Bits learned
leftToRightGMPBits :: BitLength -> Width -> Integer -> Integer
leftToRightGMPBits bits w s | w <= 0 = error "width must be positive"
leftToRightGMPBits bits w s = go (bits -1) 0
   where
    go i !acc | i < 0         = acc
    go i !acc | s `testBit` i = go (i - slideLength) (acc + learned)
      where
        Just slideLength = find (\l -> s `testBit` (i-l+1)) $ [w,w-1..1]
        learned | slideLength == 1 = 1
                | otherwise        = 2
    go i ! acc | otherwise    =  go (i-1) (acc + 1)

rightToLeft :: BitLength -> Width -> Integer -> Leak
rightToLeft bits w s | w <= 0 = error "width must be positive"
rightToLeft bits w s = V.fromList $ go 0 []
   where
    go i acc | i == bits     = acc
    go i acc | s `testBit` i = go (i+slideLength) (replicate (slideLength-1) S ++ M : acc)
      where slideLength = min w (bits - i)
    go i acc | otherwise     = go (i+1) (S : acc)

-- Bits learned
rightToLeftBits :: BitLength -> Width -> Integer -> Integer
rightToLeftBits bits w s | w <= 0 = error "width must be positive"
rightToLeftBits bits w s = go 0 0
   where
    go i !acc | i == bits     = acc
    go i !acc | s `testBit` i = go (i+slideLength) (1+acc)
      where slideLength = min w (bits - i)
    go i !acc | otherwise     = go (i+1) (1+acc)

-- A datatype to track known bits
data BitMask = BM { unknownBits :: !Integer, knownBits :: !Integer }
    deriving Eq

showBitMask bits bm = map go [(bits-1),(bits-2)..0]
  where go i | not (bm `isKnownBit` i) = '_'
             | bm `getKnownBit` i      = '1'
             | otherwise               = '0'

instance Show BitMask where
    show bm = showBitMask bits bm
      where bits = head [i | i <- [0..], 2^i > knownBits bm && 2^i > unknownBits bm]

isKnownBit :: BitMask -> Int -> Bool
bm `isKnownBit` i = not (unknownBits bm `testBit` i)

getKnownBit :: BitMask -> Int -> Bool
bm `getKnownBit` i = knownBits bm `testBit` i

possibleBits :: BitMask -> Int -> [Bool]
possibleBits bm i
    | bm `isKnownBit` i = [bm `getKnownBit` i]
    | otherwise         = [False, True]

matches :: Integer -> BitMask -> Bool
matches p bm = (p .|. unknownBits bm) == (knownBits bm .|. unknownBits bm)

isKnown :: BitMask -> Maybe Integer
isKnown bm | unknownBits bm == 0 = Just (knownBits bm)
           | otherwise           = Nothing

writeBit :: Int -> Bool -> Integer -> Integer
writeBit i True  p = p `setBit` i
writeBit i False p = p `clearBit` i

writeKnownBit :: Int -> Bool -> BitMask -> BitMask
writeKnownBit i b bm = BM (unknownBits bm `clearBit` i) (writeBit i b (knownBits bm))

nothingKnown :: BitLength -> BitMask
nothingKnown bits = BM (2^bits-1) 0

anyOddNumber :: BitLength -> BitMask
anyOddNumber bits = writeKnownBit 0 True (nothingKnown bits)

fromKnownNumber :: Integer -> BitMask
fromKnownNumber i = BM 0 i

-- can be optimized, I guess
allMatching :: BitMask -> [Integer]
allMatching bm = map sum $ mapM go [(bits-1),(bits-2)..0]
  where go i | not (bm `isKnownBit` i) = [2^i,0]
             | bm `getKnownBit` i      = [2^i]
             | otherwise               = [0]
        bits = head [i | i <- [0..], 2^i > knownBits bm && 2^i > unknownBits bm]

maskFromMatching :: [Integer] -> BitMask
maskFromMatching (m:ms) = foldl' go (fromKnownNumber m) ms
  where
    go :: BitMask -> Integer -> BitMask
    go bm i = BM (unknownBits bm .|. disagree) (knownBits bm `clear` disagree)
      where disagree = (unknownBits bm .|. knownBits bm) `xor` (unknownBits bm .|. i)

x `clear` y = complement (complement x .|. y)

countKnownBits :: BitLength -> BitMask -> Int
countKnownBits bits bm = bits - popCount (unknownBits bm)

-- A breadth first search limiting the count by n
bfSearch :: Int -> a -> (a -> [a]) -> (a -> Bool) -> Maybe Int
bfSearch limit s0 next done = go 0 (S.singleton s0)
  where
    go n _ | n >= limit = Nothing
    go n todo = case S.viewl todo of
        S.EmptyL -> Nothing
        s S.:< todo' | done s    -> Just n
                     | otherwise -> go (n+1) (todo' S.>< S.fromList (next s))

dfSearch :: Int -> a -> (a -> [a]) -> (a -> Bool) -> Maybe Int
dfSearch limit s0 next done = go 0 [s0]
  where
    go n _ | n >= limit = Nothing
    go n [] = Nothing
    go n (s:todo') | done s    = Just n
                   | otherwise = go (n+1) (next s ++ todo')

-- Currently unused
linearWorkList :: forall a. Eq a =>
    IS.IntSet -> IM.IntMap a -> (Int -> IM.IntMap a -> a) -> IM.IntMap a
linearWorkList todo s0 recalc = go todo s0
  where
    go :: IS.IntSet -> IM.IntMap a -> IM.IntMap a
    go todo s | IS.null todo         = s
    go todo s | oldValue == newValue = go todo' s
              | otherwise            = go newTodo newS
      where
        (k, todo') = IS.deleteFindMin todo
        oldValue = s IM.! k
        newValue = recalc k s
        newS = IM.insert k newValue s
        newTodo = (if (k-1) `IM.member` s then IS.insert (k-1) else id) $
                  (if (k+1) `IM.member` s then IS.insert (k+1) else id) $ todo'

type WindowLengthRanges = [(Int, (Int,Int))]

possibleWindowLengths :: BitLength -> Width -> Leak -> WindowLengthRanges
possibleWindowLengths bits w leak = IM.elems result
  where
    clamp width = max 1 (min w width)

    -- positions of the window end (from the back, i.e. bits)
    windowEnds = IM.fromList $ zip [0..] $ map ((bits - 1) - ) $
                 V.toList $ V.findIndices (== M) leak

    windows = IM.keysSet windowEnds
    maxKey = IS.findMax windows

    zipIM = IM.intersectionWith (,)

    windowMax = IM.mapWithKey (\k _ -> w_max k) windowEnds
    windowMin = IM.mapWithKey (\k _ -> w_min k) windowEnds
    result = windowEnds `zipIM` (windowMin `zipIM` windowMax)

    msbIsOne = False

    -- The first window spans precisely from the MSB to the first M
    firstWinLength =  clamp (bits - pos)
      where pos = windowEnds IM.! 0

    -- minimum length. computed from right to left
    w_min k | k == 0, msbIsOne = firstWinLength
            | k == maxKey      = 1
            | otherwise        = clamp (j + w_min_j + w -i)
        where
            i       = windowEnds IM.! k
            j       = windowEnds IM.! (k+1)
            w_min_j = windowMin IM.! (k+1)

    -- maximum length. computed from left to right
    w_max k | k == 0    = if msbIsOne then firstWinLength else w
            | otherwise = clamp (i + w_max_i - w -j)
        where
            i       = windowEnds IM.! (k-1)
            j       = windowEnds IM.! k
            w_max_i = windowMax IM.! (k-1)

countPossibilities :: BitLength -> Width -> WindowLengthRanges -> Integer
countPossibilities bits w lengths = go (bits-1) 0
  where
    vec = V.fromList lengths

    pos k = fst (vec V.! k)
    windowLengthRange k = [minWin..maxWin]
      where
        (_, (minWin, maxWin)) =  vec V.! k

    go soonestOne k | k == V.length vec = 2^0
    go soonestOne k = sum
        [ countMemo k winLength
        | winLength <- windowLengthRange k
        , let firstOneOfWindow = pos k + winLength -1
        , firstOneOfWindow <= soonestOne
        ]

    -- simple memoization
    countTable = V.fromList
        [ count k winLength | k <- [0..V.length vec -1] , winLength <- [1..w] ]
    countMemo k winLength = countTable V.!  (k * w + (winLength - 1))

    count k winLength = 2^unknown * go nextPossibleOne (k+1)
      where
        firstOneOfWindow = pos k + winLength -1
        unknown = max 0 (winLength - 2)
        nextPossibleOne = firstOneOfWindow - w

testCountPossibilities :: BitLength -> Width -> Integer -> Property
testCountPossibilities bits w s = fromIntegral (length candidates) === poss
  where
    leak = leftToRight bits w s
    windowLengths = possibleWindowLengths bits w leak
    poss = countPossibilities bits w windowLengths

    mask = getHardBits bits w windowLengths
    candidates = realCandidates bits w leak mask

getWindowStarts :: WindowLengthRanges -> IS.IntSet
getWindowStarts list = IS.fromList
    [ pos + maxWin - 1
    | (pos, (minWin, maxWin)) <- list
    ]

getHardBits :: BitLength -> Width -> WindowLengthRanges -> BitMask
getHardBits bits w final_list = go 0 (nothingKnown bits)
      where
        go i !bm | i == bits = bm
                 | i `IS.member` ones    = go (i+1) (writeKnownBit i True bm)
                 | i `IS.member` zeros   = go (i+1) (writeKnownBit i False bm)
                 | i `IS.member` inRange = go (i+1) bm
                 | otherwise             = go (i+1) (writeKnownBit i False bm)

        onesOn =     IS.fromList [ pos
                                 | (pos,(minWin, maxWin)) <- final_list ]
        onesStart =  IS.fromList [ pos + minWin - 1
                                 | (pos,(minWin, maxWin)) <- final_list
                                 , minWin == maxWin ]
        zerosAfter = IS.fromList [ j
                                 | (pos,(minWin, maxWin)) <- final_list
                                 , j <- [pos + maxWin - w .. pos - 1] ]
        inRange =    IS.fromList [ j
                                 | (pos,(minWin, maxWin)) <- final_list
                                 , j <- [pos .. pos + maxWin - 1] ]

        ones  = IS.unions [ onesOn, onesStart ]
        zeros = IS.unions [ zerosAfter ]

countHardBits :: BitLength -> Width -> WindowLengthRanges -> Int
countHardBits bits w final_list = bits - sum (map uncertain final_list)
  where uncertain (_, (minWin, maxWin)) | minWin == maxWin = max 0 (maxWin - 2)
                                        | otherwise        = max 0 (maxWin - 1)

testCountHardBits bits w final_list =
    countKnownBits bits (getHardBits bits w final_list) == countHardBits bits w final_list

knownBitsFromLeak :: BitLength -> Width -> Leak -> BitMask
knownBitsFromLeak bits w leak = getHardBits bits w $ possibleWindowLengths bits w leak

testKnownBitsFromLeak :: BitLength -> Positive Integer -> Bool
testKnownBitsFromLeak bits (Positive p) = p `matches` bm
  where
    leak = leftToRight bits 4 p
    bm = knownBitsFromLeak bits 4 leak

testKnownBitsFromLeak2 :: RSAParams -> Bool
testKnownBitsFromLeak2 params
    = dp params `matches` bitmask_dp && dq params `matches` bitmask_dq
  where
    leak_dp    = leftToRight       (bits params) 4 (dp params)
    bitmask_dp = knownBitsFromLeak (bits params) 4 leak_dp
    leak_dq    = leftToRight       (bits params) 4 (dq params)
    bitmask_dq = knownBitsFromLeak (bits params) 4 leak_dq

type Candidate = (Int, Integer, Integer, Integer, Integer)

-- | An oracle is a function that takes a possible suffix (given with length
-- and value) and indicates if this suffix is possible.
type Oracle = Int -> Integer -> Bool

selectPortion :: Double -> Int -> [Int]
selectPortion h to = go 0 0
  where
    go chosen i | i > to    = []
                | f < h     = i : go (chosen + 1) (i+1)
                | otherwise =     go  chosen      (i+1)
      where f = fromIntegral chosen / fromIntegral i

selectRandomPortion :: StdGen -> Double -> Int -> [Int]
selectRandomPortion rng h to = go rng 0
  where
    go rng i | i > to    = []
             | x < h     = i : go rng' (i+1)
             | otherwise =     go rng' (i+1)
      where (x, rng') = random rng

mkHardBitOracle :: BitLength -> Double -> Integer -> Oracle
mkHardBitOracle bits h input = oracle
  where
    rng = mkStdGen (fromIntegral (input+3))
    knownBits = IS.fromList $ selectRandomPortion rng h bits
    oracle i guess
        | i `IS.member` knownBits
        = guess `testBit` i == input `testBit` i
        | otherwise
        = True

mkBlockSumOracle :: Int -> Integer -> Oracle
mkBlockSumOracle bs input = oracle
  where
    oracle i guess | j `mod` bs == 0 = guess ==? input
                   | otherwise       = True
      where
        j = i + 1
        select k = (k `div` 2^(j-bs)) `mod` 2^bs
        (==?) = ((==) `on` popCount) `on` select

testOracle :: (Integer -> Oracle) -> Positive Integer -> Positive Int ->  Bool
testOracle mkO (Positive input) (Positive i) = o i (input `mod` 2^(i+1))
  where o = mkO input

mkLeftToRightOracle bits w input = oracle
  where
   s = leftToRight bits w input

   windowLengths = possibleWindowLengths bits w s
   bitmask = getHardBits bits w windowLengths
   windowStarts = getWindowStarts windowLengths

   shouldCheck i = (i-1) `IS.member` windowStarts

   oracle i guess
        | (guess `testBit` i) `notElem` possibleBits bitmask i
        = False
        | shouldCheck i
        , let leak' = leftToRight i w guess
        , V.drop (bits - i) s /= leak'
        = False
        | otherwise
        = True

-- Recovers the bits from the leaks for p and q
lsbRecovery :: Int -> RSAParams -> (Oracle, Oracle) -> Maybe Int
lsbRecovery limit params (op, oq) = dfSearch limit s0 next done
  where
    s0 :: Candidate
    s0 = ( 0 -- last bit guessed/set; this and all later bits are known
         , 1 -- known suffix of dp
         , 1 -- known suffix of dq
         , 1 -- known suffix of p
         , 1 -- known suffix of q
         )

    next :: Candidate -> [Candidate]
    next (i, dp, dq, p, q)
        | j == (bits params) = []
        | otherwise =
            filter prune
            [ (j, dp', dq', p', q')
            | (p',dp') <- filter (pruneP j) $ (,) <$> nextP <*> nextDP
            , (q',dq') <- filter (pruneQ j) $ (,) <$> nextQ <*> nextDQ
            ]
      where
        j = i+1
        nextDP = [dp, dp `setBit` j]
        nextDQ = [dq, dq `setBit` j]
        nextP  = [p, p `setBit` j]
        nextQ  = [q, q `setBit` j]

    pruneP i (p, dp) = and
        [ (e params*dp - 1) `mod` modulus == (kp params * (p-1)) `mod` modulus
        , op i dp
        ]
      where modulus = 2^(i+1)

    pruneQ i (q, dq) = and
        [ (e params*dq - 1) `mod` modulus == (kq params * (q-1)) `mod` modulus
        , oq i dq
        ]
      where modulus = 2^(i+1)

    prune :: Candidate -> Bool
    prune (i, dp, dq, p, q)
        = and [ (p * q) `mod` modulus == n params `mod` modulus
              ]
      where modulus = 2^(i+1)

    done :: Candidate -> Bool
    done (i, dp', dq', p', q')
        = i + 1 == bits params
        && p' * q' == n params
        {-
        && Just dp' == invertMod (e params) (p'-1)
        && Just dq' == invertMod (e params) (q'-1)
        && (dp' == dp params || traceShow ('P',p',q', dp', dq', params) True)
        && (dq' == dq params || traceShow ('Q',p',q', dp', dq', params) True)
        && (dp' `matches` bitmask_dp || traceShow (bits params, dp params) True)
        && (dq' `matches` bitmask_dq || traceShow (bits params, dq params) True)
        -}

runLSBRecovery :: Int -> Width -> RSAParams -> Maybe Int
runLSBRecovery limit w params =
    let op = mkLeftToRightOracle (bits params) w (dp params)
        oq = mkLeftToRightOracle (bits params) w (dq params)
    in lsbRecovery limit params (op, oq)

main :: IO ()
main = join . customExecParser (prefs showHelpOnEmpty) $
  info (helper <*> parser)
  (  fullDesc
  <> header "Sliding right experiment"
  <> progDesc "Slides stuff around"
  <> footer "Some parameter understand range syntax in the form 1,2,3,4 or 4,8,..,2048"
  )
  where
    parser :: Parser (IO ())
    parser =
      hsubparser (mconcat
        [ command "run_recover"        recoverCommand
        , command "left_to_right"      printLTRCommand
        , command "test_hard_rules"    testHardRulesCommand
        , command "measure_leakage"    measureLeakageCommand
        , command "collision_entropy"  measureCollisionEntropyCommand
        , command "hard_soft_scatter"  hardSoftScatterCommand
        , command "leak_preinput_size" leakPreinputCommand
        , command "entropy_scatter"    entropyScatterCommand
        ])

    bitsArg = (`div` 2) <$> -- internally, we want the key length of p and q
      option auto
      (  long "bits"
         <> metavar "BITS"
         <> help "Bit length of p and q combined (usually 1024 or 2048)"
         <> value 1024
         <> showDefault
      )

    bitsRangeArg = (map (`div` 2)) <$> -- internally, we want the key length of p and q
      option (eitherReader parseRange)
      (  long "bits"
         <> metavar "BITS,.."
         <> help "Bit length of p and q combined. Range syntax allowed"
         <> value [1024]
         <> showDefaultWith showRange
      )

    knownBitsRangeArg =
      option (eitherReader parseRange)
      (  long "known-bits"
         <> metavar "PROB,.."
         <> help "Proportion of known bits in range 0...1000. Range syntax allowed"
         <> value [0,100..1000]
         <> showDefaultWith showRange
      )

    blockSumRangeArg =
      option (eitherReader parseRange)
      (  long "sum-length"
         <> metavar "INT,.."
         <> help "Lenght of blocks of which sum is known"
         <> value [1,2,3,4]
         <> showDefaultWith showRange
      )

    widthArg = option auto
         (  long "width"
         <> metavar "WIDTH"
         <> help "Width of the window (e.g. 4 or 5)"
         <> value 4
         <> showDefault
         )

    widthRangeArg = option (eitherReader parseRange)
         (  long "width"
         <> metavar "WIDTH,.."
         <> help "Width of the window (e.g. 4 or 5). Range syntax allowed."
         <> value [4]
         <> showDefaultWith showRange
         )

    iterationsArg = option auto
         (  long "iterations"
         <> metavar "N"
         <> help "Number of iterations"
         <> value 100
         <> showDefault
         )

    limitArg = option auto
         (  long "limit"
         <> metavar "N"
         <> help "Search limit"
         <> value 100000
         <> showDefault
         )

    tsvArg = switch
         (  long "tsv"
         <> help "Tab separated output (no variance)"
         )

    printLTRCommand = info
        (pure printLTR <*> bitsArg <*> widthArg<*> bitSequenceArg)
        (progDesc "Run left-to-right and print sequence")
      where bitSequenceArg = argument auto (metavar "input_sequence")

    testHardRulesCommand = info
        (pure testHardRules <*> bitsArg <*> widthArg<*> verboseArg)
        (progDesc "Verify hard bit calculation (no output = all good)")
      where verboseArg = switch (long "verbose" <> help "Show all tested data")

    recoverCommand = info
        (pure runRecover <*> bitsArg <*> widthArg<*> iterationsArg <*> limitArg )
        (progDesc "Run key recovery, and calculate search tree size")

    measureLeakageCommand = info
        (pure measureLeakage <*> bitsRangeArg <*> widthRangeArg<*> tsvArg <*> iterationsArg)
        (progDesc "Calculate the Shannon entropy rate and hard bits")

    measureCollisionEntropyCommand = info
        (pure measureCollisionEntropy <*> bitsRangeArg <*> widthRangeArg<*> tsvArg <*> iterationsArg)
        (progDesc "Calculate the collision entropy rate")

    hardSoftScatterCommand = info
        (pure hardSoftScatter <*> bitsRangeArg <*> widthRangeArg<*> iterationsArg)
        (progDesc "Draw a hard/soft bit scatter plot")

    leakPreinputCommand = info
        (pure leakPreinput <*> widthArg)
        (progDesc "Number of input sequences for given square-and-mulitply sequence")

    entropyScatterCommand = info
        (pure runScatter <*> bitsArg <*> knownBitsRangeArg <*> widthRangeArg <*> blockSumRangeArg <*>iterationsArg <*> limitArg )
        (progDesc "Simulates recovery with various oracles of given entropy, measuing tree size")

runRecover :: BitLength -> Width -> Int -> Int -> IO ()
runRecover  bits w iter limit = do
    samples <- replicateM iter $ do
       params <- generate (genRSAParams bits)
       let res = runLSBRecovery limit w params
       putStrLn $ maybe "-" show res
       return res
    let (mean, variance) = meanVariance (V.fromList (map fromIntegral (catMaybes samples)))
    putStrLn $ "-- Mean: " ++ show mean ++ " Variance: " ++ show variance

popCountCollEntropyRate :: Int -> Double
popCountCollEntropyRate w = (/fromIntegral w) $ negate $ logBase 2 $ sum [
    (fromIntegral (binom w i) / 2^w)^2
    | i <- [0..w]
    ]

binom :: Int -> Int -> Int
binom n m = fac n `div` (fac m * fac (n-m))

fac :: Int -> Int
fac n = product [1..n]

runScatter :: BitLength -> [Int] -> [Width] -> [Int] -> Int -> Int -> IO ()
runScatter  bits knownBitsRatios widths blockSums iter limit = do
    printf "leak\th\tmean\tsolved\tcalc_mean\n"
    let todo = [ (p, o, "hard_bits")
               | ratio <- sort knownBitsRatios
               , let p = fromIntegral ratio / 1000 :: Double
               , let o = mkHardBitOracle bits p
               ] ++
               [ (h, o, "ltr")
               | w <- sortBy (flip compare) widths
               , let h = [undefined, 1, 0.786, 0.652, 0.545, 0.461, 0.395, 0.345] !! w
               , let o = mkLeftToRightOracle bits w
               ] ++
               [ (h, o, "popcnt")
               | bs <- sortBy (flip compare) blockSums
               , let h = popCountCollEntropyRate bs
               , let o = mkBlockSumOracle bs
               ]
    forM_ todo $ \(h, mkO, desc) -> do
        samples <- replicateM iter $ do
           params <- generate (genRSAParams bits)
           let op = mkO (dp params)
               oq = mkO (dq params)
               res = lsbRecovery limit params (op, oq)
           -- putStrLn $ maybe "-" show res
           return res
        let (mean, variance) = meanVariance (V.fromList (map fromIntegral (catMaybes samples)))
        let solved = fromIntegral (length (filter isJust samples)) / fromIntegral iter :: Double
        printf "%s\t%.3f\t%.1f\t%.3f\t%.1g\n"
            desc h mean solved (treeSize bits h)

treeSize :: BitLength -> Double -> Double
treeSize k h = k' * (1 + (1 - 2**(k'*(1-2*h)))/(1-2**(1-2*h)))
  where k' = 2*fromIntegral k

printLargePrime :: BitLength -> Width -> IO ()
printLargePrime bits w = generate (rndPrime bits) >>= print

printRSA :: BitLength -> Width -> IO ()
printRSA bits w = generate (genRSAParams bits) >>= print

showLeak :: Leak -> String
showLeak = map go . V.toList
  where go S = '_'
        go M = 'M'

printLTR :: BitLength -> Width -> Integer -> IO ()
printLTR bits w s = putStrLn $ showLeak $ leftToRight bits w s

realCandidates :: BitLength -> Width -> Leak -> BitMask -> [Integer]
realCandidates bits w leak mask = ms'
   where ms    = allMatching mask
         ms'   = filter (\m -> leftToRight bits w m == leak) ms

reallyKnownBits :: BitLength -> Width -> Leak -> BitMask -> BitMask
reallyKnownBits bits w leak mask = mask'
   where ms'   = realCandidates bits w leak mask
         mask' = maskFromMatching ms'

testKnownBits bits w n = do
    putStrLn (showBitMask bits (fromKnownNumber n))
    let leak = leftToRight bits w n
    putStrLn (showLeak leak)
    let mask = knownBitsFromLeak bits w leak
    putStrLn (showBitMask bits mask)
    let mask' = reallyKnownBits bits w leak mask
    putStrLn (showBitMask bits mask')

testHardRules :: BitLength -> Width -> Bool -> IO ()
testHardRules bits w verbose = forM_ [2^(bits-1)..(2^bits)-1] $ \n -> do
    let leak  = leftToRight bits w n
    let mask  = knownBitsFromLeak bits w leak
    let ms    = allMatching mask
    let ms'   = filter (\m -> leftToRight bits w m == leak) ms
    let mask' = maskFromMatching ms'

    when (verbose || mask /= mask' || n `notElem` ms) $ do
        putStr $ showBitMask bits (fromKnownNumber n) ++ ": "
        putStr $ showLeak leak ++ " "
        putStr $ showBitMask bits mask ++ " "
        putStr $ showBitMask bits mask' ++ " "
        putStr $ if mask == mask' then "✓" else "☹"
        unless (n `elem` ms) $ putStr " HORROR"
        putStr "\n"

measureCollisionEntropy :: [BitLength] -> [Width] -> Bool -> Int -> IO ()
measureCollisionEntropy bitsRange wRange tsv iter= do
    when tsv $ putStrLn $ intercalate "\t" $
        words "bits width entropy"
    sequence_ [ measureCollisionEntropyOne bits w tsv iter
              | bits <- bitsRange, w <- wRange ]

measureCollisionEntropyOne :: BitLength -> Width -> Bool -> Int -> IO ()
measureCollisionEntropyOne bits w tsv iter = do
    probs <- replicateM iter doOne
    -- There are serious problems here with numerical precision
    -- and the results seem to vary a log depending on how the logarithm
    -- is calculated
    let entropy = - (logBase 2.0 (sum probs) - fromIntegral (2*bits) - logBase 2.0 (fromIntegral iter)) / fromIntegral (2*bits)
    if tsv then
        putStrLn $ intercalate "\t" $ map show [2*bits, w] ++ map show [entropy]
    else putStrLn $ "Collision entropy rate:" ++ show entropy
  where
    stats :: V.Vector Double -> [Double]
    stats samples =
        [ minimum samples ] ++
        [ weightedAvg n 4 samples | n <- [1..3]] ++
        [ maximum samples ]
    doOne = do
        dp <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)
        dq <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)

        let sp = leftToRight bits w dp
            wlp = possibleWindowLengths bits w sp
            hardBitsp = countHardBits bits w wlp
            possp = countPossibilities bits w wlp

            sq = leftToRight bits w dq
            wlq = possibleWindowLengths bits w sq
            hardBitsq = countHardBits bits w wlq
            possq = countPossibilities bits w wlq

            totalBits = fromIntegral (2*bits)
            prob = fromIntegral possp * fromIntegral possq
        return prob

measureLeakage :: [BitLength] -> [Width] -> Bool -> Int -> IO ()
measureLeakage bitsRange wRange tsv iter= do
    when tsv $ putStrLn $ intercalate "\t" $
        words "bits width" ++
        concatMap (\n -> [n ++ "_q" ++ show i | i <- [0..4]])
            (words "rightToLeft hardBitsLearned informationGained GMP")
    sequence_ [ measureLeakageOne bits w tsv iter
              | bits <- bitsRange, w <- wRange ]

measureLeakageOne :: BitLength -> Width -> Bool -> Int -> IO ()
measureLeakageOne bits w tsv iter= do
    (samplesRTL, samplesHard, samplesReal, samplesGMP) <- V.unzip4 . V.fromList <$> replicateM iter doOne
    if tsv then do
        putStrLn $ intercalate "\t" $
            map show [2*bits, w] ++
            map show (concatMap stats [ samplesRTL, samplesHard, samplesReal, samplesGMP])
    else do
        let (mean, variance) = meanVariance samplesHard
        let good = V.sum (V.map (\n -> if n > 0.5 then 1 else 0) samplesHard) / fromIntegral iter
        putStrLn $ "Hard bits mean:  " ++ show mean ++
                   " variance: " ++ show variance ++
                   " >50%: " ++ show good
        let (mean, variance) = meanVariance samplesReal
        let good = V.sum (V.map (\n -> if n > 0.5 then 1 else 0) samplesReal) / fromIntegral iter
        putStrLn $ "Total data mean: " ++ show mean ++
                   " variance: " ++ show variance ++
                   " >50%: " ++ show good
  where
    stats :: V.Vector Double -> [Double]
    stats samples =
        [ minimum samples ] ++
        [ weightedAvg n 4 samples | n <- [1..3]] ++
        [ maximum samples ]
    doOne = do
        dp <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)
        dq <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)

        let sp = leftToRight bits w dp
            wlp = possibleWindowLengths bits w sp
            hardBitsp = countHardBits bits w wlp
            possp = countPossibilities bits w wlp

            sq = leftToRight bits w dq
            wlq = possibleWindowLengths bits w sq
            hardBitsq = countHardBits bits w wlq
            possq = countPossibilities bits w wlq

            totalBits = fromIntegral (2*bits)
            hardBitsFrac = fromIntegral (hardBitsp + hardBitsq) / totalBits

            realBitsFrac = ( totalBits - logBase 2.0 (fromIntegral possp) - logBase 2.0 (fromIntegral possq) ) / totalBits

            rtlBitsFrac = ( fromIntegral $ rightToLeftBits bits w dp + rightToLeftBits bits w dq
                          ) / totalBits

            gmpBitsFrac = ( fromIntegral $ leftToRightGMPBits bits w dp + leftToRightGMPBits bits w dq
                          ) / totalBits

        return (rtlBitsFrac, hardBitsFrac, realBitsFrac, gmpBitsFrac)

hardSoftScatter :: [BitLength] -> [Width] -> Int -> IO ()
hardSoftScatter bitsRange wRange iter = do
    putStrLn $ intercalate "\t" $ words "bits width hard soft"
    forM_ bitsRange $ \bits -> forM_ wRange $ \w -> do
        pairs <- replicateM iter (doOne bits w)
        forM_ pairs $ \(h,s) ->
            putStrLn $ intercalate "\t" (map show [bits, w] ++ map show [h,s])
  where
    doOne bits w = do
        dp <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)
        dq <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)

        let sp = leftToRight bits w dp
            wlp = possibleWindowLengths bits w sp
            hardBitsp = countHardBits bits w wlp
            possp = countPossibilities bits w wlp

            sq = leftToRight bits w dq
            wlq = possibleWindowLengths bits w sq
            hardBitsq = countHardBits bits w wlq
            possq = countPossibilities bits w wlq

            totalBits = fromIntegral (2*bits) :: Double
            hardBitsFrac = fromIntegral (hardBitsp + hardBitsq) / totalBits

            realBitsFrac = ( totalBits - logBase 2.0 (fromIntegral possp) - logBase 2.0 (fromIntegral possq) ) / totalBits

            rtlBitsFrac = ( fromIntegral $ rightToLeftBits bits w dp + rightToLeftBits bits w dq
                          ) / totalBits

            gmpBitsFrac = ( fromIntegral $ leftToRightGMPBits bits w dp + leftToRightGMPBits bits w dq
                          ) / totalBits

        return (hardBitsFrac, realBitsFrac)

leakPreinput ::  Width -> IO ()
leakPreinput w = do
    line <- getLine
    let line' = fixFirst line
    unless (all (`elem` "SM_") line') $ do
        hPutStrLn stderr "Error: Input is not a sequence of 'S', '_' or 'M'"
        exitFailure
    let bits = length line'
        sp = V.fromList (map toStep line')
        wlp = possibleWindowLengths bits w sp
            --hardBitsp = countHardBits bits w wlp
        possp = countPossibilities bits w wlp
    print possp
  where
    toStep 'M' = M
    toStep 'S' = S
    toStep '_' = S

    -- The input may explicitly say that the first character is a 1.
    -- This code knows that, but still expects there to be a S or M.
    -- Recover whether this is a S or M from looking at the following characters.
    fixFirst ('1':xs) = if 'M' `elem` take (w - 1) xs then 'S':xs else 'M':xs
    fixFirst xs = xs

showRange :: [Int] -> String
showRange = intercalate "," . map show

parseRange :: String -> Either String [Int]
parseRange s = case splitOn "," s of
    [n1,n2,"..",n3] -> enumFromThenTo <$> readEither n1 <*> readEither n2 <*> readEither n3
    [n1,"..",n3] -> enumFromTo <$> readEither n1 <*> readEither n3
    ns -> mapM readEither ns

max_stars_repo_name: nomeata/slidingright

max_issues_repo_name: nomeata/slidingright

max_forks_repo_name: nomeata/slidingright


text:
module Main where

{-# LANGUAGE ScopedTypeVariables #-}

import Data.Maybe

import Graphics.UI.WX hiding (Event)
import Reactive.Banana
import Reactive.Banana.WX

import AI.HNN.FF.Network
import Numeric.LinearAlgebra

import Data.List.Split
import Data.List (replicate)

main :: IO ()
main = start gui

gui :: IO ()
gui = do
    f       <- frame [ text := "Neural Network" ]

    lEpocs          <- staticText f [ text := "Epocs"]
    lMiddleLayers   <- staticText f [ text := "Middle & output layers"]
    lLearningRate   <- staticText f [ text := "Learning rate"]

    iEpocs          <- textCtrlRich f [ text := "1000" ]
    iMiddleLayers   <- textCtrlRich f [ text := "2,2,1"]
    iLearningRate   <- textCtrlRich f [ text := "0.8"]

    bStart  <- button f [ text := "Start" ]
    bQuit   <- button f [ text := "Quit" ]

    set bStart [ on command := parseInputAndTrain f iEpocs iMiddleLayers iLearningRate]
    set bQuit  [ on command := close f]

    set f [ layout := column 25 [ widget lEpocs, widget iEpocs,
                                  widget lMiddleLayers, widget iMiddleLayers,
                                  widget lLearningRate, widget iLearningRate,
                                  widget bStart, widget bQuit ] ]

parseInputAndTrain :: (Paint w, Textual w1, Textual w2, Textual w3) => w -> w1 -> w2 -> w3 -> IO ()
parseInputAndTrain f iEpocs iMiddleLayers iLearningRate = do
    epocsString         <- get iEpocs text
    middleLayersString  <- get iMiddleLayers text
    learningRateString  <- get iLearningRate text
    let epocs = read epocsString :: Int
    let middleLayersSplited = splitOn "," middleLayersString
    let middleLayer = map (\s -> read s :: Int) middleLayersSplited
    let learningRate = read learningRateString :: Double
    trainNeuralNetwork epocs 2 middleLayer 1 learningRate

trainNeuralNetwork :: Int -> Int -> [Int] -> Int -> Double -> IO()
trainNeuralNetwork epocs firstLayer middleLayers outputLayer learningRate =  do

    xys <- readFile "xy.txt"
    zs <- readFile "z.txt"

    let list1 = parse xys
    let list2 = parse zs
    let list11 = map fromList list1
    let list22 = map fromList list2
    let samples2 = zip list11 list22

    putStrLn "1------------------"
    n <- createNetwork firstLayer middleLayers outputLayer :: IO (Network Double)
    mapM_ (putStrLn . show) samples2
    mapM_ (putStrLn . show . output n tanh . fst) samples2
    putStrLn "2------------------"
    let n' = trainNTimes epocs learningRate tanh tanh' n samples2
    mapM_ (putStrLn . show . output n' tanh . fst) samples2

parse :: String -> [[Double]]
parse a = map (map read . words) (lines a)

oneList :: [a] -> [b] -> [(a, b)]
oneList []     _      = []
oneList (x:xs) (y:ys) = (x, y) : oneList xs ys

samples :: Samples Double
samples = [ (fromList [0, 0], fromList [0])
          , (fromList [0, 1], fromList [1])
          , (fromList [1, 0], fromList [1])
          , (fromList [1, 1], fromList [0])
          ]

ext: hs

max_stars_repo_name: tomymehdi/ai-haskell

max_issues_repo_name: tomymehdi/ai-haskell

max_forks_repo_name: tomymehdi/ai-haskell


text:
module NeuralNet.Net (
  NeuralNetDefinition,
  LayerDefinition (LayerDefinition),
  WeightsStream,
  NeuralNet,
  initNN,
  buildNNFromList,
  nnLayers,
  nnNumInputs,
  nnForward,
  nnForwardSet,
  isExampleSetCompatibleWithNN,
  isExampleSetCompatibleWithNNDef,
  updateNNLayers,
  createLogRegDefinition
) where

import NeuralNet.Matrix
import NeuralNet.Activation
import NeuralNet.Layer
import NeuralNet.Example
import Numeric.LinearAlgebra

type NumNeurons = Int
type NumInputs = Int
type NeuralNetDefinition = (NumInputs, [LayerDefinition])

type WeightsStream = [Double]

data LayerDefinition = LayerDefinition Activation NumNeurons
  deriving (Show, Eq)

data NeuralNet = NeuralNet [NeuronLayer]
  deriving (Show, Eq)

createLogRegDefinition :: NumInputs -> Activation -> NeuralNetDefinition
createLogRegDefinition n a = (n, [LayerDefinition a 1])

initNN :: WeightsStream -> NeuralNetDefinition -> NeuralNet
initNN _ (_, []) = error "No layer definitions provided"
initNN g (numInputs, layerDefs)
  | numInputs > 0 = NeuralNet (fst (initNeuronLayers g numInputs layerDefs))
  | otherwise     = error "Need positive num inputs"

initNeuronLayers :: WeightsStream -> Int -> [LayerDefinition] -> ([NeuronLayer], WeightsStream)
initNeuronLayers weightsStream _ [] = ([], weightsStream)
initNeuronLayers weightsStream numInputs (LayerDefinition a numNeurons : ds) = (NeuronLayer a w b : ls, ws2)
  where
    numSeeds = numInputs * numNeurons
    (weights, ws) = splitAt numSeeds weightsStream
    w = scale (sqrt (2.0 / fromIntegral numInputs)) ((numNeurons >< numInputs) weights)
    -- * sqrt (2 / numInputs) should be parameterisable
    b = col (replicate numNeurons 0)
    (ls, ws2) = initNeuronLayers ws numNeurons ds

buildNNFromList :: NeuralNetDefinition -> [Double] -> NeuralNet
buildNNFromList def@(numInputs, layerDefs) list
  | listSize /= defSize = error ("Wrong list size. Expected " ++ show defSize ++ " got " ++ show listSize)
  | otherwise           = NeuralNet (buildLayersFromList numInputs layerDefs list)
    where
      defSize = definitionToNNSize def
      listSize = length list

-- TODO: nnToList and include back and forward of buildNNFromList to ensure works

definitionToNNSize :: NeuralNetDefinition -> Int
definitionToNNSize (numInputs, layerDefs) = defToSizeStep numInputs layerDefs
  where
    defToSizeStep :: Int -> [LayerDefinition] -> Int
    defToSizeStep _ [] = 0
    defToSizeStep stepInputs (LayerDefinition _ numNeurons:ds) = (stepInputs * numNeurons) + numNeurons + defToSizeStep numNeurons ds

buildLayersFromList :: Int -> [LayerDefinition] -> [Double] -> [NeuronLayer]
buildLayersFromList _ [] _ = []
buildLayersFromList numInputs (LayerDefinition a numNeurons:ds) xs = layer : (buildLayersFromList numNeurons ds restXs)
  where
    numWeights = numInputs * numNeurons
    (layerXs, afterLayerXs) = splitAt numWeights xs
    (bs, restXs) = splitAt numNeurons afterLayerXs
    layer = NeuronLayer a ((numNeurons><numInputs) layerXs) (col bs)

nnLayers :: NeuralNet -> [NeuronLayer]
nnLayers (NeuralNet layers) = layers

nnForward :: NeuralNet -> [Double] -> [Double]
nnForward nn inputs = concat (toLists (forwardPropA (last setResult)))
  where
    set = createExampleSet [(inputs, error "Undefined")]
    setResult = nnForwardSet nn set

nnForwardSet :: NeuralNet -> ExampleSet -> [ForwardPropStep]
nnForwardSet nn examples = reverse foldResult
  where
    x = exampleSetX examples :: Matrix Double
    input = createInputForwardPropStep x
    foldResult = foldl (\i l -> layerForwardSet l (forwardPropA (head i)) : i) [input] (nnLayers nn)

nnNumInputs :: NeuralNet -> Int
nnNumInputs = layerNumInputs . head . nnLayers

isExampleSetCompatibleWithNNDef :: ExampleSet -> NeuralNetDefinition -> Bool
isExampleSetCompatibleWithNNDef examples def = isExampleSetCompatibleWithNN examples emptyNN
  where
    nums = replicate (definitionToNNSize def) 0
    emptyNN = buildNNFromList def nums

isExampleSetCompatibleWithNN :: ExampleSet -> NeuralNet -> Bool
isExampleSetCompatibleWithNN examples nn = nnNumInputs nn == exampleSetN examples

updateNNLayers :: NeuralNet -> [NeuronLayer] -> NeuralNet
updateNNLayers nn newLayers
  | layersShapes layers /= layersShapes newLayers = error "Must provide same shaped layers"
  | otherwise                                     = NeuralNet newLayers
    where
      layersShapes :: [NeuronLayer] -> [((Int, Int), (Int, Int))]
      layersShapes sLayers = map (\l -> (dims (layerW l), dims (layerB l))) sLayers
      layers = nnLayers nn

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net


text:
-- File created: 2009-07-15 21:15:44

{-# LANGUAGE Rank2Types #-}

module Haschoo.Evaluator.Standard.Numeric
   (procedures, isNumeric, isExact, numEq, asInt) where

import Control.Arrow     ((&&&), (***))
import Control.Monad     (ap, foldM)
import Data.Array.IArray (elems)
import Data.Array.MArray (getElems)
import Data.Char         (intToDigit)
import Data.Complex      ( Complex((:+)), mkPolar
                         , imagPart, realPart, phase, magnitude)
import Data.Ratio        (numerator, denominator, approxRational)
import Data.Function     (on)
import Numeric           (showIntAtBase, showSigned)

import qualified Haschoo.Parser as Parser
import           Haschoo.Types           (ScmValue(..), toScmMString)
import           Haschoo.Utils           (ErrOr, allM, ($<), (.:))
import           Haschoo.Evaluator.Utils (tooFewArgs, tooManyArgs, notInt)

procedures :: [(String, ScmValue)]
procedures = map (\(a,b) -> (a, ScmFunc a (return . b)))
   [ ("number?",   fmap ScmBool . scmIsNumber)
   , ("complex?",  fmap ScmBool . scmIsNumber)
   , ("real?",     fmap ScmBool . scmIsReal)
   , ("rational?", fmap ScmBool . scmIsRational)
   , ("integer?",  fmap ScmBool . scmIsInteger)

   , "exact?"   $< id &&& fmap  ScmBool        .: scmIsExact
   , "inexact?" $< id &&& fmap (ScmBool . not) .: scmIsExact

   , ("=", fmap ScmBool . scmNumEq)
   , "<"  $< id &&& fmap ScmBool .: scmCompare (<)
   , ">"  $< id &&& fmap ScmBool .: scmCompare (>)
   , "<=" $< id &&& fmap ScmBool .: scmCompare (<=)
   , ">=" $< id &&& fmap ScmBool .: scmCompare (>=)

   , ("zero?",     fmap ScmBool . scmIsZero)
   , ("positive?", fmap ScmBool . scmIsPos)
   , ("negative?", fmap ScmBool . scmIsNeg)
   , "even?" $< id &&& fmap  ScmBool        .: scmIsEven
   , "odd?"  $< id &&& fmap (ScmBool . not) .: scmIsEven

   , ("max", scmMax)
   , ("min", scmMin)

   , ("+", scmPlus)
   , ("-", scmMinus)
   , ("*", scmMul)
   , ("/", scmDiv)

   , ("abs", scmAbs)

   , ("quotient",  scmQuot)
   , ("remainder", scmRem)
   , ("modulo",    scmMod)

   , ("gcd", scmGcd)
   , ("lcm", scmLcm)

   , ("numerator",   scmNumerator)
   , ("denominator", scmDenominator)

   , ("floor",    scmFloor)
   , ("ceiling",  scmCeil)
   , ("truncate", scmTrunc)
   , ("round",    scmRound)

   , ("rationalize", scmRationalize)

   , ("exp",  scmExp)
   , ("log",  scmLog)
   , ("sin",  scmSin)
   , ("cos",  scmCos)
   , ("tan",  scmTan)
   , ("asin", scmAsin)
   , ("acos", scmAcos)
   , ("atan", scmAtan)

   , ("sqrt", scmSqrt)
   , ("expt", scmExpt)

   , ("make-rectangular", scmMakeRectangular)
   , ("make-polar",       scmMakePolar)
   , ("real-part",        scmRealPart)
   , ("imag-part",        scmImagPart)
   , ("magnitude",        scmNorm)
   , ("angle",            scmAngle)

   , ("exact->inexact", scmToInexact)
   , ("inexact->exact", scmToExact) ]

   ++ map (\(a,b) -> (a, ScmFunc a b))
   [ ("number->string", scmToString)
   , ("string->number", scmToNumber) ]

---- Predicates

scmIsNumber, scmIsReal, scmIsRational, scmIsInteger :: [ScmValue] -> ErrOr Bool
scmIsNumber [x] = Right $ isNumeric x
scmIsNumber []  = tooFewArgs  "number?"
scmIsNumber _   = tooManyArgs "number?"

scmIsReal [x] = Right $ isReal x
scmIsReal []  = tooFewArgs  "real?"
scmIsReal _   = tooManyArgs "real?"

scmIsRational [ScmRat _] = Right True
scmIsRational [x]        = Right $ isInteger x
scmIsRational []         = tooFewArgs "rational?"
scmIsRational _          = tooManyArgs "rational?"

scmIsInteger [x] = Right $ isInteger x
scmIsInteger []  = tooFewArgs  "integer?"
scmIsInteger _   = tooManyArgs "integer?"

scmIsExact :: String -> [ScmValue] -> ErrOr Bool
scmIsExact _ [x] | isNumeric x = Right $ isExact x
scmIsExact s [_]               = notNum s
scmIsExact s []                = tooFewArgs s
scmIsExact s _                 = tooManyArgs s

isExact :: ScmValue -> Bool
isExact (ScmInt _) = True
isExact (ScmRat _) = True
isExact _          = False

---- Comparison

scmNumEq :: [ScmValue] -> ErrOr Bool
scmNumEq [] = tooFewArgs "="
scmNumEq xs = allM f . (zip`ap`tail) $ xs
 where
   f (a,b) = if isNumeric a && isNumeric b
                then Right $ numEq a b
                else notNum "="

numEq :: ScmValue -> ScmValue -> Bool
numEq x y =
   case pairScmComplex x y of
        Right (ScmInt     a, ScmInt     b) -> a == b
        Right (ScmRat     a, ScmRat     b) -> a == b
        Right (ScmReal    a, ScmReal    b) -> a == b
        Right (ScmComplex a, ScmComplex b) -> a == b
        _ -> error "numEq :: the impossible happened"

scmCompare :: (forall a. Real a => a -> a -> Bool)
           -> String -> [ScmValue] -> ErrOr Bool
scmCompare p s xs@(_:_:_) = allM f . (zip`ap`tail) $ xs
 where
   f (a,b) = case liftScmRealA2 p a b of
                  Left _ -> notReal s
                  x      -> x

scmCompare _ s _ = tooFewArgs s

scmIsZero :: [ScmValue] -> ErrOr Bool
scmIsZero [ScmInt     0]    = Right True
scmIsZero [ScmRat     0]    = Right True
scmIsZero [ScmReal    0]    = Right True
scmIsZero [ScmComplex 0]    = Right True
scmIsZero [x] | isNumeric x = Right False
scmIsZero [_]               = notNum "zero?"
scmIsZero _                 = tooManyArgs "zero?"

scmIsPos, scmIsNeg :: [ScmValue] -> ErrOr Bool
scmIsPos [x] = scmCompare (>) "positive?" [x, ScmInt 0]
scmIsPos _   = tooManyArgs "positive?"

scmIsNeg [x] = scmCompare (<) "negative?" [x, ScmInt 0]
scmIsNeg _   = tooManyArgs "negative?"

scmIsEven :: String -> [ScmValue] -> ErrOr Bool
scmIsEven s [x] = if isInteger x
                     then (scmIsZero.return) =<< scmMod [x, ScmInt 2]
                     else notInt s
scmIsEven s _   = tooManyArgs s

scmMax, scmMin :: [ScmValue] -> ErrOr ScmValue
scmMax = scmMinMax max "max"
scmMin = scmMinMax min "min"

scmMinMax :: (forall a. Real a => a -> a -> a) -> String
          -> [ScmValue] -> ErrOr ScmValue
scmMinMax _ s []     = tooFewArgs s
scmMinMax f s (x:xs) = if isNumeric x then foldM go x xs else notNum s
 where
   go n a = if isNumeric a
               then case liftScmReal2 f n a of
                         Left _ -> notReal s
                         m      -> m
               else notNum s

---- +-*/

scmPlus, scmMinus, scmMul, scmDiv :: [ScmValue] -> ErrOr ScmValue
scmPlus [] = Right $ ScmInt 0
scmPlus xs = foldM go (ScmInt 0) xs
 where
   go n x = if isNumeric x then liftScmNum2 (+) n x else notNum "+"

scmMinus []                        = tooFewArgs "-"
scmMinus (x:_) | not (isNumeric x) = notNum "-"
scmMinus [x]                       = liftScmNum negate x
scmMinus (x:xs)                    = foldM go x xs
 where
   go n a = if isNumeric a then liftScmNum2 (-) n a else notNum "-"

scmMul [] = Right $ ScmInt 1
scmMul xs = foldM go (ScmInt 1) xs
 where
   go n x = if isNumeric x then liftScmNum2 (*) n x else notNum "*"

scmDiv []     = tooFewArgs "/"
scmDiv (x:xs) =
   unint x >>= case xs of
                    []  -> liftScmFrac recip
                    _:_ -> flip (foldM go) xs
 where
   go n a = unint a >>= liftScmFrac2 (/) n

   unint (ScmInt n) = Right $ ScmRat (fromInteger n)
   unint n          = if isNumeric x then Right n else notNum "/"

---- abs

scmAbs :: [ScmValue] -> ErrOr ScmValue
scmAbs [x] =
   if isReal x
      then liftScmNum abs x
      else notReal "abs"
scmAbs [] = tooFewArgs  "abs"
scmAbs _  = tooManyArgs "abs"

---- quot rem mod

scmQuot, scmRem, scmMod :: [ScmValue] -> ErrOr ScmValue
scmQuot = scmQuotRemMod quot "quotient"
scmRem  = scmQuotRemMod rem "remainder"
scmMod  = scmQuotRemMod mod "modulo"

scmQuotRemMod :: (Integer -> Integer -> Integer) -> String
              -> [ScmValue] -> ErrOr ScmValue
scmQuotRemMod f s [x,y] = do
   (e1,a) <- asInt s x
   (e2,b) <- asInt s y
   return . (if e1 && e2 then ScmInt else ScmReal . fromInteger) $ f a b
scmQuotRemMod _ s (_:_:_) = tooManyArgs s
scmQuotRemMod _ s _       = tooFewArgs  s

---- gcd lcm

scmGcd, scmLcm :: [ScmValue] -> ErrOr ScmValue
scmGcd = scmGcdLcm (\a b -> if a == 0 && b == 0 then 0 else gcd a b) 0 "gcd"
scmLcm = scmGcdLcm lcm 1 "lcm"

scmGcdLcm :: (Integer -> Integer -> Integer) -> Integer -> String
          -> [ScmValue] -> ErrOr ScmValue
scmGcdLcm f i s = fmap fin . foldM go (True,i)
 where
   go  (exact, n) = fmap ((exact &&) *** f n) . asInt s
   fin (exact, n) = if exact then ScmInt n else ScmReal (fromInteger n)

---- numerator denominator

scmNumerator, scmDenominator :: [ScmValue] -> ErrOr ScmValue
scmNumerator   = scmNumerDenom numerator   "numerator"
scmDenominator = scmNumerDenom denominator "denominator"

scmNumerDenom :: (Rational -> Integer) -> String
              -> [ScmValue] -> ErrOr ScmValue
scmNumerDenom f s [n] =
   case n of
        ScmInt      x       -> Right . ScmInt  $ f (fromInteger x)
        ScmRat      x       -> Right . ScmInt  $ f x
        ScmReal     x       -> Right . ScmReal . fromInteger $ f (realToFrac x)
        ScmComplex (x :+ 0) -> Right . ScmReal . fromInteger $ f (realToFrac x)
        _                   -> notRat s

scmNumerDenom _ s [] = tooFewArgs s
scmNumerDenom _ s _  = tooManyArgs s

---- floor ceil trunc round

scmFloor, scmCeil, scmTrunc, scmRound :: [ScmValue] -> ErrOr ScmValue
scmFloor = scmGenericRound floor    "floor"
scmCeil  = scmGenericRound ceiling  "ceiling"
scmTrunc = scmGenericRound truncate "truncate"
scmRound = scmGenericRound round    "round"

scmGenericRound :: (forall a. RealFrac a => a -> Integer) -> String
                -> [ScmValue] -> ErrOr ScmValue
scmGenericRound _ _ [   ScmInt     x]  = Right . ScmInt                $ x
scmGenericRound f _ [   ScmRat     x]  = Right . ScmInt                $ f x
scmGenericRound f _ [   ScmReal    x]  = Right . ScmReal . fromInteger $ f x
scmGenericRound _ s [n@(ScmComplex _)] =
   case asInt s n of
        Right (_,x) -> Right . ScmInt $ x
        Left _      -> notReal s
scmGenericRound _ s [_]            = notNum s
scmGenericRound _ s []             = tooFewArgs s
scmGenericRound _ s _              = tooManyArgs s

---- rationalize

scmRationalize :: [ScmValue] -> ErrOr ScmValue
scmRationalize (_:_:_:_) = tooManyArgs "rationalize"
scmRationalize [x,y]     =
   case liftScmRealFracB2 approxRational x y of
        Left _           -> notReal "rationalize"
        Right (True , n) -> Right (ScmRat  n)
        Right (False, n) -> Right (ScmReal (fromRational n))

scmRationalize _         = tooFewArgs "rationalize"

---- exp log sin cos tan asin acos atan sqrt

scmExp, scmLog,
 scmSin,  scmCos,  scmTan,
 scmAsin, scmAcos, scmAtan,
 scmSqrt :: [ScmValue] -> ErrOr ScmValue

scmExp  = scmComplex1 exp  "exp"
scmLog  = scmComplex1 log  "log"
scmSin  = scmComplex1 sin  "sin"
scmCos  = scmComplex1 cos  "cos"
scmTan  = scmComplex1 tan  "tan"
scmAsin = scmComplex1 asin "asin"
scmAcos = scmComplex1 acos "acos"

scmAtan [y,x] =
   case pairScmReal x y of
        Left  _                      -> notReal "atan"
        Right (ScmReal a, ScmReal b) -> Right (ScmReal $ phase $ a :+ b)
        Right (ScmInt a, ScmInt b)   ->
           Right (ScmReal $ phase $ fromInteger a :+ fromInteger b)
        Right (ScmRat a, ScmRat b)   ->
           Right (ScmReal $ phase $ fromRational a :+ fromRational b)
        Right _                      -> error "scmAtan :: internal error"

scmAtan xs    = scmComplex1 atan "atan" xs

scmSqrt = scmComplex1 sqrt "sqrt"

scmComplex1 :: (Complex Double -> Complex Double) -> String
            -> [ScmValue] -> ErrOr ScmValue
scmComplex1 f _ [ScmInt     x] = Right . fromComplex $ f (fromInteger  x :+ 0)
scmComplex1 f _ [ScmRat     x] = Right . fromComplex $ f (fromRational x :+ 0)
scmComplex1 f _ [ScmReal    x] = Right . fromComplex $ f (x :+ 0)
scmComplex1 f _ [ScmComplex x] = Right . fromComplex $ f x
scmComplex1 _ s [_]            = notNum s
scmComplex1 _ s []             = tooFewArgs s
scmComplex1 _ s _              = tooManyArgs s

---- expt

scmExpt :: [ScmValue] -> ErrOr ScmValue
scmExpt [x,y] =
   case pairScmComplex x y of
        Left _                              -> notNum "expt"
        Right (ScmInt     a, ScmInt      b) ->
           if b < 0
              then Right . ScmReal $ fromInteger a ^^ b
              else Right . ScmInt  $ a ^ b
        Right (ScmReal    a, ScmReal     b) -> Right . ScmReal    $ a ** b
        Right (ScmComplex a, ScmComplex  b) -> Right . ScmComplex $ a ** b
        Right (ScmRat     a, ScmRat      b) ->
           Right . ScmReal $ fromRational a ** fromRational b
        Right _                             -> error "expt :: internal error"

scmExpt (_:_:_) = tooManyArgs "expt"
scmExpt _       = tooFewArgs  "expt"

---- make-rectangular make-polar real-part imag-part magnitude angle

scmMakeRectangular, scmMakePolar :: [ScmValue] -> ErrOr ScmValue
scmMakeRectangular = scmMakeComplex (:+)    "make-rectangular"
scmMakePolar       = scmMakeComplex mkPolar "make-polar"

scmMakeComplex :: (Double -> Double -> Complex Double) -> String
               -> [ScmValue] -> ErrOr ScmValue
scmMakeComplex f s [x,y] =
   case pairScmReal x y of
        Left _                       -> notReal s
        Right (ScmReal a, ScmReal b) -> Right . ScmComplex $ f a b
        Right (ScmRat  a, ScmRat  b) ->
           Right . ScmComplex $ (f `on` fromRational) a b
        Right (ScmInt  a, ScmInt  b) ->
           Right . ScmComplex $ (f `on` fromInteger) a b
        Right _                      ->
           error $ s ++ " :: internal error"

scmMakeComplex _ s (_:_:_) = tooManyArgs s
scmMakeComplex _ s _       = tooFewArgs  s

scmRealPart, scmImagPart, scmNorm, scmAngle :: [ScmValue] -> ErrOr ScmValue
scmRealPart = scmComplexPart realPart  Right                  "real-part"
scmImagPart = scmComplexPart imagPart  (liftScmNum $ const 0) "imag-part"
scmNorm     = scmComplexPart magnitude (liftScmNum abs)       "magnitude"
scmAngle    = scmComplexPart phase     (Right . angle)        "angle"
 where
   angle (ScmReal n) = f ScmReal n
   angle (ScmRat  n) = f ScmRat  n
   angle (ScmInt  n) = f ScmInt  n
   angle _           = error "angle :: the impossible happened"

   f c n = if n >= 0 then c 0 else ScmReal pi

scmComplexPart :: (Complex Double -> Double)
               -> (ScmValue -> ErrOr ScmValue)
               -> String
               -> [ScmValue] -> ErrOr ScmValue
scmComplexPart f _ _ [ScmComplex x]    = Right $ ScmReal (f x)
scmComplexPart _ g _ [x] | isNumeric x = g x
scmComplexPart _ _ s [_]               = notNum s
scmComplexPart _ _ s []                = tooFewArgs s
scmComplexPart _ _ s _                 = tooManyArgs s

---- exact->inexact inexact->exact

scmToInexact, scmToExact :: [ScmValue] -> ErrOr ScmValue
scmToInexact [ScmInt x]         = Right $ ScmReal (fromInteger  x)
scmToInexact [ScmRat x]         = Right $ ScmReal (fromRational x)
scmToInexact [x] | isNumeric x  = Right x
scmToInexact [_]                = notNum      "exact->inexact"
scmToInexact []                 = tooFewArgs  "exact->inexact"
scmToInexact _                  = tooManyArgs "exact->inexact"

scmToExact [ScmReal    x]     =
   let r = toRational x
    in Right $ if denominator r == 1
                  then ScmInt (numerator r)
                  else ScmRat r
scmToExact [ScmComplex x]     =
   case imagPart x of
        0 -> Right $ ScmReal (realPart x)
        _ -> fail $ "inexact->exact :: implementation restriction: " ++
                    "can't exactify complex numbers"
scmToExact [x] | isNumeric x  = Right x
scmToExact [_]                = notNum      "exact->inexact"
scmToExact []                 = tooFewArgs  "exact->inexact"
scmToExact _                  = tooManyArgs "exact->inexact"

---- number->string string->number

scmToString, scmToNumber :: [ScmValue] -> IO (ErrOr ScmValue)

-- TODO: scmShow should use this or vice versa
scmToString (x:xs) | isNumeric x =
   case xs of
        []             -> g x 10
        [ScmInt radix] ->
           if radix `elem` [2,8,10,16]
              then g x (fromInteger radix)
              else return. fail $
                      "number->string :: invalid radix " ++ show radix

        [_] -> return$ notInt      "number->string"
        _   -> return$ tooManyArgs "number->string"
 where
   g n radix = case f n radix of
                    Left  err -> return$ fail err
                    Right s   -> fmap Right $ toScmMString s

   f (ScmInt i) radix = Right $ showInt radix i
   f (ScmRat r) radix = Right $
      if denominator r == 1
         then showInt radix (numerator r)
         else concat [ showInt radix (numerator r)
                     , "/"
                     , showInt radix (denominator r) ]

   f (ScmReal r) radix =
      if radix == 10
         then Right $
                 case () of
                      _| isNaN      r -> "+nan.#"
                       | isInfinite r -> (if r < 0 then '-' else '+') : "inf.#"

                       -- R5RS 6.2.6 specifies that "the result contains a
                       -- decimal point" if inexact and a decimal-point
                       -- containing result can make it work
                       | otherwise    -> show r
         else fail "number->string :: nondecimal radix for inexact real"

   f (ScmComplex (a :+ b)) radix =
      -- As in the ScmReal case
      if radix == 10
         then Right . concat $ let bs = either undefined id $ f (ScmReal b) 10
                                in [ either undefined id $ f (ScmReal a) 10
                                   , if head bs `elem` "-+" then [] else "+"
                                   , bs
                                   , "i" ]

         else fail "number->string :: nondecimal radix for inexact complex"

   f _ _ = error "number->string :: the impossible happened"

   showInt radix i = showSigned (showIntAtBase radix intToDigit) 0 i ""

scmToString [] = return$ tooFewArgs "number->string"
scmToString _  = return$ notNum     "number->string"

scmToNumber (ScmString  s : xs) = return$ toNumHelper xs     (elems s)
scmToNumber (ScmMString s : xs) = fmap   (toNumHelper xs) (getElems s)
scmToNumber [] = return$ tooFewArgs "string->number"
scmToNumber _  = return$ fail "Nonstring argument to string->number"

toNumHelper :: [ScmValue] -> String -> ErrOr ScmValue
toNumHelper xs s =
   case xs of
        []             -> Right $ f s 10
        [ScmInt radix] ->
           if radix `elem` [2,8,10,16]
              then Right $ f s (fromInteger radix)
              else fail $ "string->number :: invalid radix " ++ show radix

        [_]            -> notInt       "string->number"
        _              -> tooManyArgs  "string->number"
 where
   f str radix =
       case Parser.runParser (Parser.number radix) "" str of
            Right n | isNumeric n -> n
            _                     -> ScmBool False

-------------

notNum, notReal, notRat :: String -> ErrOr a
notNum  = fail . ("Nonnumeric argument to " ++)
notReal = fail . ("Nonreal argument to " ++)
notRat  = fail . ("Nonrational argument to " ++)

isNumeric, isInteger, isReal :: ScmValue -> Bool
isNumeric (ScmInt     _) = True
isNumeric (ScmRat     _) = True
isNumeric (ScmReal    _) = True
isNumeric (ScmComplex _) = True
isNumeric _              = False

isInteger (ScmInt     _)      = True
isInteger (ScmRat     x)      = denominator x == 1
isInteger (ScmReal    x)      = x == fromInteger (round x)
isInteger (ScmComplex (a:+b)) = b == 0 && a == fromInteger (round a)
isInteger _                   = False

isReal (ScmComplex (_ :+ b)) = b == 0
isReal x                     = isNumeric x

asInt :: String -> ScmValue -> ErrOr (Bool, Integer)
asInt s x | not (isInteger x) = notInt s
asInt _ (ScmInt     x)        = Right (True,  x)
asInt _ (ScmRat     x)        = Right (True,  numerator x)
asInt _ (ScmReal    x)        = Right (False, round x)
asInt _ (ScmComplex (x:+0))   = Right (False, round x)
asInt _ _                     = error "This can't happen"

fromComplex :: Complex Double -> ScmValue
fromComplex (x :+ 0) = ScmReal x
fromComplex x        = ScmComplex x

-- Lift ScmValue's numeric types to a common type
--
-- Versions without 'A' at the end also return results in the correct
-- constructor of the two given. Versions with 'B' just return whether the
-- result should be exact or not.
--
-- Ugly and verbose... too lazy to metaize these {{{

--- Num
liftScmNum :: (forall a. Num a => a -> a) -> ScmValue -> ErrOr ScmValue
liftScmNum f (ScmInt     a) = Right . ScmInt     $ f a
liftScmNum f (ScmRat     a) = Right . ScmRat     $ f a
liftScmNum f (ScmReal    a) = Right . ScmReal    $ f a
liftScmNum f (ScmComplex a) = Right . ScmComplex $ f a
liftScmNum _ _ = fail "liftScmNum :: internal error"

liftScmNum2 :: (forall a. Num a => a -> a -> a)
            -> (ScmValue -> ScmValue -> ErrOr ScmValue)
liftScmNum2 f (ScmInt     a) (ScmInt     b) = Right . ScmInt    $ f a b
liftScmNum2 f (ScmRat     a) (ScmRat     b) = Right . ScmRat    $ f a b
liftScmNum2 f (ScmReal    a) (ScmReal    b) = Right . ScmReal   $ f a b
liftScmNum2 f (ScmComplex a) (ScmComplex b) = Right . ScmComplex$ f a b

-- Int+{Rat,Real,Complex}
liftScmNum2 f (ScmInt     a) (ScmRat     b) = Right . ScmRat    $ f (fInt a) b
liftScmNum2 f (ScmRat     a) (ScmInt     b) = Right . ScmRat    $ f a (fInt b)
liftScmNum2 f (ScmInt     a) (ScmReal    b) = Right . ScmReal   $ f (fInt a) b
liftScmNum2 f (ScmReal    a) (ScmInt     b) = Right . ScmReal   $ f a (fInt b)
liftScmNum2 f (ScmInt     a) (ScmComplex b) = Right . ScmComplex$ f (fInt a) b
liftScmNum2 f (ScmComplex a) (ScmInt     b) = Right . ScmComplex$ f a (fInt b)

-- Rat+{Real,Complex}
liftScmNum2 f (ScmRat     a) (ScmReal    b) = Right . ScmReal   $ f (fRat a) b
liftScmNum2 f (ScmReal    a) (ScmRat     b) = Right . ScmReal   $ f a (fRat b)
liftScmNum2 f (ScmRat     a) (ScmComplex b) = Right . ScmComplex$ f (fRat a) b
liftScmNum2 f (ScmComplex a) (ScmRat     b) = Right . ScmComplex$ f a (fRat b)

-- Real+Complex
liftScmNum2 f (ScmReal    a) (ScmComplex b) = Right . ScmComplex$ f (a :+ 0) b
liftScmNum2 f (ScmComplex a) (ScmReal    b) = Right . ScmComplex$ f a (b :+ 0)

liftScmNum2 _ _ _ = fail "liftScmNum2 :: internal error"

--- Frac
liftScmFrac :: (forall a. Fractional a => a -> a) -> ScmValue -> ErrOr ScmValue
liftScmFrac f (ScmRat     a) = Right . ScmRat     $ f a
liftScmFrac f (ScmReal    a) = Right . ScmReal    $ f a
liftScmFrac f (ScmComplex a) = Right . ScmComplex $ f a
liftScmFrac _ _ = fail "liftScmFrac :: internal error"

liftScmFrac2 :: (forall a. Fractional a => a -> a -> a)
             -> (ScmValue -> ScmValue -> ErrOr ScmValue)
liftScmFrac2 f (ScmRat     a) (ScmRat     b) = Right . ScmRat    $ f a b
liftScmFrac2 f (ScmReal    a) (ScmReal    b) = Right . ScmReal   $ f a b
liftScmFrac2 f (ScmComplex a) (ScmComplex b) = Right . ScmComplex$ f a b

-- Rat+{Real,Complex}
liftScmFrac2 f (ScmRat     a) (ScmReal    b) = Right . ScmReal   $ f (fRat a) b
liftScmFrac2 f (ScmReal    a) (ScmRat     b) = Right . ScmReal   $ f a (fRat b)
liftScmFrac2 f (ScmRat     a) (ScmComplex b) = Right . ScmComplex$ f (fRat a) b
liftScmFrac2 f (ScmComplex a) (ScmRat     b) = Right . ScmComplex$ f a (fRat b)

-- Real+Complex
liftScmFrac2 f (ScmReal    a) (ScmComplex b) = Right . ScmComplex$ f (a :+ 0) b
liftScmFrac2 f (ScmComplex a) (ScmReal    b) = Right . ScmComplex$ f a (b :+ 0)

liftScmFrac2 _ _ _ = fail "liftScmFrac2 :: internal error"

--- Real
liftScmReal2 :: (forall a. Real a => a -> a -> a)
             -> (ScmValue -> ScmValue -> ErrOr ScmValue)
liftScmReal2 f x y =
   case pairScmReal x y of
        Right (ScmReal a, ScmReal b) -> Right . ScmReal $ f a b
        Right (ScmRat  a, ScmRat  b) -> Right . ScmRat  $ f a b
        Right (ScmInt  a, ScmInt  b) -> Right . ScmInt  $ f a b
        Left s                       -> Left s
        Right _                      -> fail "liftScmReal2 :: internal error"

liftScmRealA2 :: (forall a. Real a => a -> a -> b)
              -> (ScmValue -> ScmValue -> ErrOr b)
liftScmRealA2 f x y =
   case pairScmReal x y of
        Right (ScmReal a, ScmReal b) -> Right $ f a b
        Right (ScmRat  a, ScmRat  b) -> Right $ f a b
        Right (ScmInt  a, ScmInt  b) -> Right $ f a b
        Left s                       -> Left s
        Right _                      -> fail "liftScmRealA2 :: internal error"

--- RealFrac
liftScmRealFracB2 :: (forall a. RealFrac a => a -> a -> b)
                  -> (ScmValue -> ScmValue -> ErrOr (Bool, b))
liftScmRealFracB2 f x y =
   case pairScmReal x y of
        Right (ScmReal a, ScmReal b) -> Right (False, f a b)
        Right (ScmRat  a, ScmRat  b) -> Right (True,  f a b)
        Right (ScmInt  a, ScmInt  b) ->
           Right (True,  f (fInt a) (fInt b :: Rational))
        Left s                       -> Left s
        Right _                      ->
           fail "liftScmRealFracB2 :: internal error"

pairScmReal :: ScmValue -> ScmValue -> ErrOr (ScmValue, ScmValue)
pairScmReal (ScmComplex a) _ | imagPart a /= 0 = fail "pairScmReal :: complex"
pairScmReal _ (ScmComplex b) | imagPart b /= 0 = fail "pairScmReal :: complex"

pairScmReal a@(ScmInt     _) b@(ScmInt     _) = Right (a,b)
pairScmReal a@(ScmRat     _) b@(ScmRat     _) = Right (a,b)
pairScmReal a@(ScmReal    _) b@(ScmReal    _) = Right (a,b)
pairScmReal   (ScmComplex a)   (ScmComplex b) =
   Right $ ((,) `on` ScmReal . realPart) a b

-- Int+{Rat,Real,Complex}
pairScmReal (ScmInt     a) (ScmRat     b) = Right $ ((,)`on`ScmRat)  (fInt a) b
pairScmReal (ScmRat     a) (ScmInt     b) = Right $ ((,)`on`ScmRat)  a (fInt b)
pairScmReal (ScmInt     a) (ScmReal    b) = Right $ ((,)`on`ScmReal) (fInt a) b
pairScmReal (ScmReal    a) (ScmInt     b) = Right $ ((,)`on`ScmReal) a (fInt b)
pairScmReal (ScmInt     a) (ScmComplex b) =
   Right $ ((,)`on`ScmReal) (fInt a) (realPart b)
pairScmReal (ScmComplex a) (ScmInt     b) =
   Right $ ((,)`on`ScmReal) (realPart a) (fInt b)

-- Rat+{Real,Complex}
pairScmReal (ScmRat     a) (ScmReal    b) = Right $ ((,)`on`ScmReal) (fRat a) b
pairScmReal (ScmReal    a) (ScmRat     b) = Right $ ((,)`on`ScmReal) a (fRat b)
pairScmReal (ScmRat     a) (ScmComplex b) =
   Right $ ((,)`on`ScmReal) (fRat a) (realPart b)
pairScmReal (ScmComplex a) (ScmRat     b) =
   Right $ ((,)`on`ScmReal) (realPart a) (fRat b)

-- Real+Complex
pairScmReal (ScmReal    a) (ScmComplex b) =
   Right $ ((,)`on`ScmReal) a (realPart b)
pairScmReal (ScmComplex a) (ScmReal    b) =
   Right $ ((,)`on`ScmReal) (realPart a) b

pairScmReal _ _ = fail "pairScmReal :: internal error"

pairScmComplex :: ScmValue -> ScmValue -> ErrOr (ScmValue, ScmValue)
pairScmComplex a@(ScmInt     _) b@(ScmInt     _) = Right (a,b)
pairScmComplex a@(ScmRat     _) b@(ScmRat     _) = Right (a,b)
pairScmComplex a@(ScmReal    _) b@(ScmReal    _) = Right (a,b)
pairScmComplex a@(ScmComplex _) b@(ScmComplex _) = Right (a,b)

-- Int+{Rat,Real,Complex}
pairScmComplex (ScmInt     a) (ScmRat     b) = Right$((,)`on`ScmRat) (fInt a) b
pairScmComplex (ScmRat     a) (ScmInt     b) = Right$((,)`on`ScmRat) a (fInt b)
pairScmComplex (ScmInt     a) (ScmReal    b) =
   Right $ ((,)`on`ScmReal) (fInt a) b
pairScmComplex (ScmReal    a) (ScmInt     b) =
   Right $ ((,)`on`ScmReal) a (fInt b)
pairScmComplex (ScmInt     a) (ScmComplex b) =
   Right $ ((,)`on`ScmComplex) (fInt a) b
pairScmComplex (ScmComplex a) (ScmInt     b) =
   Right $ ((,)`on`ScmComplex) a (fInt b)

-- Rat+{Real,Complex}
pairScmComplex (ScmRat     a) (ScmReal    b) =
   Right $ ((,)`on`ScmReal) (fRat a) b
pairScmComplex (ScmReal    a) (ScmRat     b) =
   Right $ ((,)`on`ScmReal) a (fRat b)
pairScmComplex (ScmRat     a) (ScmComplex b) =
   Right $ ((,)`on`ScmComplex) (fRat a) b
pairScmComplex (ScmComplex a) (ScmRat     b) =
   Right $ ((,)`on`ScmComplex) a (fRat b)

-- Real+Complex
pairScmComplex (ScmReal    a) (ScmComplex b) =
   Right $ ((,)`on`ScmComplex) (a:+0) b
pairScmComplex (ScmComplex a) (ScmReal    b) =
   Right $ ((,)`on`ScmComplex) a (b:+0)

pairScmComplex _ _ = fail "pairScmComplex :: internal error"

fInt :: Num a => Integer -> a
fInt = fromInteger

fRat :: Fractional a => Rational -> a
fRat = fromRational

-- }}}

max_stars_repo_name: Deewiant/haschoo

max_issues_repo_name: Deewiant/haschoo

max_forks_repo_name: Deewiant/haschoo


text:
-- ch3 example 3.17

import Statistics.Distribution
import Statistics.Distribution.Hypergeometric

h :: HypergeometricDistribution
h = hypergeometric m l k
  where
    m =  4 -- r
    l = 20 -- N of population size
    k =  5 -- n of samplie size

p :: Int -> Double
p = probability h

main :: IO ()
main = do
  putStrLn $ "p2 + p3 + p3 = " ++ show (p 2 + p 3 + p 4)
  putStrLn $ "1 - p0 - p1 =  " ++ show (1 - p 0 - p 1)
  putStrLn $ "1 - P(X<=1) =  " ++ show (complCumulative h 1)

  putStrLn $ "Indeed p0 and p1 are: " ++ show (p 0) ++ " and " ++ show (p 1)
  putStrLn $ "The mean = " ++ show (mean h) ++ " variance = " ++ show (variance h)
          ++ " and stdDev = " ++ show (stdDev h)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-|
Module: Class.Mappable
Description: Map functions over data as polymorphically as possible.
License: 0BSD
Maintainer: wygulmage@users.noreply.github.com
Stability: Experimental
Portability: GHC

'Mappable' lets you 'map' over various structures. It's intended to unify all the functions named @map@ ('Data.List.map', 'Data.Set.map', 'Data.Text.map', etc.). To do this it uses type families.
'Data.Functor.Functor' and 'Data.MonoTraversable.MonoFunctor' provide simpler interfaces. But, 'Functor' 'fmap' requires a fully polymorphic container, while 'Data.MonoTraversable.MonoFunctor' 'Data.MonoTraversable.omap' can only map with endofunctions (and neither provides 'mapCoerce' yet).
-}

{-# LANGUAGE NoImplicitPrelude
           , ScopedTypeVariables
           , TypeFamilies
           , UndecidableInstances
           , FlexibleContexts
           , FlexibleInstances
           , DefaultSignatures
           , RankNTypes
           , MultiParamTypeClasses
           , GeneralizedNewtypeDeriving
   #-}

module Class.Mappable (
Mappable (..), omap,
ValidMorph,
WrapMappable (..),
WrapFunctor (..),
) where

import Control.Arrow (Arrow)
import Control.Applicative
import Control.Monad (Monad)
import Control.Monad.ST (ST)
import qualified Control.Exception as Exception
import Data.Ord
import Data.Semigroup
import qualified Data.Semigroup as Semigroup
import qualified Data.Monoid as Monoid
import qualified Data.Functor as Data
import Data.Bool (Bool)
import Data.Char
import Data.Int (Int)
import Data.Word (Word8)
import qualified Data.ByteString as Bytes
import qualified Data.ByteString.Lazy as Lazy.Bytes
import qualified Data.Text as Text
import qualified Data.Text.Lazy as Lazy.Text
import qualified Data.Sequence as Seq
import qualified Data.Tree as Rose
import qualified Data.Set as Set
import qualified Data.IntSet as IntSet
import qualified Data.Map as Map
import qualified Data.IntMap as IntMap
import Data.Complex
import Data.Tagged
import Data.Maybe (Maybe)
import Data.Either
import Data.Proxy
import Data.Functor.Identity
import qualified Data.List.NonEmpty as NonEmpty
import Data.Coerce (Coercible, coerce)
import Data.Kind

import GHC.Arr (Array)
import GHC.Conc (STM)

import Prelude (IO)

import Text.ParserCombinators.ReadP (ReadP)
import Text.ParserCombinators.ReadPrec (ReadPrec)

newtype WrapFunctor m a = WrapFunctor (m a)
   deriving (Data.Functor, Applicative, Alternative, Monad)

newtype WrapMappable t = WrapMappable t

type family ObjectDefault (ma :: Type) where ObjectDefault (m (a :: Type)) = a
{- ^ the rightmost parameter of a type -}

type family MorphDefault (x :: Type) (y :: Type) where
   MorphDefault (polymorphic a) b = polymorphic b
   MorphDefault monomorphic b     = monomorphic
{- ^ Try to replace the rightmost parameter of a type. If the type is polymorphic, the rightmost parameter is replaced; otherwise the original type is returned. -}

class    (Object (Morph mb a) ~ a)=> ValidMorph (mb :: Type) (a :: Type)
instance (Object (Morph mb a) ~ a)=> ValidMorph (mb :: Type) (a :: Type)

class (Morph mb (Object mb) ~ mb)=> Mappable (mb :: Type) where
{-^ @Mappable@ extends the power of 'Functor' to monomorphic containers and containers that constrain their contents.
To create an instance of @Mappable@ for a 'Functor', write e.g. @instance Mappable [b]@.
To create an instance of @Mappable@ for a monomorphic type, write e.g.
@instance Mappable 'Text.Text' where
   type Object 'Text.Text' = 'Char'
   'map' = 'Text.map'
@
To create an instance for a type that constraints its objects, you need to include a definition of 'mapCoerce'.
@instance (Ord b)=> Mappable ('Set.Set' b) where
   'map' = 'Set.map'
   'mapCoerce' = 'map' 'coerce'
@
To define an instance where the 'Object' is phantom, use a kind annotation: @instance Mappable ('Proxy' (o :: 'Type'))@.

Laws:
   1. for all f g. @map (g '$!') . map f@ === @map ((g '$!') . f)@, as long as equal ('==') objects produced by @f@ are indistinguishable to 'g'.
   2. @map 'id'@ === 'id'
   3. 'mapCoerce' === @'map' 'coerce'@
The first law forces @g@'s argument to accommodate both strict and lazy containers, and has the condition on equal objects to accommodate both 'Arg' and 'Set.Set'.
-}
   type Morph mb (a :: Type) :: Type
   {-^ Try to change the type of an object in a container. @Morph [a] b = [b]@; @Morph 'Text.Text' b = 'Text.Text'. -}
   type Morph mb a = MorphDefault mb a
   type Object mb :: Type
   {-^ The type of the objects in the container. By default it's the type's rightmost parameter. -}
   type Object mb = ObjectDefault mb

   map ::
      forall a. (ValidMorph mb a)=>
      (a -> Object mb) -> Morph mb a -> mb
   {-^ Map a function over all the 'Object's in a container. -}

   mapCoerce ::
      forall a. (ValidMorph mb a, Coercible a (Object mb))=>
      Morph mb a -> mb
   {-^ Coerce all the 'Object's in a container to a new type. -}
   -- mapCoerce = map coerce

   default map ::
      (ValidMorph mb a, Morph mb a ~ m a, m (Object mb) ~ mb, Data.Functor m)=>
      (a -> Object mb) -> Morph mb a -> mb
   map = Data.fmap
   {-# INLINE map #-}

   default mapCoerce ::
      (ValidMorph mb a, Coercible (Morph mb a) mb)=> Morph mb a -> mb
   mapCoerce = coerce
   {-# INLINE mapCoerce #-}

omap :: (Mappable t)=> (Object t -> Object t) -> t -> t
omap = map
{-# INLINE omap #-}

------ Instances ------

-- TODO: Add instances for lifted `Functor`s (e.g. `Ap m b`). Should these rely on the `Functor` instance or the `Mappable` instance? In favor of `Mappble` is being able to use constrained containers like `Set`.

--- State Functors ---
instance Mappable (IO b)
instance Mappable (ST s b)
instance Mappable (STM b)

--- Empty Containers ---
instance Mappable (Proxy (b :: Type))

--- Identity Containers ---
instance Mappable (Identity b)
instance Mappable (Down b)
instance Mappable (Dual b)
instance Mappable (Semigroup.First b)
instance Mappable (Semigroup.Last b)
instance Mappable (Max b)
instance Mappable (Min b)
instance Mappable (Product b)
instance Mappable (Sum b)

instance Mappable Semigroup.All where
   type Object Semigroup.All = Bool
   map = coerce

instance Mappable Semigroup.Any where
   type Object Semigroup.Any = Bool
   map = coerce

--- Affine Containers ---
instance Mappable (Maybe b)
instance Mappable (Monoid.First b)
instance Mappable (Monoid.Last b)

--- Fixed-size Vectors ---
instance Mappable (Complex b)

--- Sequences ---
instance Mappable [b]
instance Mappable (ZipList b)

instance Mappable (NonEmpty.NonEmpty b)

instance Mappable (Seq.Seq b)
instance Mappable (Seq.ViewL b)
instance Mappable (Seq.ViewR b)

instance Mappable (Rose.Tree b)

instance Mappable (Array i b)

instance Mappable Bytes.ByteString where
   type Object Bytes.ByteString = Word8
   map = Bytes.map

instance Mappable Lazy.Bytes.ByteString where
   type Object Lazy.Bytes.ByteString = Word8
   map = Lazy.Bytes.map

instance Mappable Text.Text where
   type Object Text.Text = Char
   map = Text.map

instance Mappable Lazy.Text.Text where
   type Object Lazy.Text.Text = Char
   map = Lazy.Text.map

--- Sets ---
instance (Ord b)=> Mappable (Set.Set b) where
   map = Set.map
   mapCoerce = map coerce
{-^ WARNING: @map (g $!) . map f@ is not equivalent to @map ((g $!) . f)@ when @g@ can distinguish between different but equal ('==') values produced by @f@. For example, @map (\ ('Arg' () n) -> n) . map ('Arg' ())@ is likely to give a very different result than @map ((\ ('Arg' () n) -> n) . 'Arg' ())@. -}

instance Mappable IntSet.IntSet where
   type Object IntSet.IntSet = Int
   map = IntSet.map

--- Bifunctors ---
instance Mappable (Const c (b :: Type))

instance Mappable (Tagged c b)

instance Mappable (Either c b)

instance Mappable (c, b)
instance Mappable (Arg c b)

--- Functions ---
instance Mappable (a -> b)

instance Mappable (Exception.Handler b)

instance Mappable (ReadP b)
instance Mappable (ReadPrec b)

--- Maps ---
instance Mappable (Map.Map i b)
instance Mappable (IntMap.IntMap b)

--- Wrapped Instances ---
instance (Arrow p)=> Mappable (WrappedArrow p c b) where
   mapCoerce = map coerce

instance (Data.Functor m)=> Mappable (WrapFunctor m a) where
   mapCoerce = map coerce

instance (Mappable t)=> Mappable (WrapMappable t) where
   type Object (WrapMappable t) = Object t
   type Morph (WrapMappable t) a = WrapMappable (Morph t a)
   map f (WrapMappable xs) = WrapMappable (map f xs)
   mapCoerce (WrapMappable xs) = WrapMappable (mapCoerce xs)

max_stars_repo_name: wygulmage/class-mappable

max_issues_repo_name: wygulmage/class-mappable

max_forks_repo_name: wygulmage/class-mappable