Implementation of GF(2^16)

F16.hs

{-# LANGUAGE BangPatterns #-}

module F16 where

import Data.Int
import Data.Bits
import Data.List
import Data.Ratio
import Numeric
import Data.Ord
import Debug.Trace

newtype F16 = F16 Integer

bitsPerDigit = 4
ones = iterate (\x -> (shiftL x bitsPerDigit) .|. 1) 0

masklo = ones !! 16
maskhi = shiftL masklo (16*bitsPerDigit)

instance Eq F16 where
  F16 a == F16 b = a == b

reduce :: Integer -> Integer
reduce z = if zz == 0 then (z .&. masklo) else reduce x
  where zz = z .&. maskhi
        r1 = shiftR zz ((16-5)*bitsPerDigit)
        r2 = shiftR zz ((16-3)*bitsPerDigit)
        r3 = shiftR zz ((16-1)*bitsPerDigit)
        r4 = shiftR zz (16*bitsPerDigit)
        r = xor r1 (xor r2 (xor r3 r4))
        x = (xor (xor z zz) r) .&. (maskhi .|. masklo)

instance Num F16 where
  (+) (F16 a) (F16 b) = F16 ((a+b) .&. masklo)
  (*) (F16 a) (F16 b) = F16 (reduce (a*b))
  (-)                 = (+)
  negate              = id
  abs = error "abs not implemented for F16"
  signum = error "signum not implemented for F16"
  fromInteger a = F16 (foldl' go 0 [0..15])
    where go s k = if testBit a k then s .|. (bit (k*bitsPerDigit)) else s

instance Fractional F16 where
  recip (F16 0) = error "division by zero"
  recip x = x ^ 65534
  fromRational r = fromInteger p * (recip (fromInteger q))
    where p = numerator r
          q = denominator r

instance Enum F16 where
  toEnum = fromIntegral
  fromEnum (F16 a) = go 1 0 a
    where go _  s 0 = s
          go !b !s a = if testBit a 1
                         then go (2*b) (s+b) (shiftR a 1)
                         else go (2*b) s     (shiftR a 1)

toCoeffs :: Integer -> [Int]
toCoeffs a = map fromEnum terms
  where
    terms = map go bits 
    bits = takeWhile (/=0) $ iterate (\x -> shiftR x bitsPerDigit) a
    go a = testBit a 0

showF16 :: Integer -> String
showF16 a = intercalate " + " terms
  where go (a,k)  = if testBit a 0 then xterm k else ""
        xterm 0 = "1"
        xterm k = "x^" ++ show (k :: Int)
        bits = takeWhile (/=0) $ iterate (\x -> shiftR x bitsPerDigit) a
        terms = filter (not.null) $ map go $ zip bits [0..] :: [String]

instance Show F16 where
  show (F16 a) = showF16 a -- everyOther (showOct a "")
    where
      everyOther (a:_:rest) = [a] ++ everyOther rest
      everyOther as = as

-- find generator
findgen = [ g | g <- map fromInteger [1..32767], let o = order g, o == 65535 ]

order :: F16 -> Int
order g = 1 + (length $ takeWhile test $ iterate (*g) g)
  where test 0 = False -- should never happen
        test 1 = False
        test _ = True

TestF16.hs

module TestF16 where

import Math.Polynomial
import Math.Core.Field hiding (F16, powers)
import Control.Monad
import Numeric

import F16

r16 = poly LE [(1::F2), 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
r4  = poly LE [(1::F2), 1, 0, 0, 1]

test1 =
  let p = poly LE [ 1, 1 ]
      a = remPoly (powPoly p 15) r4
  in a

test2 coeffs =
  let p = poly LE coeffs
      op x = remPoly (multPoly p x) r4
  in mapM_ print $ take 16 $ iterate op one 

polyPowers coeffs =
  let p = poly LE coeffs
      op x = remPoly (multPoly p x) r16
  in iterate op one

compF16 :: F16 -> Poly F2 -> Bool
compF16 (F16 x) p =
  let cs16 = map fromIntegral (toCoeffs x)
      csp = polyCoeffs LE p
  in cs16 == csp

-- polynomial power mod
powerMod m p 0 = one
powerMod m p k =
  let (q,r) = quotRem k 2
      s = powerMod m p q
      s' = remPoly (multPoly s s) m
      s'' = if r == 0 then s' else remPoly (multPoly s' p) m
  in s''

comparePowers :: Int -> F16 -> IO ()
comparePowers n a@(F16 x) = do
  let coeffs = map fromIntegral (toCoeffs x)
      p1 = [ powerMod r16 (poly LE coeffs) k | k <- [0..] ] -- take n (polyPowers coeffs)
      p2 = take n (iterate (*a) 1)
  forM_ (zip3 [(0::Int)..] p1 p2) $ \(k,a,p) -> do
    putStrLn $ "n = " ++ show k ++ ": " ++ show (compF16 p a)

-- test n powers of an F16 value
test3 :: Int -> F16 -> IO Bool
test3 n a@(F16 x) = do
  let 
      coeffs = map fromIntegral (toCoeffs x)
      polys = [ powerMod r16 (poly LE coeffs) k | k <- [0..] ]
      f16s = take n (iterate (*a) 1)
      tests = zip3 [0..] polys f16s

      loop [] = return True
      loop (t:ts) = do b <- runTest t
                       if b then loop ts else return False
      runTest (k,a,p) = do
        let ok = compF16 p a
        if ok then return True
              else do putStrLn $ "Failed on k = " ++ show k
                      return False
  loop tests