diff --git a/src/Classical/cyclic_code.jl b/src/Classical/cyclic_code.jl index b71307d..5dcb4b6 100644 --- a/src/Classical/cyclic_code.jl +++ b/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 # """ diff --git a/src/Classical/cyclotomic.jl b/src/Classical/cyclotomic.jl index 79c9a09..465c4cc 100644 --- a/src/Classical/cyclotomic.jl +++ b/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 \ No newline at end of file diff --git a/src/Classical/weight_dist.jl b/src/Classical/weight_dist.jl index 1e40c67..daa0aa5 100644 --- a/src/Classical/weight_dist.jl +++ b/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 diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index 621e998..de47291 100644 --- a/src/CodingTheory.jl +++ b/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 diff --git a/src/Quantum/GeneralizedToricCode.jl b/src/Quantum/GeneralizedToricCode.jl index be5312a..da28b86 100644 --- a/src/Quantum/GeneralizedToricCode.jl +++ b/src/Quantum/GeneralizedToricCode.jl @@ -34,12 +34,12 @@ end Return the bivariate bicycle code defined by the residue ring elements `a` and `b`. """ -function BivariateBicycleCode(a::T, b::T) where T <: Union{MPolyQuoRingElem{FqMPolyRingElem}, - MPolyQuoRingElem{fpMPolyRingElem}} +function BivariateBicycleCode(a::T, b::T) where T <: Union{MPolyQuoRingElem{FqMPolyRingElem}, MPolyQuoRingElem{fpMPolyRingElem}} + # pass in quo ring elems then make sure what goes to CSS are the lifted elems R = parent(a) R == parent(b) || throw(DomainError("Polynomials must have the same parent.")) - F = base_ring(base_ring(a)) + F = base_ring(base_ring(a)) # currenty MPolyQuoRing{fpMPolyRingElem} but constructor wants PolyRing{<:FinFieldElem} where FinFieldElem is AbstractAlgebra.FinFieldElem order(F) == 2 || throw(DomainError("This code family is currently only defined over binary fields.")) length(symbols(parent(a))) == 2 || throw(DomainError("Polynomials must be over two variables.")) g = gens(modulus(R)) @@ -68,7 +68,10 @@ function BivariateBicycleCode(a::T, b::T) where T <: Union{MPolyQuoRingElem{FqMP Q, _ = quo(R, I) k_dim = 2 * vector_space_dimension(Q) - return BivariateBicycleCode(R, F, 2 * l * m, k_dim, a, b, l, m) + R_base = base_ring(R) + a_lift = lift(a) + b_lift = lift(b) + return BivariateBicycleCode(R_base, F, 2 * l * m, k_dim, a_lift, b_lift, l, m) end """ @@ -276,7 +279,7 @@ function CSSCode(S::BivariateBicycleCode) y = identity_matrix(S.F, S.l) ⊗ matrix(S.F, [mod1(i + 1, S.m) == j ? 1 : 0 for i in 1:S.m, j in 1:S.m]) A = zero_matrix(S.F, S.l * S.m, S.l * S.m) - for ex in exponents(lift(S.a)) + for ex in exponents(lift(S.f)) # iszero(ex[1]) || iszero(ex[2]) || throw(ArgumentError("Polynomial `a` must not have any `xy` terms")) power, which = findmax(ex) if which == 1 @@ -287,7 +290,7 @@ function CSSCode(S::BivariateBicycleCode) end B = zero_matrix(S.F, S.l * S.m, S.l * S.m) - for ex in exponents(lift(S.b)) + for ex in exponents(lift(S.g)) # iszero(ex[1]) || iszero(ex[2]) || throw(ArgumentError("Polynomial `b` must not have any `xy` terms")) power, which = findmax(ex) if which == 1 @@ -297,7 +300,6 @@ function CSSCode(S::BivariateBicycleCode) end end - return CSSCode(hcat(A, B), hcat(transpose(B), transpose(A))) end diff --git a/src/Quantum/types.jl b/src/Quantum/types.jl index 9adbb9a..f653814 100644 --- a/src/Quantum/types.jl +++ b/src/Quantum/types.jl @@ -365,12 +365,12 @@ struct FiniteMonomialCode <: AbstractMonomialCode end struct BivariateBicycleCode <: AbstractBivariateBicycleCode - R::CTPolyRing + R::fpMPolyRing F::CTFieldTypes n::Int k::Int - f::CTPolyRingElem - g::CTPolyRingElem + f::fpMPolyRingElem + g::fpMPolyRingElem l::Int m::Int end diff --git a/test/Classical/cyclic_code_test.jl b/test/Classical/cyclic_code_test.jl index 312b567..966348b 100644 --- a/test/Classical/cyclic_code_test.jl +++ b/test/Classical/cyclic_code_test.jl @@ -1,4 +1,5 @@ -@testitem "Classical/cyclic_code.jl" begin +@testset "Classical/cyclic_code.jl" begin + using Test using Oscar, CodingTheory @testset "Cyclic codes" begin @@ -42,7 +43,9 @@ x = gen(R) @test generator_polynomial(C) == 2 + x + x^2 + x^3 @test dimension(C) == 10 - @test_broken minimum_distance(C) == 3 + @test minimum_distance(C) == 3 + # @test minimum_distance(C) == 3 + C = BCHCode(3, 13, 3, 1) @test defining_set(C) == [1, 2, 3, 5, 6, 9] @test generator_polynomial(C) == 1 + 2x + x^2 + 2x^3 + 2x^4 + 2x^5 + x^6