Cyclic

src/Classical/cyclic_code.jl

@@ -45,16 +45,26 @@ function CyclicCode(q::Int, n::Int, cosets::Vector{Vector{Int}})
 
     def_set = sort!(reduce(vcat, cosets))
     k = n - length(def_set)
-    com_cosets = complement_qcosets(q, n, cosets)
+    com_cosets = complement_qcosets(q, n, cosets) 
     g = _generator_polynomial(R, β, def_set)
-    h = _generator_polynomial(R, β, reduce(vcat, com_cosets))
-    e = _idempotent(g, h, n)
+    if isempty(com_cosets)
+        h = nothing
+        e = 0 
+    else
+        h = _generator_polynomial(R, β, reduce(vcat, com_cosets)) 
+        e = _idempotent(g, h, n)
+    end
+
+    if isempty(com_cosets)
+        return
+    end
     G = _generator_matrix(E, n, k, g)
     H = _generator_matrix(E, n, n - k, reverse(h))
     G_stand, H_stand, P, rnk = _standard_form(G)
     # HT will serve as a lower bound on the minimum weight
     # take the weight of g as an upper bound
-    δ, b, HT = find_delta(n, cosets)
+    δ, b = find_delta(n, cosets)
+    HT = -1
     ub = wt(G[1, :])
 
     # verify
@@ -240,7 +250,8 @@ function BCHCode(q::Int, n::Int, δ::Int, b::Int = 0)
     G_stand, H_stand, P, rnk = _standard_form(G)
     # HT will serve as a lower bound on the minimum weight
     # take the weight of g as an upper bound
-    δ, b, HT = find_delta(n, cosets)
+    δ, b = find_delta(n, cosets)
+    HT = -1
     upper = wt(G[1, :])
 
     # verify
@@ -613,13 +624,45 @@ end
 function _idempotent(g::FqPolyRingElem, h::FqPolyRingElem, n::Int)
     # solve 1 = a(x) g(x) + b(x) h(x) for a(x) then e(x) = a(x) g(x) mod x^n - 1
     d, a, b = gcdx(g, h)
-    return mod(g * a, gen(parent(g))^n - 1)
+    @assert d==1
+    return mod(g * a, gen(parent(g))^n - 1) 
 end
 
 # TODO: these
 # MattsonSolomontransform(f, n)
 # inverseMattsonSolomontransform
 
+function find_delta(n::Int, cosets::Vector{Vector{Int}})
+    n <= 0 && throw(ArgumentError("n must be positive"))
+
+    cosets = collect(Iterators.flatten(cosets))   
+    present = falses(n) # eg if [0] in cosets then present[1] = true
+    for e in cosets
+        present[mod(e, n) + 1] = true
+    end
+
+    curr_run_len = 0
+    best_run_len = curr_run_len
+    curr_offset = 1
+    best_offset = curr_offset 
+    for i in 1:(2n)
+        idx = mod(i - 1, n) + 1
+        if present[idx]
+            curr_run_len += 1
+            if curr_run_len > best_run_len
+                best_offset = curr_offset
+                best_run_len = min(curr_run_len, n)  
+            end
+        else
+            curr_offset = idx + 1
+            curr_run_len = 0
+        end
+    end
+    @assert curr_run_len < n
+    return best_run_len+1, best_offset
+end
+
+#=
 """
     find_delta(n::Int, cosets::Vector{Vector{Int}})
 
@@ -700,6 +743,7 @@ function find_delta(n::Int, cosets::Vector{Vector{Int}})
 
     return δ, offset, currbound
 end
+=#
 
 """
     dual_defining_set(def_set::Vector{Int}, n::Int)
@@ -880,6 +924,83 @@ Return `true` if the BCH code is antiprimitive.
 """
 is_antiprimitive(C::AbstractBCHCode) = C.n == Int(order(C.F)) + 1
 
+function print_all_cyclotomic_cosets(n::Int, q::Int)
+    println("All cyclotomic cosets for $(n) modulo $(q):")
+    rng = [i for i in 1:n]
+    flat=false
+    qcosets = defining_set(rng, q, n, flat) 
+    qcosets = unique(qcosets)
+    rows = Vector()
+    for i in 1:length(qcosets) 
+      cosets_one_removed = copy(qcosets)
+      deleteat!(cosets_one_removed, i);
+      Cd = CyclicCode(q, n, cosets_one_removed) 
+      d = dimension(Cd)
+      g = generator_polynomial(Cd)
+      e = idempotent(Cd)
+      push!(rows, [d,i,g,e,sort(cosets_one_removed)])
+    end
+    w1 = maximum(textwidth(string(r[2])) for r in rows) + 1
+    w2 = maximum(textwidth(string(r[3])) for r in rows) + 1
+    w3 = maximum(textwidth(string(r[4])) for r in rows) + 1
+    Printf.@printf("%-*s %-*s  %-*s   %s\n", w1, "i", w2, "g", w3, "e", "cosets")
+    for (a, b, c, d, e) in rows
+        Printf.@printf("%-*s %-*s  %-*s   %s\n", w1, b, w2, c, w3, d, e)
+    end
+end
+
+function print_all_cyclic_codes(n::Int, q::Int, def_set=false)
+    println("All cyclic codes of length $(n):")
+    rng = [i for i in 1:n]
+    flat=false
+    qcosets = defining_set(rng, q, n, flat) 
+    qcosets = unique(qcosets)
+
+    combs = sort(collect(Combinatorics.combinations(qcosets)))
+    rows = Vector()
+    for i in 1:length(combs) 
+        comb = combs[i]
+        if isempty(comb) || sort(collect(Iterators.flatten(comb))) == collect(0:n-1) # CyclicCode constructor breaks on this case
+            continue
+        end
+      C = CyclicCode(q, n, comb) 
+      g = generator_polynomial(C)
+      e = idempotent(C)
+      d = dimension(C)
+      if !def_set
+          push!(rows, [d,g,e])
+      else
+          dset = defining_set(C)
+          if C.δ > 2
+            push!(rows, [d,g,e,dset,repr(C.δ)])
+          else
+            push!(rows, [d,g,e,dset,"-"])
+          end
+      end
+    end
+    sort!(rows, by=first)
+    w1 = maximum(textwidth(string(r[1])) for r in rows) + 1
+    w2 = maximum(textwidth(string(r[2])) for r in rows)
+    if def_set
+        w3 = maximum(textwidth(string(r[3])) for r in rows)
+        w4 = maximum(textwidth(string(r[4])) for r in rows)
+    end
+    if !def_set
+      Printf.@printf("%*s   %-*s   %s\n", w1, "dim ", w2, "gen poly", "idempotent")
+    else
+      Printf.@printf("%*s   %-*s  %-*s  %-*s  %s\n", w1, "dim ", w2, "gen poly", w3, "idempotent", w4, "def set", "δ")
+    end
+    if !def_set
+      for (a, b, c) in rows
+          Printf.@printf("%*d   %-*s   %s\n", w1, a, w2, b, c)
+      end
+    else
+      for (a, b, c, d, e) in rows
+          Printf.@printf("%*d   %-*s  %-*s  %-*s  %s\n", w1, a, w2, b, w3, c, w4, d, e)
+      end
+    end
+end
+
 # "Schur products of linear codes: a study of parameters"
 # Diego Mirandola
 # """

src/Classical/cyclotomic.jl

@@ -209,3 +209,79 @@ function dual_qcosets(q::Int, n::Int, qcosets::Vector{Vector{Int64}})
     end
     return comp_cosets
 end
+
+function _coerce_Kx_to_Ky_x(p::PolyRingElem{T}, A) where {T <: RingElement}
+    # coerce p(x) in K[x] into A = Ky[x] (where Ky = K[y]) by coefficient embedding.
+    Kx = parent(p)
+    xA = gen(A)
+    coeffs = Oscar.coefficients(p)  
+    q = zero(A)
+    for i in 0:length(coeffs)
+        q += A(coeffs[i]) * xA^i
+    end
+    return q
+end
+
+function _composed_product(
+    f::PolyRingElem{T},
+    g::PolyRingElem{T}
+) where {T <: RingElement}
+    # computes the composed-product of two univariate polynomials using the resultant 
+
+    # the multivariate resultant called here only accepts AbstractAlgebra.Generic.Poly type as input
+    Kx = parent(f)
+    parent(g) === Kx || throw(ArgumentError("f and g must have the same parent K[x]."))
+    K = base_ring(Kx)
+    d = degree(g)
+    d < 0 && throw(ArgumentError("g must be nonzero."))
+    Kx, x = polynomial_ring(K, :x)
+    A,  y = polynomial_ring(Kx, :y)  # resultant(f, g) now eliminates y
+    fA = _coerce_Kx_to_Ky_x(f, A)
+    gA = _coerce_Kx_to_Ky_x(g, A)
+    fy = evaluate(fA, y)
+
+    fy = (fy + zero(A)) # converts the type of fy to AbstractAlgebra.Generic.Poly
+    # gy_scaled = x^d * g(y/x) = sum_{i=0}^d a_i * y^i * x^(d-i)
+    gy_scaled = zero(A)
+    for i in 0:d
+        ai = coeff(gA, i)           
+        gy_scaled += ai * y^i * x^(d - i)
+    end
+    res = resultant(fy, gy_scaled)
+    return res
+end
+
+"""
+    construct_field(f::PolyRingElem{T}, g::PolyRingElem{T}) where {T <: RingElement}
+
+constructs the smallest finite field extension where f and g split simultaneously
+"""
+function construct_field(f::PolyRingElem{T}, g::PolyRingElem{T}) where {T <: RingElement}
+    xfacs = [x[1] for x in factor(f)]
+    yfacs = [y[1] for y in factor(g)]
+    comp_prods = [_composed_product(f1,f2) for f1 in xfacs, f2 in yfacs]
+    comp_prods = [begin
+        facs = collect(factor(c))                
+        facs[argmax(degree(f[1]) for f in facs)][1]
+    end for c in comp_prods] # select irreducible factor of highest degree
+    dgs = unique([degree(c) for c in comp_prods])
+
+    m = 1 # lcm of the degrees
+    for dg in dgs
+        m = lcm(m, dg)
+    end
+    # |K| = 2^k where k is the smallest integer with 2^k=1 (mod lcm(d_i))
+    ZZm, _ = residue_ring(Nemo.ZZ, m) 
+    two = ZZm(2)
+    one = ZZm(1)
+    i = 0
+    for i in 0:m
+        if two^i == one
+            break
+        end
+    end
+    if i == m
+        throw(Error("failed to construct finite field"))
+    end
+    return Oscar.GF(2, m, :α), m, comp_prods
+end

src/Classical/weight_dist.jl

@@ -1723,3 +1723,60 @@ function minimum_distance(C::AbstractLinearCode; alg::Symbol = :trellis, sect::B
     #     Leon(C)
     end
 end
+
+# function verify_confinement(C::AbstractLinearCode, max_weight::Int, verbose::Bool = false)
+
+#     generator_matrix(C, true) # ensure G_stand exists
+#     if _has_empty_vec(C.G, :cols) 
+#         throw(ArgumentError("Codes with standard form of generator matrix having 0 columns not supported"))
+#     end
+#     # generate if not pre-stored
+#     parity_check_matrix(C)
+
+#     k, n = size(C.G)
+#     println(typeof(C.G))
+
+#     A = deepcopy(_Flint_matrix_to_Julia_int_matrix(C.G)') 
+#     A_mats_trunc = deepcopy(A[k + 1 : n, :]) 
+
+#     founds = [copy(zeros(UInt16, C.n - C.k)) for _ in 1:2^max_weight]
+#     count = UInt128(0)
+
+#     for r in 1:max_weight
+#         # iteration begins with a single matrix multiplication of the generator matrix by first_vec
+#         init_rank = 1 
+#         first_vec = zeros(Int, k)
+#         if init_rank == 1
+#             for i in 1:r
+#                 first_vec[i] = 1
+#             end
+#         else
+#             CodingTheory._subset_unrank_to_vec!(init_rank, UInt64(r), first_vec)
+#         end
+#         c_itr = zeros(UInt16, C.n - C.k)
+#         is_first = true
+
+#         for u in SubsetGrayCode(k, r, len, init_rank)
+#             show_progress && ProgressMeter.next!(prog_bar)
+#             if is_first 
+#                 LinearAlgebra.mul!(c_itr, A_mats_trunc, first_vec)
+#                 @inbounds @simd for j in eachindex(c_itr)
+#                     c_itr[j] %= p
+#                 end
+#                 is_first = false
+#             else
+#                 for ci in u 
+#                     if ci != -1
+#                         @simd for i in eachindex(c_itr)
+#                             @inbounds c_itr[i] = xor(c_itr[i], A_mats_trunc[i, ci])
+#                         end
+#                     end
+#                 end
+#             end
+#             w = r + sum(c_itr) 
+#             verbose && @assert w != 0
+#             count = add!(count, count, 1)
+#             founds[count] =  
+#         end # message loop
+#     end # weight loop
+# end

src/CodingTheory.jl

@@ -21,6 +21,8 @@ using Distributions
 using ProgressMeter
 using DocStringExtensions
 
+import Combinatorics
+import Printf 
 import LinearAlgebra: tr, Adjoint, transpose, kron, diagm
 import Oscar: dual, factor, transpose, order, polynomial, nrows, ncols, degree,
     lift, quo, vector_space, dimension, extend, support, complement,
@@ -54,8 +56,7 @@ const CTPolyRingElem = PolyRingElem{<:CTFieldElem}
 const CTGroupAlgebra = GroupAlgebraElem{fpFieldElem, GroupAlgebra{fpFieldElem, FinGenAbGroup, FinGenAbGroupElem}}
 const CTChainComplex = Union{ComplexOfMorphisms{AbstractAlgebra.FPModule{fpFieldElem}}} # residue and group algebras later
 const CTPolyMatrix = Union{AbstractAlgebra.Generic.MatSpaceElem{fpPolyRingElem}, AbstractAlgebra.Generic.MatSpaceElem{FqPolyRingElem}}
-const CTLRPolyElem = AbstractAlgebra.Generic.LaurentMPolyWrap{fpFieldElem, fpMPolyRingElem,
-AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}}
+const CTLRPolyElem = AbstractAlgebra.Generic.LaurentMPolyWrap{fpFieldElem, fpMPolyRingElem, AbstractAlgebra.Generic.LaurentMPolyWrapRing{fpFieldElem, fpMPolyRing}}
 
 include("Classical/types.jl")
 export AbstractCode, AbstractNonadditiveCode, AbstractNonlinearCode, AbstractAdditiveCode,
@@ -257,7 +258,8 @@ export defining_set, splitting_field, polynomial_ring, primitive_root, offset,
     design_distance, qcosets, qcosets_reps, generator_polynomial, parity_check_polynomial,
     idempotent, is_primitive, is_narrowsense, is_reversible, find_delta, dual_defining_set,
     CyclicCode, BCHCode, ReedSolomonCode, complement, ==, ∩, +, QuadraticResidueCode,
-    zeros, BCH_bound, is_degenerate, nonzeros, is_antiprimitive, FireCode
+    zeros, BCH_bound, is_degenerate, nonzeros, is_antiprimitive, FireCode, 
+    print_all_cyclotomic_cosets, print_all_cyclic_codes
 
 #############################
 # Classical/quasi-cyclic_code.jl