## 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)
-}

module Main

where

import Debug.Trace

import System.Random.Mersenne

import Data.List
import Test.QuickCheck

import Test.QuickCheck
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
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)
++ (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]
```