CRC.hs

{- Computes the probability of observing a Schrodinger CRC using
QuickCheck-generated random data.  Requires a patch to QuickCheck available at
http://www.cs.indiana.edu/~lepike/support/QcBiphase/pike-quotient-qc .

Build with
> ghc --make CRC.hs

Run with
> ./CRC

Change the value of main to your desired computation.

Lee Pike <lee-pike-@-gmail-.-com> (remove '-'s)
Copyright, 2009
BSD3 License
-}

module Main

where

import Debug.Trace

import System.Random.Mersenne

import Data.List
import Test.QuickCheck

import Test.QuickCheck
import Test.QuickCheck.Monadic
import Test.QuickCheck.Gen

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

-- Helper for Mersenne randoms for ints, generates a random in [low, high]
randomHlp :: (MTRandom a, Integral a) => (a, a) -> IO a
randomHlp (low, high) = do r <- randomIO::IO Double -- generates in [0, 1) 
                           return $ low + (round (r * (fromIntegral $ high - low)))

type Bits = Integer 

-- |Bitvector type (however, only lists of 0s and 1s are valid).
type BV = [Bits]

-- |Changes a bitstring (list of 0's and 1's), with the most sig bit at the
-- |head, to its corresponding natural number representation.
fromBV :: BV -> Integer
fromBV y = foldl (+) 0 $ zipWith (*) (reverse y) $ map (\x -> 2^x) $ [0..]

-- |Converts a natural into a bitstring (list of 0's and 1's), with the most sig
-- |bit at the head.
toBV :: Integer -> BV
toBV x = reverse $ toBVHelp 1 x
    where toBVHelp i rst =
              let n = 2^i
                  r = rem rst n
                  d = div rst n
              in if d == 0
                   then if r == 0 then [0] else [1]
                   else if d == 1 && r == 0
                          then [0,1]
                          else if r == 0
                                 then 0:(toBVHelp (i+1) rst)
                                 else 1:(toBVHelp (i+1) (rst-r))

-- Sanity check.
prop_bitconvert x = x >= 0 ==> (fromBV . toBV) x == x

-- |Specific CRC polynomials (see
-- |<http://en.wikipedia.org/wiki/Cyclic_redundancy_check>).
ccitt8, crc4itu,crc8atm, parity :: BV
ccitt8 = [1,1,0,0,0,1,1,0,1]
ccitt16 = [1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1] --0x8005 x^16 + x^15 + x^2 + 1
crc4itu = [1,0,0,1,1]
crc8atm = [1,0,0,0,0,0,1,1,1]
crc8 = [1,1,1,0,1,0,1,0,1]
crc80x17 = [1,0,0,0,0,0,1,1,1]
parity = [1,1]
usb5 = [1,0,0,1,0,1] -- Standard and optimal for words of 11 bits.  Koopman &
                     -- Chakravarty "CRC Polynomial Selection for Embedded
                     -- Networks," DSN 2004.
can = [1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1] -- For the CAN bus.  HD=6 up to 112-bit data words

-- |Removes the leading 0s from a bitvector, but don't return an empty bitvector.
shrinkBV :: BV -> BV
shrinkBV bv = let bv' = dropWhile (== 0) bv
              in if null bv' then [0] else bv'

-- |Takes a bitvector and the CRC polynomial and returns the bitvector divided
-- by the polynomial over GF(2).  For a polynomial of degree n, add n zeros.
crc :: BV -- ^ bitvector
    -> BV -- ^ polynomial
    -> BV
crc bv p = crc' (bv ++ (replicate (length p - 1) 0)) p

crc' bv p = let bv' = shrinkBV bv
                xor = let z = p ++ [0,0..]
                      in  zipWith (\a b -> (a+b) `mod` 2) bv' z
            in if length bv' < length p then reverse $ take (length p - 1) (reverse bv)
                 else crc' xor p 

-- |CRC check where the bitvector and polynomial are represented as naturals.
crcFromInt :: Integer -- ^ bitvector
       -> Integer -- ^ polynomial
       -> BV
crcFromInt nbv np = crc (toBV nbv) (toBV np)

-- Sanity check for parity.
prop_parity n = crcFromInt n (fromBV parity) == crc (toBV n ++ [1,1]) parity

-- | Takes a bitvector bv, a Hamming Distance hd, and returns a random bitvector
-- | such that some random i number of 0s in bv are flipped to 1, where i >= hd
-- | (and less than the total number of 0s!).  This function should only be
-- | called by bvs with at least hd 0s.
rndBv0to1 :: BV -> Int -> IO BV
rndBv0to1 bv hd =
    let idxs0s = elemIndices 0 bv
    in do i <- randomHlp (hd, length idxs0s) -- total number of bits to flip
          idxs <- whichBitsFlipped i idxs0s -- which bits are we flipping?
          return $ replaceBits bv idxs
              where whichBitsFlipped j xs =
                      if j == 0 then return []
                        else do r <- randomHlp (0, (length xs) - 1)
                                rst <- whichBitsFlipped (j-1) (delete (xs !! r) xs)
                                return $ (xs !! r):rst

-- | Takes a bitvector bv and a list of indices ls and returns a bitvector such
-- | that for each index in ls, the bit at that index is flipped.
replaceBits :: BV -> [Int] -> BV
replaceBits bv ls = let ls' = sort ls
                    in replaceBits' 0 ls'
    where replaceBits' k indicies =
              case indicies of []   -> drop k bv
                               i:is ->    (drop k $ take i bv)
                                       ++ (bitFlip $ bv !! i):(replaceBits' (i+1) is)

-- |Flip a bit.
bitFlip :: Bits -> Bits
bitFlip x = if x == 1 then 0 else 1

-- |Takes a bitvector bv and polynomial, and computes the CRC.  It then randomly
-- |flips from 0s to 1s [hamming distance (HD) <= i <= # of zeros] in word +
-- |crc.  It returns false if the corrupted CRC is valid for the corrupted
-- |message.
schrodingerCRC :: BV -- ^ bitvector
               -> BV -- ^ crc polynomial
               -> Int -- ^ hamming distance
               -> Bool -- ^ show output
               -> IO Bool
schrodingerCRC bv poly hd out =
    let fcs = crc bv poly -- get the FCS on the original
        wrdfcs = bv ++ fcs
    in if (length $ filter (== 0) wrdfcs) < hd
       then do if out then putStrLn $ "0s less than Hamming Distance: " ++ (show wrdfcs)
                 else return ()
               return True
       else do wrdfcs' <- rndBv0to1 wrdfcs 0 -- flip some bits from 0 to 1
               let rcvfcs = drop (length bv) wrdfcs'
               let bv' = take (length bv) wrdfcs'
               let fcs' = crc bv' poly
               if out
                 then do putStrLn $ "Original bitvector: " ++ (show bv)
                                      ++ "  Original FCS: " ++ (show fcs)
                         putStrLn $ "Received bitvector: "
                                      ++ (show bv') ++ "  Received FCS: " ++ (show rcvfcs)
                         putStrLn $ "Actual FCS: " ++ (show fcs')
                 else return ()
               return $ rcvfcs /= fcs' || bv' == bv

prop_noschrodinger poly i hd out =
  monadicIO $ do rnd <- pick $ genBVs i
                 b <- run $ schrodingerCRC (toBV rnd) poly hd out
                 assert b

-- Generates a random BV up to some length i.
genBVs :: Integer -> Gen Integer
genBVs i = choose (0::Integer, (2^i-1))

-- Execute
main :: IO ()
-- ~3.29e-2 == 0.0329 
main = quickCheckQuotientWith (stdArgs {maxSuccess = 1000}) (prop_noschrodinger usb5 11 3 True)

-- result: ~39/1000000 == 3.9e-5
--main = quickCheckQuotientWith (stdArgs {maxSuccess = 1000000}) (prop_noschrodinger can 64 6 False)

-------------------
-- Functions for calculating the probability of bit corruptions.

-- Factorial.
fac :: Integer -> Integer
fac n = product [1..n]

top :: Integer -> Integer
top n = round $ (fromIntegral n)/2

prob :: Double -> Integer -> Integer -> Double
prob p flipped n = p^flipped * (1-p)^(n - flipped)

choose' :: Integer -> Integer -> Double
choose' flipped n =
    fromRational $ (fromIntegral (fac n)) / (fromIntegral ((fac flipped) * fac (n - flipped)))

summation :: Double -> Integer -> Integer -> Double
summation p hd n =
    sum $ map f [hd..t]
        where t = top n
              f = \x -> (prob p x t) * (choose' x t)

-- binomial formula.  Takes a n trials, prob per trial, and number of outcomes you want.
binomial n p x = let fi i   = fromInteger i
                     choose = (fi (fac n)) / fi ((fac x) * (fac (n - x)))
                     px     = p ** (fi x)
                     notpx  = (1 - p) ** (fi (n - x))
                 in  choose * px * notpx

-- prob X <= x
noMore n p x = foldl1' (+) $ map (binomial n p) [0..x]

-- Prob X >= x
atLeast n p x = foldl1' (+) $ map (binomial n p) [x..n]

---