use Conway polynomials for Galois fields

Project.toml

@@ -3,7 +3,7 @@ uuid = "a6d4fa9c-9e0b-4795-89f3-f481b7b5e384"
 license = "MIT"
 authors = ["Klaus Crusius <[email protected]>"]
 repo = "https://github.com/KlausC/CommutativeRings.jl"
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -12,7 +12,7 @@ Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
-Polynomials = "1,2,3"
+Polynomials = "1,2,3,4"
 Primes = "0.5"
 julia = "1"
 

src/CommutativeRings.jl

@@ -10,7 +10,8 @@ export Ring, RingInt, FractionRing, QuotientRing, Polynomial
 export ZZ, QQ, ZZmod, Frac, Quotient, UnivariatePolynomial, MultivariatePolynomial
 export GaloisField, FSeries
 
-export PowerSeries, O, precision, absprecision
+export SpecialPowerSeries, PowerSeries, O, precision, absprecision
+export Conway
 
 export Hom, Ideal
 
@@ -78,6 +79,8 @@ include("linearalgebra.jl")
 include("rationalcanonical.jl")
 include("powerseries.jl")
 include("specialseries.jl")
+include("conway.jl")
 
 using .SpecialPowerSeries
+using .Conway
 end # module

src/conway.jl

@@ -0,0 +1,181 @@
+module Conway
+
+export conway, quasi_conway
+
+using Primes
+using ..CommutativeRings
+using CommutativeRings: intpower, _isprimitive, coeffs
+
+const CONWAY_POLYNOMIALS = Dict{Tuple{Int,Int},Any}()
+
+function store_conway!(row::Vector)
+    p, n, c = row
+    key = (p, n)
+    @assert c[n+1] == 1
+    c = c[1:n]
+    ix = findlast(!iszero, c)
+    if ix !== nothing && ix > 0
+        resize!(c, ix)
+    end
+    d = CONWAY_POLYNOMIALS
+    d[key] = c
+end
+
+function store_conway!(poly::P) where P<:UnivariatePolynomial{<:ZZmod}
+    c = value.(coeffs(poly))
+    p = characteristic(P)
+    map!(x -> x >= 0 ? x : p + x, c, c)
+    n = deg(poly)
+    store_conway!([p, n, c])
+end
+
+"""
+    conway(p, n)
+
+Return the conway-polynomial of degree `n` over `ZZ/p`, as far as available.
+
+For `n in (1,2)` the polynomials are calculated and memoized in the data cache.
+"""
+function conway(p::Integer, n::Integer, X::Symbol = :x)
+    isprime(p) || throw(ArgumentError("p is not prime ($p)"))
+    n > 0 || throw(ArgumentError("n is not positive ($n)"))
+
+    c = get(CONWAY_POLYNOMIALS, (p, n), missing)
+    if ismissing(c) && n <= 2
+        poly = quasi_conway(p, n, X)
+        store_conway!(poly)
+        return poly
+    end
+    ismissing(c) && return c
+    B = (ZZ/p)[X]
+    B(c) + monom(B, n)
+end
+
+"""
+    quasi_conway(p, m, X::Symbol, nr::Integer=1, factors=nothing)
+
+Return the first irreducible polynomial over `ZZ/p` of degree `m`.
+
+The polynomials are scanned in the canonical order for Conway polynomials.
+
+Variable symbol is `X`. If given, `factors` must be the prime factorization of `m`.
+If `nr > 0` is given, the `nr+1`^st of found polynomials is returned.
+"""
+function quasi_conway(p::Integer, m::Integer, X::Symbol = :x, nr = 0, factors = nothing)
+    Z = ZZ / p
+    P = Z[X]
+    fact = factors === nothing ? factor(intpower(p, m) - 1) : factors
+    mm = prod(fact)
+    x = monom(P)
+    g = generator(Z)
+    s = (-1)^(m - 1)
+    nx = max(nr, 0)
+    # find the next irreducible, for which x is primitive (drop first nr-1)
+    for poly in Monic(P, m - 1)
+        gen = (poly(-x) * x - g) * s # observation for all calculated Conway polynomials
+        if isirreducible(gen) && _isprimitive((x, gen), mm, fact)
+            nx == 0 && return gen
+            nx -= 1
+        end
+    end
+    text = "no irreducible polynomial of degree $m found with generator '$X' (nr = $nr)"
+    throw(ArgumentError(text))
+end
+
+"""
+    has_conway_property1(poly)
+
+Check if `poly` is monomial.
+"""
+function has_conway_property1(poly::T) where {P,S<:ZZmod{P},T<:UnivariatePolynomial{S}}
+    ismonic(poly)
+end
+
+"""
+    has_conway_property2(poly)
+
+Check if an irreducible polynomial `poly` over `ZZ/p` of degree `n` has property 2 (is primitive).
+
+That means that the standard generator `x` in the quotient ring (field) `(ZZ/p) / poly`
+is also generator in the multiplicative group of the ring.
+
+If `m` is the order of this group (`p^n - 1` if a field), it sufficient to check, if
+`x ^ (m / s)` != 1 for all prime factors of `m`. (`y ^ m == 1` for all `y != 0`)
+
+The function returns `missing` if `p^n - 1` has no known prime factors.
+"""
+function has_conway_property2(
+    poly::T,
+    factors = nothing,
+) where {P,S<:ZZmod{P},T<:UnivariatePolynomial{S}}
+    p = P
+    fact = factors === nothing ? Primes.factor(intpower(p, m) - 1) : factors
+    mm = prod(fact)
+    isirreducible(poly) || throw(ArgumentError("polynomial is not irreducible"))
+    x = monom(T)
+    _isprimitive((x, poly), mm, fact)
+end
+
+"""
+    has_conway_property3(poly::(ZZ/p)[:x]})
+
+Check if a polynomial `poly` over `ZZ/p` of degree `n` has property 3.
+
+That means, for all divisors `1 <= m < n` we have
+`conway(p, m)(x^w) == 0 mod poly` with `w = (p^n - 1) / (p^m - 1)`.
+
+The definition is recursive and requires the knowledge of Conway polynomials of lesser degree.
+
+If such polynomials are unknown, `missing` is returned.
+"""
+function has_conway_property3(poly::T) where {P,S<:ZZmod{P},T<:UnivariatePolynomial{S}}
+    n = deg(poly)
+    p = P
+    Q = T / poly
+    for m in factors(n)
+        if m < n
+            nm = n ÷ m
+            q = big(p)^m
+            # w = (p^n - 1) / (p^m - 1); Z = x^w mod poly
+            X = monom(Q)
+            Y = X^q
+            Z = X * Y
+            for i = 2:nm-1
+                Y = Y^q
+                Z *= Y
+            end
+            w = (q^nm - 1) ÷ (q - 1)
+            @assert X^w == Z
+            c = conway(p, m)
+            ismissing(c) && return missing
+            !iszero(c(Z)) && return false
+        end
+    end
+    return true
+end
+
+function readparse(fn::AbstractString, action::Function)
+    open(fn, "r") do io
+        readparse(io, action)
+    end
+end
+function readparse(io::IO, action::Function)
+    while !eof(io)
+        line = readline(io)
+        if line[1] == '[' && line[end] == ','
+            row = eval(Meta.parse(line[1:end-1]))
+            action(row)
+        end
+    end
+end
+
+function __init__()
+    # file imported from
+    # `https://www.math.rwth-aachen.de/~Frank.Luebeck/data/ConwayPol/CPimport.txt.gz`
+    # at 2023-11-11
+    file = joinpath(@__DIR__, "..", "data", "CPimport.txt")
+    empty!(CONWAY_POLYNOMIALS)
+    readparse(file, r -> store_conway!(r))
+end
+
+end # module

src/galoisfields.jl

@@ -25,11 +25,12 @@ function GF(n::Integer, k::Integer = 1; mod = nothing, nr = 0, maxord = 2^20)
     f = Primes.factor(n)
     length(f) == 1 || throw(ArgumentError("$n is not p^r with p prime and r >= 1"))
     p, r = f.pe[1]
-    _GF(p, r * k; mod, nr, maxord)
+    _GF(p, r * k, mod, nr, maxord)
 end
-function _GF(p::Integer, r::Integer; nr::Integer = 0, mod = nothing, maxord = 2^20)
-    r == 1 || mod === nothing || throw(ArgumentError("given modulus requires prime base"))
-    mm = intpower(p, mod === nothing ? r : deg(mod)) - 1
+function _GF(p::Integer, r::Integer, mod, nr::Integer, maxord::Int)
+    modpol = mod isa UnivariatePolynomial
+    r == 1 || !modpol || throw(ArgumentError("given modulus requires prime base"))
+    mm = intpower(p, !modpol ? r : deg(mod)) - 1
     fact = Primes.factor(mm)
     Q = GFImpl(p, r, fact; nr, mod)
     ord = order(Q)
@@ -250,7 +251,8 @@ end
 Return the generator element of this implementation of Galois field.
 """
 function generator(::Type{G}) where {Id,G<:GaloisField{Id,<:Integer}}
-    G(2, NOCHECK)
+    id = min(2, order(G) - 1)
+    G(id, NOCHECK)
 end
 function generator(::Type{G}) where {Id,G<:GaloisField{Id,<:Quotient}}
     G(monom(Quotient(G)))
@@ -362,12 +364,12 @@ function tovalue(::Type{<:GaloisField{Id,V}}, num::Integer) where {Id,V<:Quotien
 end
 
 function tonumber(a::Quotient{<:UnivariatePolynomial}, p::Integer)
-    s = 0
     u = a.val
+    s = zero(intpower(p, deg(u)+1))
     for c in reverse(u.coeff)
         s = s * p + c.val
     end
-    s * p^ord(u)
+    s * intpower(p, ord(u))
 end
 
 function toquotient(g::G) where {Id,T<:Integer,Q,G<:GaloisField{Id,T,Q}}
@@ -382,7 +384,7 @@ end
 function toquotient(
     a::Integer,
     ::Type{Q},
-) where {Z,P<:UnivariatePolynomial{Z,:α},Q<:Quotient{P}}
+) where {Z,P<:UnivariatePolynomial{Z},Q<:Quotient{P}}
     p = characteristic(Q)
     r = dimension(Q)
     ord = order(Q)
@@ -429,34 +431,32 @@ function GFImpl(
 )
     isprime(p) || throw(ArgumentError("base $p must be prime"))
     m > 0 || throw(ArgumentError("exponent m=$m must be positive"))
-    Z = ZZ / p
-    m == 1 && mod === nothing && return Z
-    P = Z[:α]
-    if mod === nothing
-        fact = factors === nothing ? Primes.factor(intpower(p, m) - 1) : factors
-        mm = prod(fact)
-        x = monom(P)
-        nx = max(nr, 0)
-        # find the next irreducible, for which x is primitive (drop first nr-1)
-        for gen in irreducibles(P, m)
-            if _isprimitive((x, gen), mm, fact)
-                nx == 0 && return P / gen
-                nx -= 1
-            end
+
+    m == 1 && mod === nothing && return ZZ / p
+    if mod === :conway
+        gen = Conway.conway(p, m, :γ)
+        if !ismissing(gen)
+            return typeof(gen) / gen
         end
-        throw(
-            ArgumentError(
-                "no irreducible polynomial of degree $m found with generator p(γ) = γ (nr = $nr)",
-            ),
-        )
-    else
+        mod = nothing
+    end
+
+    if isnothing(mod)
+        poly = Conway.quasi_conway(p, m, :α, nr, factors)
+        return typeof(poly) / poly
+
+    elseif mod isa UnivariatePolynomial
         m == 1 || throw(ArgumentError("given mod requires prime base"))
+        Z = ZZ / p
+        P = Z[:α]
         # do not check if x is primitive here
         gen = P(Z.(mod.coeff), ord(mod))
         if isirreducible(gen)
             return P / gen
         end
-        throw(ArgumentError("given polynomial $gen is not irreducible in $P"))
+        throw(ArgumentError("given polynomial $gen is not irreducible over $Z"))
+    else
+        throw(ArgumentError("modulus '$mod' not supported"))
     end
 end
 
@@ -465,9 +465,9 @@ function Base.show(io::IO, g::G) where G<:GaloisField
     m = dimension(G)
     p = characteristic(G)
     cc = toquotient(g).val
-    print(io, '{', cc[m-1])
+    print(io, '{', cc[m-1].val)
     for k = m-2:-1:0
-        print(io, ':', cc[k])
+        print(io, ':', cc[k].val)
     end
     print(io, '%', p, '}')
 end
@@ -563,7 +563,7 @@ function *(
 end
 
 function monom(::Type{Q}, k::Integer, lc = 1) where {P<:UnivariatePolynomial,Q<:Quotient{P}}
-    Q(monom(P, k, lc))
+    k < dimension(Q) ? Q(monom(P, k, lc)) : lc == 1 ? Q(monom(P))^k : Q(monom(P))^k * lc
 end
 
 """

src/intfactorization.jl

@@ -370,7 +370,7 @@ calculate array of bounds `b` with `abs(v[i]) <= b[i+1] for i = 0:m`.
 Algorithm see TAoCP 2.Ed 4.6.2 Exercise 20.
 """
 function coeffbounds(u::UnivariatePolynomial{ZZ{T},X}, m::Integer) where {T<:Integer,X}
-    I = widen(T)
+    W = widen(T)
     n = deg(u)
     0 <= m <= n || throw(ArgumentError("required m ∈ [0,deg(u)] but $m ∉ [0,$n]"))
     accuracy = 100 # use fixed point decimal arithmetic with accuracy 0.01 for the norm
@@ -379,15 +379,15 @@ function coeffbounds(u::UnivariatePolynomial{ZZ{T},X}, m::Integer) where {T<:Int
     iszero(u0) && throw(ArgumentError("required u(0) != 0"))
     rua = norm(value.(u.coeff)) * accuracy
     ua = Integer(ceil(rua))
-    v = zeros(I, m + 1)
+    v = zeros(W, m + 1)
     if m > 0
         v[m+1], v[1] = un, u0
     else
         v[1] = min(u0, un) # abs(content(u)) would be possible but not worth the effort
     end
     u0 *= accuracy
     un *= accuracy
-    bk = I(1)
+    bk = W(1)
     for j = m-1:-1:1
         bj = bk
         bk = bk * j ÷ (m - j)

src/typevars.jl

@@ -64,5 +64,5 @@ sintern(m::IdentSymbols) = m
 sintern(m::BigInt) = Symbol(m)
 sintern(a::IdentSymbols...) = Symbol(tuple(a...))
 sintern(a::AbstractVector{<:IdentSymbols}) = Symbol(a)
-sintern(a::Tuple{Vararg{<:IdentSymbols}}) = Symbol(a)
+sintern(a::Tuple{Vararg{T}}) where T<:IdentSymbols = Symbol(a)
 sintern(a::Any) = Symbol(Base.hash(a))

src/zzmod.jl

@@ -76,7 +76,7 @@ function value(a::ZZmod{X,T}) where {X,T}
     m = modulus(a)
     m1 = m >> 1
     m2 = m - m1
-    v <= m2 ? S(v) : S(v - m2) - S(m1)
+    v < m2 ? S(v) : S(v - m2) - S(m1)
 end
 
 # get type variable
@@ -166,7 +166,7 @@ _unsigned(x::Integer) = unsigned(x)
 function Base.show(io::IO, a::ZZmod)
     v = a.val
     m = modulus(a)
-    if m >= 100 && v > m ÷ 2
+    if m > 10 && v > m ÷ 2
         print(io, '-', signed(m - v), '°')
     else
         print(io, signed(v), '°')