galois-field-0.2.0: src/PolynomialRing.hs

module PolynomialRing
  ( Polynomial(..)
  , cut
  , polyInv
  , polyMul
  , polyQR
  , toPoly
  ) where

import Protolude

import Test.Tasty.QuickCheck (Arbitrary(..))
import Text.PrettyPrint.Leijen.Text (Pretty(..))

import GaloisField (GaloisField(..))

newtype Polynomial k = X [k]
  deriving (Eq, Generic, NFData, Show)

instance (Arbitrary k, GaloisField k) => Arbitrary (Polynomial k) where
  arbitrary = toPoly <$> arbitrary

instance GaloisField k => Num (Polynomial k) where
  X xs + X ys   = X (polyAdd xs ys)
  {-# INLINE (+) #-}
  X xs * X ys   = X (polyMul xs ys)
  {-# INLINE (*) #-}
  X xs - X ys   = X (polySub xs ys)
  {-# INLINE (-) #-}
  negate (X xs) = X (map negate xs)
  {-# INLINE negate #-}
  fromInteger n = X (let m = fromInteger n in if m == 0 then [] else [m])
  {-# INLINABLE fromInteger #-}
  abs           = panic "not implemented."
  signum        = panic "not implemented."

instance GaloisField k => Pretty (Polynomial k) where
  pretty (X xs) = pretty xs

polyAdd :: GaloisField k => [k] -> [k] -> [k]
polyAdd xs     []     = xs
polyAdd []     ys     = ys
polyAdd (x:xs) (y:ys) = let z  = x + y
                            zs = polyAdd xs ys
                        in if z == 0 && null zs then [] else z : zs
{-# INLINE polyAdd #-}

polyMul :: GaloisField k => [k] -> [k] -> [k]
polyMul _      [] = []
polyMul []     _  = []
polyMul (x:xs) ys = let ws = map (* x) ys
                        zs = polyMul xs ys
                    in if null xs then ws else polyAdd ws (0 : zs)
{-# INLINE polyMul #-}

polySub :: GaloisField k => [k] -> [k] -> [k]
polySub xs     []     = xs
polySub []     ys     = map negate ys
polySub (x:xs) (y:ys) = let z  = x - y
                            zs = polySub xs ys
                        in if z == 0 && null zs then [] else z : zs
{-# INLINE polySub #-}

polyQR :: forall k . GaloisField k => [k] -> [k] -> ([k], [k])
polyQR xs ys = polyGCD ([], xs)
  where
    m :: Int
    m = length ys
    polyGCD :: ([k], [k]) -> ([k], [k])
    polyGCD qr@(qs, rs)
      | n < 0     = qr
      | otherwise = polyGCD (polyAdd qs ts, polySub rs (polyMul ts ys))
      where
        n :: Int
        n = length rs - m
        ts :: [k]
        ts = replicate n 0 ++ [last rs / last ys]
        last :: [k] -> k
        last []     = 0
        last [z]    = z
        last (_:zs) = last zs
{-# INLINE polyQR #-}

polyInv :: forall k . GaloisField k => [k] -> [k] -> Maybe [k]
polyInv [x] _ = Just [recip x]
polyInv xs ps = case extGCD (ps, xs) of
  ([y], (ys, _)) -> Just (map (/ y) ys)
  _              -> Nothing
  where
    extGCD :: ([k], [k]) -> ([k], ([k], [k]))
    extGCD (y, []) = (y, ([], [1]))
    extGCD (y, x)  = (g, (polySub t (polyMul s q), s))
      where
        (q, r)      = polyQR y x
        (g, (s, t)) = extGCD (x, r)
{-# INLINE polyInv #-}

cut :: GaloisField k => [k] -> [k]
cut = reverse . dropWhile (== 0) . reverse
{-# INLINABLE cut #-}

toPoly :: GaloisField k => [k] -> Polynomial k
toPoly = X . cut
{-# INLINABLE toPoly #-}