Small experimental Haskell library for computations on finite fields

Galois.hs

-- this is file Galois.hs
-- Load it into GHCi

module Galois
  (addGF
  , subGF
  , scaleGF
  , shiftGF
  , mulGF
  , expGF
  , invGF
  , isPrimitiveGF
  , randomGF
  , randomPrimitiveGF
  , smallPrimitiveGF
  , logGF
  , x4
  , y4
  , x7
  , y7
  , x8
  , y8
  , p4
  , p7
  , p8
  , gf4
  , gf7
  , gf8
  ) where

import System.Random
import Data.List
import PolyUtils
import Math.NumberTheory.Primes.Factorisation (factorise)
import Math.NumberTheory.Moduli.Chinese (chineseRemainder)

data GF = GF{basePrime::Integer
            ,power::Integer
            ,irreducible::[Integer]
       } deriving(Eq, Show)

data GFElement
  = GFElement{field::GF
             ,poly::[Integer]
             } deriving (Eq, Show)

addGF :: GFElement -> GFElement -> GFElement
addGF x y
  | field x == field y = GFElement fx z
  | otherwise = error "Elements are from different fields"
  where
    fx = field x
    p = basePrime fx
    z = zipWith (\i j -> mod (i+j) p) (poly x) (poly y)

(&+) :: GFElement -> GFElement -> GFElement
x &+ y = addGF x y 

subGF :: GFElement -> GFElement -> GFElement
subGF x y
  | field x == field y = GFElement fx z
  | otherwise = error "Elements are from different fields"
  where
    fx = field x
    p = basePrime fx
    z = zipWith (\i j -> mod (i-j) p) (poly x) (poly y)

(&-) :: GFElement -> GFElement -> GFElement
x &- y = addGF x y 

scaleGF :: Integer -> GFElement -> GFElement
scaleGF k x = GFElement fx z
  where
    fx = field x
    p = basePrime fx
    z = map (\i -> mod (i*k) p) (poly x)

shiftGF::GFElement -> GFElement
shiftGF x = GFElement fx y
  where
    fx = field x
    px = poly x
    h = head px
    p = basePrime fx
    r = irreducible fx
    y = zipWith (\i j -> mod (i-h*j) p) ((tail px)++[0]) (tail r)

zeroGF::GF -> GFElement
zeroGF f = GFElement f z
  where
    d = power f
    z = take (fromIntegral d) [0,0..]

oneGF::GF -> GFElement
oneGF f = GFElement f z
  where
    d = power f
    z = (take (fromIntegral d-1) [0,0..]) ++ [1]

-- mulGF multiplies u and v by taking all shifted versions of v (that is, X^k*v
-- for all k from 0 to the power of the field), scaling each shift by the
-- corresponding element of u and adding them.

mulGF :: GFElement -> GFElement -> GFElement
mulGF x y
  | field x == field y = u
  | otherwise = error "Elements are from different fields"
  where
    fx = field x
    d = power fx
    fz = zeroGF fx
    allshifts = reverse (take (fromIntegral d) (iterate shiftGF y))
    u = foldl (addGF) fz (zipWith (\z w -> scaleGF z w) (poly x) allshifts)

(&*) :: GFElement -> GFElement -> GFElement
x &* y = mulGF x y 

exp2GF :: GFElement -> Integer -> GFElement
exp2GF x k 
  | k == 0 = oneGF $ field x
  | even k = mulGF z z
  | otherwise = mulGF (mulGF z z) x
  where z = exp2GF x (div k 2)

expGF :: GFElement -> Integer -> GFElement
expGF x k
  | k == 0    = oneGF $ field x
  | otherwise = mulGF t $ expGF (mulGF x x) (div k 2)
  where
    t = if (mod k 2) == 1 then x
        else oneGF $ field x

(&^) :: GFElement -> Integer -> GFElement
x &^ k = expGF x k


invGF :: GFElement -> GFElement
invGF x
  | x == zeroGF fx = error "The zero element is not invertible"
  | otherwise = expGF x k
  where
    fx = field x
    p = basePrime fx
    d = power fx
    k = p^d - 2

-- in a field of order n = p^d, a field element x is primitive if
--    x^((n-1)/r) /= 1
-- for all prime divisors r of n-1

isPrimitiveGF :: GFElement -> Bool
isPrimitiveGF x =
  all (/= (oneGF fx)) (map (\y -> expGF x (div n y)) nd)
    where
      fx = field x
      p = basePrime fx
      d = power fx
      n = p^d - 1
      nd = nub $ primeFactors n

-- discrete logs by Baby step, Giant step:

bgLogGF :: GFElement -> GFElement -> Integer
bgLogGF g h
  | not $ isPrimitiveGF g = error("g is not a primitive")
  | h == zeroGF f       = error("zero is not in the multiplicative group")
  | field g /= field h  = error("Elements are not in the same field")
  | otherwise           = k
  where
    f = field g
    m = (basePrime f)^(power f) - 1
    w = ceiling $ sqrt $ fromIntegral m
    a = map (\j -> expGF g $ w*j) [0,1..w-1]
    b = map (\i -> mulGF h $ expGF (invGF g) i) [0,1..w-1]
    ab = head $ intersect a b
    i  = case (elemIndex ab a) of
           Just k -> toInteger k
           Nothing -> -1
    j  = case (elemIndex ab b) of
           Just k -> toInteger k
           Nothing -> -1
    k = i*w+j

-- subLogGF is a helper function for the Pohlig-Hellman algorithm, it finds logs
-- in a smaller multiplicative subgroup, where the order m is given

subLogGF :: GFElement -> GFElement -> Integer -> Integer
subLogGF g h m
--  | not $ isPrimitiveGF g = error("g is not a primitive element")
  | h == zeroGF (field g) = error("zero is not in the multiplicative group")
  | field g /= field h  = error("Elements are not in the same field")
  | otherwise           = k
  where
    w = ceiling $ sqrt $ fromIntegral m
    a = map (\j -> expGF g $ w*j) [0,1..w-1]
    b = map (\i -> mulGF h $ expGF (invGF g) i) [0,1..w-1]
    ab = head $ intersect a b    -- by definition and construction this exists
    i  = case (elemIndex ab a) of
           Just k -> toInteger k
           Nothing -> -1
    j  = case (elemIndex ab b) of
           Just k -> toInteger k
           Nothing -> -1
    k = i*w+j

-- logGF finds the discrete logarithm by the Pohlig-Hellman algorithm.  This is
-- very efficient if the order of the multiplicative subgroup has small factors.
-- Even so, it can be slow: for a field of size 5^42, computing logs took over
-- 13 minutes.

logGF :: GFElement -> GFElement -> (Maybe Integer)
logGF g h
  | not $ isPrimitiveGF g = error("First element is not a primitive")
  | h == zeroGF f       = error("zero is not in the multiplicative group")
  | field g /= field h  = error("Elements are not in the same field")
  | otherwise           = k
  where
    f = field g
    k = chineseRemainder (zip logs pps)
    logs = zipWith3 (\x y z -> subLogGF x y (z-1)) gs hs pps
    gs = map (\x -> expGF g x) nps
    hs = map (\x -> expGF h x) nps
    n = (basePrime f)^(power f) - 1
    pps = map (\x -> (fst x)^(snd x)) $ factorise n
    nps = map (\x -> div n x) pps
    eliminate = maybe 0 id

-- finds the multiplicative order of an element by iterating over the list of
-- divisors of the size of the cyclic group.  Note: it's very slow!

ord2GF :: GFElement -> Integer
ord2GF x = head $ dropWhile (notOne) divs
  where
    divs = divisors $ (basePrime $ field x)^(power $ field x)-1
    notOne z = (x &^ z) /= (oneGF $ field x)

-- This is a slightly faster "ord" that computes the powers of x on the go,
-- rather than computing the powers to all divisors.  Powers are computed as
-- powers of lower values; for example, x^12 = (x^4)^3 where x^4 has already
-- been computed.  This should reduce the total time.

ordGF :: GFElement -> Integer
ordGF x = divs !! (eliminate $ elemIndex og xs3)
  where
    og = oneGF $ field x
    xs2 = scanl (&^) x (head $ group fs)
    divs = divisors $ (basePrime $ field x)^(power $ field x)-1
    fs = primeFactors $ (basePrime $ field x)^(power $ field x) - 1
    fps = map (scanl (*) 1) (tail $ group fs)
    xs3 = foldl (cartprodGF) xs2 fps
    cartprodGF u v = [i &^ j | i <- u, j <- v]
    eliminate = maybe 0 id

ord3GF :: GFElement -> Integer
ord3GF x = fst $ head $ dropWhile (\z -> snd z /= og) fx
  where
    og = oneGF $ field x
    xs2 = scanl (&^) x (head $ group fs)
    divs = divisors $ (basePrime $ field x)^(power $ field x)-1
    fs = primeFactors $ (basePrime $ field x)^(power $ field x) - 1
    fps = map (scanl (*) 1) (tail $ group fs)
    xs3 = foldl (cartprodGF) xs2 fps
    fx = zip divs xs3
    cartprodGF u v = [i &^ j | i <- u, j <- v]
    
-- creates a random element of the field

randomGF :: GF -> Int -> GFElement
randomGF f seed = GFElement f (randomList d p seed)
  where
    p = basePrime f
    d = power f

-- smallPrimitiveGF finds values by starting at one and working up

smallPrimitiveGF :: GF -> GFElement
smallPrimitiveGF f = head $ dropWhile (\x -> not $ isPrimitiveGF x) gs
  where
    p = basePrime f
    d = power f
    ns = map (\x -> toBase x p) [1,2..]
    gs = map(\x -> GFElement f $ replicate (fromIntegral d - length x) 0 ++ x) ns

-- finds a primitive value by the simple heuristic of producing random values
-- until one is primitive.  Since this is a random operation it needs to be
-- seeded

randomPrimitiveGF :: GF -> Int -> GFElement
randomPrimitiveGF f seed
  | isPrimitiveGF p = p
  | otherwise     = randomPrimitiveGF f (seed+1)
  where
    p = randomGF f seed
  
-- RandomList creates pseudo-random numbers, given the parameters and a seed,

randomList :: Integer -> Integer -> Int -> [Integer]
randomList n k seed = rs::[Integer]
  where
    rs = take (fromIntegral n) $ randomRs (0,k-1) $ mkStdGen seed

-- exPolyGCD computes the GCD of two polynomials a,b modulo p and returns a
-- tuple containing the GCD g, and the Bezout polynomials s,t such that
-- a s + b t = g.
{-
exPolyGCD :: [Integer] -> [Integer] -> [Integer] -> ([Integer],[Integer],[Integer])
exPolyGCD [0] b = ([b], [0], [1])
exPolyGCD a b = let (g, s, t) = exPolyGCD (polyMod b a) a p
                in (g, polyMul (polySub t (polyDiv b a)) s, s)
-}
-- polyAdd and polySub adds and subtracts two polynomials, first ensuring that
-- they are of the same length

polyAdd :: [Integer] -> [Integer] -> Integer -> [Integer]
polyAdd a b p = zipWith (\x y -> mod (x+y) p) a2 b2
  where
    ad = length a
    bd = length b
    a2 | ad >= bd = a
       | otherwise = (replicate (bd - ad) $ toInteger 0) ++ a
    b2 | bd >= ad = b
       | otherwise = (replicate (ad - bd) $ toInteger 0) ++ b

polySub :: [Integer] -> [Integer] -> Integer -> [Integer]
polySub a b p = zipWith (\x y -> mod (x-y) p) a2 b2
  where
    ad = length a
    bd = length b
    a2 | ad >= bd = a
       | otherwise = (replicate (bd - ad) $ toInteger 0) ++ a
    b2 | bd >= ad = b
       | otherwise = (replicate (ad - bd) $ toInteger 0) ++ b

-- polyScale a k p produces ka mod p where k is a scalar and p is a polynomial

polyScale :: [Integer] -> Integer -> Integer -> [Integer]
polyScale a k p = map (\x -> mod (x*k) p) a

{-
polyMul :: [Integer] -> [Integer] -> Integer -> [Integer]
polyMul (a:as) b p
  | all (== 0) a = [0]
  | otherwise    = 
  where 
rMul :: [Integer] -> [Integer] -> [Integer]
rMul (x:xs) y
-}

p4,p7,p8::[Integer]
x4,y4,x7,y7,x8,y8::GFElement
gf4,gf7,gf8,gf12::GF

gf4 = GF 2 4 p4
gf8 = GF 2 8 p8
gf7 = GF 7 4 p7
gf12 = GF 7 12 p12

p4 = [1,0,0,1,1]
x4 = GFElement gf4 [1,0,1,0]
y4 = GFElement gf4 [0,1,1,1]

p8 = [1,0,0,0,1,1,0,1,1]
x8 = GFElement gf8 [0,1,0,1,0,0,1,1]
y8 = GFElement gf8 [1,1,0,0,1,0,1,0]

p7 = [1,3,5,6,2]
x7 = GFElement gf7 [2,3,5,6]
y7 = GFElement gf7 [5,3,0,2]

p12 = [1,0,0,0,2,5,3,2,4,0,5,0,3]
x12 = GFElement gf12 [0,0,0,0,0,0,0,0,0,1,4,6]
y12 = GFElement gf12 [3,5,2,2,1,1,5,0,0,4,6,1]

{-
Some fields for which the order of the multiplicative cyclic group is smooth -
all prime factors are small:

2^36 - 1 = [3,3,3,5,7,13,19,37,73,109]
3^16 - 1 = [2,2,2,2,2,2,5,17,41,193]
3^20 - 1 = [2,2,2,2,5,5,11,11,61,1181]
3^48 - 1 = [2,2,2,2,2,2,5,7,13,17,41,73,97,193,577,769,6481]
3^54 - 1 = [2,2,2,7,13,19,37,109,433,757,8209,19441,19927]
5^18 - 1 = [2,2,2,3,3,3,7,19,31,829,5167]
5^30 - 1 = [2,2,2,3,3,7,11,31,61,71,181,521,1741,7621]
5^42 - 1 = [2,2,2,3,3,7,7,29,31,43,127,379,449,7603,19531,519499]
7^14 - 1 = [2,2,2,2,3,29,113,911,4733]
7^24 - 1 = [2,2,2,2,2,2,3,3,5,5,13,19,43,73,181,193,409,1201]
11^6 - 1 = [2,2,2,3,3,5,7,19,37]
11^12 - 1 = [2,2,2,2,3,3,5,7,13,19,37,61,1117]
11^24 - 1 = [2,2,2,2,2,3,3,5,7,13,19,37,61,1117,7321,10657,20113]

-}

galois.txt

To use these functions, you need:

1.  to have them both loadable (in the same directory)
2.  to have download and installed the arithmoi package: https://hackage.haskell.org/package/arithmoi
 
You define a field by giving its prime p , its power d, and an irreducible polynomial modulo p:
 
>> gf = GF 2 8 [1,0,0,0,1,1,1,0,1]
 
Polynomials as given as lists of their coefficients, starting with the highest power, so that the constant term is on the left.
Field parameters can be obtained as

>> basePrime gf
>> power gf
>> irreducible gf

Currently there are no tests to check whether p is in fact prime, or if the polynomial is irreducible.  You can create elements either by the list of coefficients:

>> x = GFElement gf [1,0,0,1,0,0,1,1]

or by choosing one at random:

>> y = randomGF gf 100

Note that the randomGF function requires a seed.  To obtain a different random polynomial you need a different seed.  Future versions will include either the state monad Control.Monad.State or the random monad Control.Monad.Random wihch will alleviate always having to provide a seed.

Arthmetic is provided by addGF, subGF, mulGF, expGF and their infix operator versions &+, &-, &*, &^:

>> x &+ y
>> x &* y
>> x &^ 50

Other functions are logGF, which finds the logarithm using a primitive element:

>> p = smallPrimitiveGF gf
>> logGF p x
>> q = randomPrimitiveGF gf 100
>> logGF q y

and ordGF which finds the multiplicative order:

>> ordGF x

None of these functions are particularly optimized, and many run very slowly on large fields.  For a field with p = 5, d = 42, computing logs took over 13 minutes!  There is plenty of room for improvement.

Other functions can be found in the two files Galois.hs and PolyUtils.hs

You can find irreducible polynomials (the Conway polynomials) at http://www.math.rwth-aachen.de/~Frank.Luebeck/data/ConwayPol/index.html
Note that the polynomials there are given in order of increasing powers; the lists of coefficients thus need to be reversed before using them here.

PolyUtils.hs

 module PolyUtils
   (addPoly
   ,subPoly
   ,mulPoly
   ,scalePoly
   ,rightTrim
   ,leftTrim
   ,degPoly
   ,divPoly
   ,remPoly
   ,gcdPoly
   ,expPoly
   ,toBase
   ,modInv
   ,primeFactors
   ,divisors
   ) where
  

{-
This is s small collection of routines for managing polynomials modulo a prime p.

Polynomials are treated as lists of coefficients of descending powers, so that

  9X^6 + 7X^3 + 4X^2 + 3X + 2

has the form

  [9,0,0,7,4,3,2]

To work with modular polynomials, use the ModPoly type with its constructor MP:

  p1 = MP [4,3,2,0,1]::ModPoly 7
  p2 = MP [3,2,1,4]::ModPoly 7

  p1 @+ p2      addition
  p1 @* p2      multiplication
  divp p1 p2    quotient

etc

-}

{-# LANGUAGE TypeOperators, DataKinds #-}
{-# LANGUAGE FlexibleContexts, TypeFamilies, UndecidableInstances #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}   -- needed for "forall"

import Data.List
import System.Random
import Math.NumberTheory.GCD (extendedGCD)
import Math.NumberTheory.ArithmeticFunctions (divisorsList)
import Math.NumberTheory.Primes.Factorisation (factorise)
import Math.NumberTheory.Moduli.Chinese (chineseRemainder)

primeFactors :: Integer -> [Integer]
primeFactors x = foldl (++) [] $  map (\z -> replicate (snd z) (fst z)) $ factorise x

divisors :: Integer -> [Integer]
divisors x = map (toInteger) $ divisorsList x

modInv :: Integer -> Integer -> Integer
modInv a b = toInteger $ (x `mod` b)
  where
    (g,x,y) = extendedGCD a b

-- toBase n k returns the list of "digits" of n in base k

toBase :: Integer -> Integer -> [Integer]
toBase n k
  | div n k == 0 = [mod n k]
  | otherwise    = (toBase (div n k) k)++[mod n k]

data MP = MP{modulus::Integer
            ,coeffs::[Integer]
            } deriving (Eq, Show)

-- adds lists starting from the right ends
rightAdd :: (Integral a) => [a] -> [a] -> [a]
rightAdd x [] = x
rightAdd [] y = y
rightAdd x y = (rightAdd (init x) (init y)) ++ [(last x)+(last y)]

addPoly :: MP -> MP -> MP
addPoly x y
  | modulus x /= modulus y = error "Can't add polynomials with different moduli"
  | otherwise = MP p $ map (`mod` p) (rightAdd xc yc)
  where
    p = modulus x
    xc = coeffs x
    yc = coeffs y

(@+) :: MP -> MP -> MP 
x @+ y  = addPoly x y 

-- subtracts lists starting from the right ends
rightSub :: (Integral a) => [a] -> [a] -> [a]
rightSub x [] = x
rightSub [] y = y
rightSub x y = (rightSub (init x) (init y)) ++ [(last x)-(last y)]

subPoly :: MP -> MP -> MP
subPoly x y 
  | modulus x /= modulus y = error "Can't add polynomials with different moduli"
  | otherwise = MP p $ map (`mod` p) (rightSub xc yc)
  where
    p = modulus x
    xc = coeffs x
    yc = coeffs y
    
(@-) :: MP -> MP -> MP
x @- y = subPoly x y

scalePoly :: Integer -> MP -> MP
scalePoly k x = MP p $ map (\z -> z*k `mod` p) (coeffs x)
  where p = modulus x

listMul :: (Integral a) => [a] -> [a] -> [a]
listMul x [] = []
listMul [] y = []
listMul x y = rightAdd (map (*(last y)) x) ((listMul (init y) x)++[0])

mulPoly :: MP -> MP -> MP
mulPoly x y
  | modulus x /= modulus y = error "Can't add polynomials with different moduli"
  | otherwise = MP p (map (`mod` p) (listMul xc yc))
  where
    p = modulus x
    xc = coeffs x
    yc = coeffs y
    
(@*) :: MP -> MP -> MP
x @* y  = mulPoly x y

-- leftTrim and rightTrim trim leading and trailing zeros respectively

leftTrim :: MP -> MP
leftTrim x = MP (modulus x) $ dropWhile (==0) (coeffs x)

rightTrim :: MP -> MP
rightTrim x = MP (modulus x) $ reverse $ dropWhile (==0) $ reverse (coeffs x)

-- degree of polynomial: zeros at the front are ignored

degPoly :: MP -> Int
degPoly x = if w <= 1 then 0
            else w-1
  where
    w = length $ coeffs $ leftTrim x

{-
-- quoRem :: Integral a => [a] -> [a] -> ([a],[a])
-- quoRem p q
--   | (degp q) < (degp p) = ([],q)
--   | otherwise = ([w]++ (fst $ quoRem $ p (tail q)),tail $ (q @- (scalep w p)))
--   where
--     w = div (head q) (head p)
-}

divp :: Integer -> [Integer] -> [Integer] -> [Integer]
divp p (x:xs) (y:ys)
  | (length (x:xs)) < (length (y:ys)) = []
  | ys == []                      = map (*(modInv y p)) (x:xs) 
  | otherwise                     = [w] ++ (divp p v (y:ys))
  where
    w = x * (modInv y p) `mod` p
    v = frontSub xs (map (*w) ys)

divPoly :: MP -> MP -> MP
divPoly x y 
  | modulus x /= modulus y = error "Can't add polynomials with different moduli"
  | otherwise = MP p $ map (`mod` p) (divp p xc yc)
  where
    p = modulus x
    xc = coeffs x
    yc = coeffs y
    
(@/) :: MP -> MP -> MP
x @/ y  = divPoly x y

-- frontSub is a helper function for divp; it subtracts lists starting from the left

frontSub :: Integral a => [a] -> [a] -> [a]
frontSub x [] = x
frontSub [] y = map (*(-1)) y
frontSub (x:xs) (y:ys) = [x-y] ++ (frontSub xs ys)

remp :: Integer -> [Integer] -> [Integer] -> [Integer]
remp p x y = dropWhile (==0) $ map (`mod` p) $ rightSub x (listMul y (divp p x y))

remPoly :: MP -> MP -> MP
remPoly x y
  | modulus x /= modulus y = error "Can't add polynomials with different moduli"
  | otherwise = MP p (map (`mod` p) (remp p xc yc))
  where
    p = modulus x
    xc = coeffs x
    yc = coeffs y
    
gcdp :: Integer -> [Integer] -> [Integer] -> [Integer]
gcdp p [] y = toMon y
gcdp p x [] = toMon x
gcdp p x y
  | x == y    = toMon x
  | otherwise = gcdp p y (remp p x y)

gcdPoly :: MP -> MP -> MP
gcdPoly x y
  | modulus x /= modulus y = error "Can't add polynomials with different moduli"
  | otherwise = MP p (map (`mod` p) (gcdp p xc yc))
  where
    p = modulus x
    xc = coeffs x
    yc = coeffs y
    
toMon :: Integral a => [a] -> [a]
toMon x = map (`div` head x) x

-- computes x^k mod y where x,y polynomials, k an Int
expRemP :: Integer -> [Integer] -> Integer -> [Integer] -> [Integer]
expRemP p x 0 y = [1]
expRemP p x k y
  | even k = remp p (listMul z z) y
  | otherwise = remp p (listMul x $ listMul z z) y
  where z = expRemP p x (div k 2) y

expPoly :: MP -> Integer -> MP -> MP
expPoly x k y
  | modulus x /= modulus y = error "Can't add polynomials with different moduli"
  | otherwise = MP p (map (`mod` p) (expRemP p xc k yc))
  where
    p = modulus x
    xc = coeffs x
    yc = coeffs y

-- isIrreducible checks a polynomial; for, say p = [1,2,3,4,5]::[Integer/7] you enter
-- isIrreducible p 7:: (KnownNat m)  => ModPoly m -> Bool

isIrreducible :: MP -> Bool
isIrreducible x
  | w /= x1 = False
  | any (\x->length (coeffs x)>1) ws = False
  | otherwise = True
  where
    x1 = MP p [1,0]
    p = modulus x
    d = toInteger $ degPoly x
    w = expPoly x1 (p^d) x
    fs = nub $ primeFactors d
    ws = map (\z -> gcdPoly (expPoly x1 (p^(div d z)) x @- x1) x) fs


p1 = MP 7 [5,4,3,2]
p2 = MP 7 [6,1,3,5]
p3 = MP 7 [1,2,3,4,3,2,1]

-- Now in highest power order (constant term at far right) x6/p =
-- [5,6,1,3] with remainder [4,2]

x6 = MP 7 [1,0,0,0,0,0,0]
x7 = MP 7 [1,0,0,0,0,0,0,0]
p4 = MP 7 [1,3,2,4]

-- x7/p4 = [1,4,0,2,6] with remainder [6,1,4]

z = MP 7 [1,0]

-- these next are all  irreducible

ir = MP 7 [1,4,1,2,6,4]
ir6 = MP 7 [1,4,1,1,1,6,4]
ir6b = MP 7 [1,2,5,5,2,4,5]
ir12 = MP 7 [1,0,0,0,2,5,3,2,4,0,5,0,3]