From 15d35ae6bfd68774a6f9afd12976c24c848ef504 Mon Sep 17 00:00:00 2001 From: Jake Ireland Date: Tue, 5 Nov 2024 20:28:55 +1300 Subject: [PATCH] Apply formatting suggestions --- .JuliaFormatter.toml | 26 ++ dev/.JuliaFormatter.toml | 0 docs/make.jl | 10 +- examples/bounds.jl | 128 +++--- examples/get_codewords_time_analysis.jl | 48 ++- examples/plot_random.jl | 18 +- examples/plot_random_core.jl | 285 +++++++------- examples/plot_random_from_file.jl | 22 +- examples/random_numbers_core.jl | 43 ++- examples/sort_bounds.jl | 170 ++++---- examples/turtle.jl | 136 ++++--- justfile | 21 +- perf/runbenchmarks.jl | 222 +++++++++-- src/CodingTheory.jl | 78 +++- src/abstract_types.jl | 79 ++-- src/algebra.jl | 151 ++++---- src/bounds.jl | 251 ++++++------ src/distance.jl | 192 +++++---- src/levenshtein.jl | 37 +- src/messages.jl | 452 ++++++++++++---------- src/powers.jl | 2 +- src/rref.jl | 101 +++-- src/utils.jl | 112 +++--- src/wip/primepowers.jl | 336 ++++++++-------- src/wip/primes.jl | 363 +++++++++-------- test/runtests.jl | 492 +++++++++++++++++------- 26 files changed, 2186 insertions(+), 1589 deletions(-) create mode 100644 .JuliaFormatter.toml delete mode 100644 dev/.JuliaFormatter.toml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..3d80988 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,26 @@ +indent = 4 +margin = 92 +always_for_in = true +for_in_replacement = "in" +whitespace_typedefs = true +whitespace_ops_in_indices = true +remove_extra_newlines = true +import_to_using = true +pipe_to_function_call = true +short_to_long_function_def = false +force_long_function_def = false +long_to_short_function_def = false +always_use_return = true +whitespace_in_kwargs = true +annotate_untyped_fields_with_any = false +format_docstrings = false +conditional_to_if = false +normalize_line_endings = "unix" +trailing_comma = true +trailing_zero = true # TODO +join_lines_based_on_source = false +indent_submodule = false +separate_kwargs_with_semicolon = false +surround_whereop_typeparameters = true +short_circuit_to_if = false +disallow_single_arg_nesting = true diff --git a/dev/.JuliaFormatter.toml b/dev/.JuliaFormatter.toml deleted file mode 100644 index e69de29..0000000 diff --git a/docs/make.jl b/docs/make.jl index 8375201..0dc42e8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,6 @@ using CodingTheory, Documenter -Documenter.makedocs( +Documenter.makedocs(; clean = true, doctest = true, modules = Module[CodingTheory], @@ -8,14 +8,10 @@ Documenter.makedocs( highlightsig = true, sitename = "CodingTheory Documentation", expandfirst = [], - pages = [ - "Index" => "index.md", - ] + pages = ["Index" => "index.md"], ) -deploydocs(; - repo = "github.com/jakewilliami/CodingTheory.jl.git", -) +deploydocs(; repo = "github.com/jakewilliami/CodingTheory.jl.git") # deploydocs( # target = "build", diff --git a/examples/bounds.jl b/examples/bounds.jl index b3d7a52..b62f378 100755 --- a/examples/bounds.jl +++ b/examples/bounds.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia -tauto --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - using CSV, DataFrames, StatsPlots, ProgressMeter include(joinpath(dirname(dirname(@__FILE__)), "src", "CodingTheory.jl")) @@ -16,29 +10,29 @@ function integer_search(stop_at::Integer)::Array{Array{Number, 1}} A = [] upper_bound = 10 increment_bound = 10 - processed = [] - i = 0 + processed = [] + i = 0 p = Progress(stop_at, 1) # minimum update interval: 1 second - + while true for q in 1:upper_bound, n in 1:upper_bound, d in 1:upper_bound - # skip configurations that have already been processed - # arelessthan(upper_bound - increment_bound, q, n, d) && continue - (q, n, d) ∈ processed && continue + # skip configurations that have already been processed + # arelessthan(upper_bound - increment_bound, q, n, d) && continue + (q, n, d) ∈ processed && continue # note that we actually don't want to filter out non-linear codes - # distance shouldn't be larger than the block length; filter trivial distances of one; we filter out combinations when all are equal; filter out codes that are generalised hamming codes; filter out even distances - # if d < n && ! isone(q) #&& ! CodingTheory.allequal(q, n, d) && ! ishammingperfect(q, n, d) && ! iseven(d) && ! isprimepower(q) - q, n, d = big(q), big(n), big(d) - hb = hamming_bound(q, n, d, no_round) - sb = singleton_bound(q, n, d, no_round) - - # if isinteger(hb) #&& ! isone(hb) && ! isone(sb) && hb ≤ sb - push!(processed, (q, n, d)) - # i += 1; println("$i:\t$q, $n, $d") - push!(A, [q, n, d, hb, sb]) - isequal(length(A), stop_at) && return A - next!(p) - # end + # distance shouldn't be larger than the block length; filter trivial distances of one; we filter out combinations when all are equal; filter out codes that are generalised hamming codes; filter out even distances + # if d < n && ! isone(q) #&& ! CodingTheory.allequal(q, n, d) && ! ishammingperfect(q, n, d) && ! iseven(d) && ! isprimepower(q) + q, n, d = big(q), big(n), big(d) + hb = hamming_bound(q, n, d, no_round) + sb = singleton_bound(q, n, d, no_round) + + # if isinteger(hb) #&& ! isone(hb) && ! isone(sb) && hb ≤ sb + push!(processed, (q, n, d)) + # i += 1; println("$i:\t$q, $n, $d") + push!(A, [q, n, d, hb, sb]) + isequal(length(A), stop_at) && return A + next!(p) + # end # end end upper_bound += increment_bound @@ -47,41 +41,52 @@ end function make_integer_csv(stop_at::Integer) A = integer_search(stop_at) - D = DataFrame(q = Number[], n = Number[], d = Number[], hamming_bound = Number[], singleton_bound = Number[]) + D = DataFrame(; + q = Number[], + n = Number[], + d = Number[], + hamming_bound = Number[], + singleton_bound = Number[], + ) for i in A push!(D, i) end - - data_file_desktop = joinpath(homedir(), "Desktop", "bound_integers_$(stop_at).csv") - data_file_other = joinpath(dirname(dirname(@__FILE__)), "other", "bound_integers_$(stop_at).csv") + + data_file_desktop = joinpath(homedir(), "Desktop", "bound_integers_$(stop_at).csv") + data_file_other = joinpath( + dirname(dirname(@__FILE__)), "other", "bound_integers_$(stop_at).csv" + ) CSV.write(data_file_desktop, D) CSV.write(data_file_other, D) - println("\nWrote data to $(data_file_other).") + return println("\nWrote data to $(data_file_other).") end function bound_comparison(stop_at::Integer)::Tuple{Int, Array{Array{Number, 1}}} A = Array{AbstractFloat}[] - processed = [] + processed = [] upper_bound = 10 increment_bound = 10 - i = 0 - + i = 0 + while true @showprogress for q in 1:upper_bound, n in 1:upper_bound, d in 1:upper_bound - # skip configurations that have already been processed - # arelessthan(upper_bound - increment_bound, q, n, d) && continue - (q, n, d) ∈ processed && continue - # remove q non prime powers + # skip configurations that have already been processed + # arelessthan(upper_bound - increment_bound, q, n, d) && continue + (q, n, d) ∈ processed && continue + # remove q non prime powers if isprimepower(big(q)) hb = hamming_bound(q, n, d) sb = singleton_bound(q, n, d) - # filter out the more trivial bounds that are one (or if distance is one); filter out distances that are more than or equal to the block length, as they don't make sense in our context; filter out even distances - # if d < n || 1 ∉ (hb, sb, d) - if d < n && ! isone(d) && ! CodingTheory.allequal(q, n, d) && ! iszero(mod(d, 2)) - push!(processed, (q, n, d)) - # i += 1; println("$i:\t$q, $n, $d") + # filter out the more trivial bounds that are one (or if distance is one); filter out distances that are more than or equal to the block length, as they don't make sense in our context; filter out even distances + # if d < n || 1 ∉ (hb, sb, d) + if d < n && + !isone(d) && + !CodingTheory.allequal(q, n, d) && + !iszero(mod(d, 2)) + push!(processed, (q, n, d)) + # i += 1; println("$i:\t$q, $n, $d") comparison = hb > sb ? 1 : (sb > hb ? -1 : 0) push!(A, [q, n, d, hb, sb, comparison]) isequal(length(A), stop_at) && return stop_at, A @@ -93,27 +98,32 @@ function bound_comparison(stop_at::Integer)::Tuple{Int, Array{Array{Number, 1}}} end function plot_bound_comparison(stop_at::Integer) - number_of_bounds, A = bound_comparison(stop_at) - D = DataFrame(q = Number[], n = Number[], d = Number[], hamming_bound = Number[], singleton_bound = Number[], comparison = Number[]) - - for i in A - push!(D, i) - end - - data_file_desktop = joinpath(homedir(), "Desktop", "bound_comparison_$(stop_at).csv") - data_file_other = joinpath(dirname(dirname(@__FILE__)), "other", "bound_comparison_$(stop_at).csv") - - CSV.write(data_file_desktop, D) - CSV.write(data_file_other, D) - println("Data written to $(data_file_desktop) and $(data_file_other)") - - `Rscript --vanilla "~"/projects/CodingTheory.jl/other/regression-tree.R $data_file_other $stop_at` - println("Regression tree saved at $(joinpath(dirname(dirname(@__FILE__)), "other", "Rplots_$(stop_at).pdf"))") -end + number_of_bounds, A = bound_comparison(stop_at) + D = DataFrame(; + q = Number[], + n = Number[], + d = Number[], + hamming_bound = Number[], + singleton_bound = Number[], + comparison = Number[], + ) + for i in A + push!(D, i) + end + data_file_desktop = joinpath(homedir(), "Desktop", "bound_comparison_$(stop_at).csv") + data_file_other = joinpath( + dirname(dirname(@__FILE__)), "other", "bound_comparison_$(stop_at).csv" + ) + CSV.write(data_file_desktop, D) + CSV.write(data_file_other, D) + println("Data written to $(data_file_desktop) and $(data_file_other)") + `Rscript --vanilla "~"/projects/CodingTheory.jl/other/regression-tree.R $data_file_other $stop_at` + return println("Regression tree saved at $(joinpath(dirname(dirname(@__FILE__)), "other", "Rplots_$(stop_at).pdf"))",) +end # make_integer_csv constructs a csv file which looks for hamming bounds (given certain conditions) that are integers before rounding make_integer_csv(stop_at) # takes a command line argument diff --git a/examples/get_codewords_time_analysis.jl b/examples/get_codewords_time_analysis.jl index 1e73373..29c4fae 100755 --- a/examples/get_codewords_time_analysis.jl +++ b/examples/get_codewords_time_analysis.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - using StatsPlots, ProgressMeter, CSV, DataFrames include(joinpath(dirname(dirname(@__FILE__)), "src", "CodingTheory.jl")) @@ -16,7 +10,7 @@ function graphing_time(stop_at::Integer) stop_n_at = 8 num_of_datapoints = stop_at data = Matrix{Float64}(undef, num_of_datapoints, 3) # needs 3 columns: q, n, and time - + @showprogress for i in 1:num_of_datapoints q = rand(1:stop_q_at) n = rand(1:stop_n_at) @@ -24,26 +18,46 @@ function graphing_time(stop_at::Integer) data[i, :] .= [q, n, time_took] end - save_path = joinpath(dirname(@__DIR__), "other", "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.pdf") + save_path = joinpath( + dirname(@__DIR__), + "other", + "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.pdf", + ) theme(:solarized) - - plot_points = scatter(data[:,1:2], data[:,3], fontfamily = font("Times"), ylabel = "Time elapsed during calculation [seconds]", clabel = "Value of q or n", label = ["q" "n"], smooth = true)#, xlims = (0, stop_n_at)) # smooth = true - + + plot_points = scatter( + data[:, 1:2], + data[:, 3]; + fontfamily = font("Times"), + ylabel = "Time elapsed during calculation [seconds]", + clabel = "Value of q or n", + label = ["q" "n"], + smooth = true, + )#, xlims = (0, stop_n_at)) # smooth = true + savefig(plot_points, save_path) println("Plot saved at $save_path.") - - D = DataFrame(q = Number[], n = Number[], time = Number[]) + + D = DataFrame(; q = Number[], n = Number[], time = Number[]) for i in eachrow(data) push!(D, i) end - - data_file_desktop = joinpath(homedir(), "Desktop", "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.csv") - data_file_other = joinpath(dirname(dirname(@__FILE__)), "other", "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.csv") + + data_file_desktop = joinpath( + homedir(), + "Desktop", + "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.csv", + ) + data_file_other = joinpath( + dirname(dirname(@__FILE__)), + "other", + "get_codewords_time_analysis,n=$stop_n_at,m=$stop_q_at,i=$num_of_datapoints.csv", + ) CSV.write(data_file_desktop, D) CSV.write(data_file_other, D) - println("Wrote data to $(data_file_other).") + return println("Wrote data to $(data_file_other).") end graphing_time(stop_at) diff --git a/examples/plot_random.jl b/examples/plot_random.jl index 422d674..c8ec839 100755 --- a/examples/plot_random.jl +++ b/examples/plot_random.jl @@ -1,17 +1,11 @@ -#!/usr/bin/env bash - #= - exec julia -tauto --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - include("plot_random_core.jl") if isempty(ARGS) || length(ARGS) < 4 - throw(error(""" - Please specify the number of datapoints you want to calculate from the command line, along with q, n, and d. For example - ./examples/plot_random.jl 6 7 3 1000 10000 # q, n, d, num_data_points, selection_set_size - That is respectively: the size of the alphabet; the block length of the words; the distance of the code; the number of datapoints in the code; and the selection set size of each "random" sample. - """)) + throw(error(""" + Please specify the number of datapoints you want to calculate from the command line, along with q, n, and d. For example + ./examples/plot_random.jl 6 7 3 1000 10000 # q, n, d, num_data_points, selection_set_size + That is respectively: the size of the alphabet; the block length of the words; the distance of the code; the number of datapoints in the code; and the selection set size of each "random" sample. + """,),) end global q = parse(Int, ARGS[1]) @@ -24,4 +18,4 @@ global upper_bound_adjustment = 10 global 🐖 = length(get_codewords_greedy(q, n, d)) global 🍖 = hamming_bound(q, n, d) global 🐺 = singleton_bound(q, n, d) -graphing(q, n, d, stop_at, m = m) +graphing(q, n, d, stop_at; m = m) diff --git a/examples/plot_random_core.jl b/examples/plot_random_core.jl index f171610..bd018cd 100755 --- a/examples/plot_random_core.jl +++ b/examples/plot_random_core.jl @@ -1,147 +1,174 @@ -#!/usr/bin/env bash - #= - exec julia -tauto --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - ENV["GKS_ENCODING"] = "utf-8" include("random_numbers_core.jl") -function add_plot_info(plt::AbstractPlot; extra::AbstractString=string(), counter::Integer=0) +function add_plot_info( + plt::AbstractPlot; extra::AbstractString = string(), counter::Integer = 0 +) x = xlims(Plots.current())[1] y = ylims(Plots.current())[2] factor1 = y > 1 ? 2 : 1.6 factor2 = y > 1 ? 4 : 2.5 # plot!(annotations = (max(x * 2, x+100), (y - 10 * counter) / factor1, # text(extra, :left, 13))) - annotate!(plt, max(x * 2, x + 100), ((y - 10 * counter) / factor1), extra) + return annotate!(plt, max(x * 2, x + 100), ((y - 10 * counter) / factor1), extra) end -function graphing(q::Integer, n::Integer, d::Integer, data::Union{AbstractArray, DataFrame}, stop_at::Integer; m::Integer=10_000) - save_path = joinpath(dirname(@__DIR__), "other", "random_number_analysis", "random_number_in_code,q=$(q),n=$(n),d=$(d),i=$(stop_at),m=$(m).pdf") - - if data isa DataFrame - data = Matrix(data) - end - - 🐖_obtained = ifelse(🐖 ∈ data, true, false) - - ## Plot points - theme(:solarized) - - x_max = 0 - if 🍖 - maximum(data) ∈ [0:upper_bound_adjustment...] - x_max = 🍖 - end - if 🐺 - maximum(data) ∈ [0:upper_bound_adjustment...] - x_max = 🐺 - end - - bin_adjustment = ifelse(🐖_obtained, 0.5, minimum(abs.(extrema(vcat(data, x_max)) .- 🐖)) + 0.5) # ensure the bins encompass the greedy - bins = minimum(data) - bin_adjustment : max(maximum(data), x_max) + bin_adjustment - +function graphing( + q::Integer, + n::Integer, + d::Integer, + data::Union{AbstractArray, DataFrame}, + stop_at::Integer; + m::Integer = 10_000, +) + save_path = joinpath( + dirname(@__DIR__), + "other", + "random_number_analysis", + "random_number_in_code,q=$(q),n=$(n),d=$(d),i=$(stop_at),m=$(m).pdf", + ) + + if data isa DataFrame + data = Matrix(data) + end + + 🐖_obtained = ifelse(🐖 ∈ data, true, false) + + ## Plot points + theme(:solarized) + + x_max = 0 + if 🍖 - maximum(data) ∈ [0:upper_bound_adjustment...] + x_max = 🍖 + end + if 🐺 - maximum(data) ∈ [0:upper_bound_adjustment...] + x_max = 🐺 + end + + bin_adjustment = ifelse( + 🐖_obtained, 0.5, minimum(abs.(extrema(vcat(data, x_max)) .- 🐖)) + 0.5 + ) # ensure the bins encompass the greedy + bins = (minimum(data) - bin_adjustment):(max(maximum(data), x_max) + bin_adjustment) + plt = histogram( - data, - fontfamily = font("Times"), - xlabel = "|C|", - ylabel = "Frequency of Times Found", - label = "", - title = "(q, n, d) = ($q, $n, $d): $(format(stop_at, commas = true)) Random Samples, with Selection Size $(format(m, commas = true))", - # xticks = minimum(data) : max(maximum(data), x_max), - bins = bins - ) - - ## Annotate plot - annotation_counter = 0 - if isempty(data) - add_plot_info(plt, extra="🐖 = $(🐖)", counter=annotation_counter); annotation_counter += 1 - add_plot_info(plt, extra="🍖 = $(🍖)", counter=annotation_counter); annotation_counter += 1 - add_plot_info(plt, extra="🐺 = $(🐺)", counter=annotation_counter); annotation_counter += 1 - savefig(plt, save_path) - println("Plot saved at $save_path.") - return nothing - end - - 🐖_idx = searchsortedlast(bins, 🐖) # <- the index of that value - 🐖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐖_idx] # height of bar of choice - annotate!(plt, 🐖, 🐖_y_top, "🐖") - - annotation_factor = 10 - x_annotation_factor = abs(bins[end] - bins[1]) * 0.05 - y_annotation_factor = abs(bins[end] - bins[1]) * 1.2 - 🍖_y, 🐺_y = 0, 0 - 🍖_annotated, 🐺_annotated = false, false - - if abs(x_max - 🍖) ≤ upper_bound_adjustment && abs(x_max - 🐺) ≤ upper_bound_adjustment # the singleton AND hamming bounds can fit on the plot - 🍖_idx = searchsortedlast(bins, 🍖) # <- the index of that value - 🍖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🍖_idx] # height of bar of choice - 🐺_idx = searchsortedlast(bins, 🐺) # <- the index of that value - 🐺_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐺_idx] # height of bar of choice - - # adjust emoticon positioning if overlapping values. - if 🐖 == 🍖 - 🍖_y = 🍖_y_top + y_annotation_factor - end - - if 🐖 == 🐺 - if 🐖 == 🍖 # then all of 🍖, 🐖, and 🐺 are the same - 🍖_y = 🍖_y_top + y_annotation_factor - 🐺_y = 🐺_y_top + 2*y_annotation_factor - else - 🐺_y = 🐺_y_top + y_annotation_factor - end - end - - annotate!(plt, 🍖, 🍖_y, "🍖") - annotate!(plt, 🐺, 🐺_y, "🐺") - - 🍖_annotated = true - 🐺_annotated = true - end - - if abs(x_max - 🍖) ≤ upper_bound_adjustment && ! 🍖_annotated # only the hamming bound fits on the plot - 🍖_idx = searchsortedlast(bins, 🍖) # <- the index of that value - 🍖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🍖_idx] # height of bar of choice - - # adjust emoticon positioning if overlapping values. - if 🐖 == 🍖 - 🍖_y = 🍖_y_top + y_annotation_factor - end - - annotate!(plt, 🍖, 🍖_y, "🍖") - - 🍖_annotated = true - end - - if abs(x_max - 🐺) ≤ upper_bound_adjustment && ! 🐺_annotated # only the singleton bound fits on the plot - 🐺_idx = searchsortedlast(bins, 🐺) # <- the index of that value - 🐺_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐺_idx] # height of bar of choice - - # adjust emoticon positioning if overlapping values. - if 🐖 == 🐺 - 🐺_y = 🐺_y_top + y_annotation_factor - end - - annotate!(plt, 🐺, 🐺_y, "🐺") - - 🐺_annotated = true - end - - if ! 🍖_annotated - annotate!(plt, maximum(filter(!isnan, plt.series_list[1].plotattributes[:x])) - x_annotation_factor, maximum(filter(!isnan, plt.series_list[1].plotattributes[:y])) + 0, "🍖 = $(🍖)") - # add_plot_info(plt, extra="🍖 = $(🍖)", counter=annotation_counter); annotation_counter += 1 - end - - if ! 🐺_annotated - annotate!(plt, maximum(filter(!isnan, plt.series_list[1].plotattributes[:x])) - x_annotation_factor, maximum(filter(!isnan, plt.series_list[1].plotattributes[:y])) - y_annotation_factor, "🐺 = $(🐺)") - # add_plot_info(plt, extra="🐺 = $(🐺)", counter=annotation_counter); annotation_counter += 1 - end + data; + fontfamily = font("Times"), + xlabel = "|C|", + ylabel = "Frequency of Times Found", + label = "", + title = "(q, n, d) = ($q, $n, $d): $(format(stop_at, commas = true)) Random Samples, with Selection Size $(format(m, commas = true))", + # xticks = minimum(data) : max(maximum(data), x_max), + bins = bins, + ) + + ## Annotate plot + annotation_counter = 0 + if isempty(data) + add_plot_info(plt; extra = "🐖 = $(🐖)", counter = annotation_counter) + annotation_counter += 1 + add_plot_info(plt; extra = "🍖 = $(🍖)", counter = annotation_counter) + annotation_counter += 1 + add_plot_info(plt; extra = "🐺 = $(🐺)", counter = annotation_counter) + annotation_counter += 1 + savefig(plt, save_path) + println("Plot saved at $save_path.") + return nothing + end + + 🐖_idx = searchsortedlast(bins, 🐖) # <- the index of that value + 🐖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐖_idx] # height of bar of choice + annotate!(plt, 🐖, 🐖_y_top, "🐖") + + annotation_factor = 10 + x_annotation_factor = abs(bins[end] - bins[1]) * 0.05 + y_annotation_factor = abs(bins[end] - bins[1]) * 1.2 + 🍖_y, 🐺_y = 0, 0 + 🍖_annotated, 🐺_annotated = false, false + + if abs(x_max - 🍖) ≤ upper_bound_adjustment && abs(x_max - 🐺) ≤ upper_bound_adjustment # the singleton AND hamming bounds can fit on the plot + 🍖_idx = searchsortedlast(bins, 🍖) # <- the index of that value + 🍖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🍖_idx] # height of bar of choice + 🐺_idx = searchsortedlast(bins, 🐺) # <- the index of that value + 🐺_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐺_idx] # height of bar of choice + + # adjust emoticon positioning if overlapping values. + if 🐖 == 🍖 + 🍖_y = 🍖_y_top + y_annotation_factor + end + + if 🐖 == 🐺 + if 🐖 == 🍖 # then all of 🍖, 🐖, and 🐺 are the same + 🍖_y = 🍖_y_top + y_annotation_factor + 🐺_y = 🐺_y_top + 2 * y_annotation_factor + else + 🐺_y = 🐺_y_top + y_annotation_factor + end + end + + annotate!(plt, 🍖, 🍖_y, "🍖") + annotate!(plt, 🐺, 🐺_y, "🐺") + + 🍖_annotated = true + 🐺_annotated = true + end + + if abs(x_max - 🍖) ≤ upper_bound_adjustment && !🍖_annotated # only the hamming bound fits on the plot + 🍖_idx = searchsortedlast(bins, 🍖) # <- the index of that value + 🍖_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🍖_idx] # height of bar of choice + + # adjust emoticon positioning if overlapping values. + if 🐖 == 🍖 + 🍖_y = 🍖_y_top + y_annotation_factor + end + + annotate!(plt, 🍖, 🍖_y, "🍖") + + 🍖_annotated = true + end + + if abs(x_max - 🐺) ≤ upper_bound_adjustment && !🐺_annotated # only the singleton bound fits on the plot + 🐺_idx = searchsortedlast(bins, 🐺) # <- the index of that value + 🐺_y_top = plt.series_list[1].plotattributes[:y][1:6:end][🐺_idx] # height of bar of choice + + # adjust emoticon positioning if overlapping values. + if 🐖 == 🐺 + 🐺_y = 🐺_y_top + y_annotation_factor + end + + annotate!(plt, 🐺, 🐺_y, "🐺") + + 🐺_annotated = true + end + + if !🍖_annotated + annotate!( + plt, + maximum(filter(!isnan, plt.series_list[1].plotattributes[:x])) - + x_annotation_factor, + maximum(filter(!isnan, plt.series_list[1].plotattributes[:y])) + 0, + "🍖 = $(🍖)", + ) + # add_plot_info(plt, extra="🍖 = $(🍖)", counter=annotation_counter); annotation_counter += 1 + end + + if !🐺_annotated + annotate!( + plt, + maximum(filter(!isnan, plt.series_list[1].plotattributes[:x])) - + x_annotation_factor, + maximum(filter(!isnan, plt.series_list[1].plotattributes[:y])) - + y_annotation_factor, + "🐺 = $(🐺)", + ) + # add_plot_info(plt, extra="🐺 = $(🐺)", counter=annotation_counter); annotation_counter += 1 + end savefig(plt, save_path) println("Plot saved at $save_path.") - return nothing + return nothing end -graphing(q::Integer, n::Integer, d::Integer, stop_at::Integer; m::Integer=10_000) = - graphing(q, n, d, obtain_data(q, n, d, stop_at, m = m), stop_at, m = m) +function graphing(q::Integer, n::Integer, d::Integer, stop_at::Integer; m::Integer = 10_000) + return graphing(q, n, d, obtain_data(q, n, d, stop_at; m = m), stop_at; m = m) +end diff --git a/examples/plot_random_from_file.jl b/examples/plot_random_from_file.jl index 75a84c3..a2e09a9 100755 --- a/examples/plot_random_from_file.jl +++ b/examples/plot_random_from_file.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia -tauto --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - #= Example usage: ``` @@ -15,13 +9,13 @@ include("plot_random_core.jl") # reading from pre-obtained data if isempty(ARGS) || length(ARGS) < 5 - throw(error(""" - Please specify the dataset from the command line, along with q, n, and d. For example - ``` - ./examples/plot_random_from_file.jl /Users/jakeireland/Desktop/random_number_in_code,q=21,n=4,d=3,i=1000,m=10000.csv 6 7 3 10000# dataset, q, n, d, selection_set_size - ``` - That is respectively: the size of the alphabet; the block length of the words; the distance of the code; the number of datapoints in the code; and the selection set size of each "random" sample. - """)) + throw(error(""" + Please specify the dataset from the command line, along with q, n, and d. For example + ``` + ./examples/plot_random_from_file.jl /Users/jakeireland/Desktop/random_number_in_code,q=21,n=4,d=3,i=1000,m=10000.csv 6 7 3 10000# dataset, q, n, d, selection_set_size + ``` + That is respectively: the size of the alphabet; the block length of the words; the distance of the code; the number of datapoints in the code; and the selection set size of each "random" sample. + """,),) end global data = ARGS[1] # e.g. "/Users/jakeireland/Desktop/random_number_in_code,q=21,n=4,d=3,i=1000,m=10000.csv" @@ -37,4 +31,4 @@ global 🐺 = singleton_bound(q, n, d) # data = CSV.read("/Users/jakeireland/projects/CodingTheory.jl/other/random_number_analysis/random_number_in_code,q=$(q),n=$(n),d=$(d),i=$(stop_at),m=$(m).csv") # data = CSV.read("/Users/jakeireland/projects/CodingTheory.jl/other/random_number_analysis/mildly_interesting_ones/[2]-10-3-2/random_number_in_code,q=10,n=3,d=2,i=1000,m=10000.csv") data = CSV.read(data) -graphing(q, n, d, data, nrow(data), m = m) +graphing(q, n, d, data, nrow(data); m = m) diff --git a/examples/random_numbers_core.jl b/examples/random_numbers_core.jl index f45d082..fe5ef56 100755 --- a/examples/random_numbers_core.jl +++ b/examples/random_numbers_core.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia -tauto --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - println("Loading required libraries...\n") using Base.Threads: @threads @@ -16,22 +10,29 @@ plotly() include(joinpath(dirname(dirname(@__FILE__)), "src", "CodingTheory.jl")) using .CodingTheory -function obtain_data(q::Integer, n::Integer, d::Integer, stop_at::Integer; m::Integer=10_000) - n < d && throw(error("You cannot have a distance greater than the block length.")) +function obtain_data( + q::Integer, n::Integer, d::Integer, stop_at::Integer; m::Integer = 10_000 +) + n < d && throw(error("You cannot have a distance greater than the block length.")) data = [] - data_path = joinpath(dirname(@__DIR__), "other", "random_number_analysis", "random_number_in_code,q=$(q),n=$(n),d=$(d),i=$(stop_at),m=$(m).csv") - - ## Obtain data - p = Progress(stop_at, 1) # minimum update interval: 1 second + data_path = joinpath( + dirname(@__DIR__), + "other", + "random_number_analysis", + "random_number_in_code,q=$(q),n=$(n),d=$(d),i=$(stop_at),m=$(m).csv", + ) + + ## Obtain data + p = Progress(stop_at, 1) # minimum update interval: 1 second @floop ThreadedEx() for _ in 1:stop_at - push!(data, length(get_codewords_random(q, n, d, m = m))) - next!(p) # increment progress bar + push!(data, length(get_codewords_random(q, n, d; m = m))) + next!(p) # increment progress bar end - - df = DataFrame(size_of_code = data) - CSV.write(data_path, df) - - println("Wrote data to $(data_path)") - - return data + + df = DataFrame(; size_of_code = data) + CSV.write(data_path, df) + + println("Wrote data to $(data_path)") + + return data end diff --git a/examples/sort_bounds.jl b/examples/sort_bounds.jl index c1da4cd..85b896d 100755 --- a/examples/sort_bounds.jl +++ b/examples/sort_bounds.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - include(joinpath(dirname(@__DIR__), "src", "CodingTheory.jl")) using .CodingTheory @@ -13,17 +7,22 @@ using Formatting # using Lazy: @> const ruled_out = Tuple[ - # (3, 11, 5), # n words ≈ 186 ≠ 729 - # (21, 4, 3), # n words ≈ 359 ≠ 2401 - # (6, 7, 3), # n words ≈ 2995 ≠ 7776 - # (2, 23, 7), # this is one of the Golay codes, silly me - ] +# (3, 11, 5), # n words ≈ 186 ≠ 729 +# (21, 4, 3), # n words ≈ 359 ≠ 2401 +# (6, 7, 3), # n words ≈ 2995 ≠ 7776 +# (2, 23, 7), # this is one of the Golay codes, silly me +] function load_df(pathname::AbstractString) # use bounds.jl to get the hamming bounds df = DataFrame(CSV.read(pathname)) println("The complete list has $(format(nrow(df), commas = true)) observations.") - insert!(df, ncol(df) + 1, CodingTheory.sizeof_perfect_code.(df.q, df.n, df.d), :size_of_perfect_code) + insert!( + df, + ncol(df) + 1, + CodingTheory.sizeof_perfect_code.(df.q, df.n, df.d), + :size_of_perfect_code, + ) df = round.(BigInt, df) return df @@ -31,53 +30,52 @@ end function filter_df( df::DataFrame; - q_not_one::Bool=true, - d_not_one::Bool=true, - d_not_two::Bool=true, - d_not_n::Bool=true, - d_leq_n::Bool=true, - not_hamming_perfect::Bool=true, - d_not_even::Bool=true, - q_prime_power::Bool=false, - q_not_prime_power::Bool=true, - hamming_bound_isinteger::Bool=true, - hamming_bound_is_lower_bound::Bool=true - ) - + q_not_one::Bool = true, + d_not_one::Bool = true, + d_not_two::Bool = true, + d_not_n::Bool = true, + d_leq_n::Bool = true, + not_hamming_perfect::Bool = true, + d_not_even::Bool = true, + q_prime_power::Bool = false, + q_not_prime_power::Bool = true, + hamming_bound_isinteger::Bool = true, + hamming_bound_is_lower_bound::Bool = true, +) if q_not_one df = @where(df, :q .≠ 1) end - + if d_not_one df = @where(df, :d .≠ 1) end - + if d_not_two df = @where(df, :d .≠ 2) end - + if d_not_n df = @where(df, :d .≠ :n) end - + if d_leq_n df = @where(df, :d .≤ :n) end - + if not_hamming_perfect - df = @where(df, .! ishammingperfect.(:q, :n, :d)) + df = @where(df, .!ishammingperfect.(:q, :n, :d)) end - + if d_not_even - df = @where(df, .! iseven.(:d)) + df = @where(df, .!iseven.(:d)) end - + if q_prime_power df = @where(df, CodingTheory.isprimepower.(:q)) end - + if q_not_prime_power - df = @where(df, .! CodingTheory.isprimepower.(:q)) + df = @where(df, .!CodingTheory.isprimepower.(:q)) end if hamming_bound_isinteger @@ -85,60 +83,56 @@ function filter_df( end if hamming_bound_is_lower_bound - df = @where(df, isequal.(:hamming_bound, :smallest_bound)) + df = @where(df, isequal.(:hamming_bound, :smallest_bound)) # df = @where(df, isequal.(:hamming_bound, min.(:hamming_bound, :singleton_bound))) end - + # Filter out those (q, n, d) configurations already ruled out. # df = @where(df, tuple.(:q, :n, :d) .∉ Ref(ruled_out)) - + return df end function view_df( df::DataFrame; - nrows::Integer=10, - q_not_one::Bool=true, - d_not_one::Bool=true, - d_not_two::Bool=true, - d_not_n::Bool=true, - d_leq_n::Bool=true, - not_hamming_perfect::Bool=true, - d_not_even::Bool=true, - q_prime_power::Bool=false, - q_not_prime_power::Bool=true, - hamming_bound_isinteger::Bool=true, - hamming_bound_is_lower_bound::Bool=true - ) + nrows::Integer = 10, + q_not_one::Bool = true, + d_not_one::Bool = true, + d_not_two::Bool = true, + d_not_n::Bool = true, + d_leq_n::Bool = true, + not_hamming_perfect::Bool = true, + d_not_even::Bool = true, + q_prime_power::Bool = false, + q_not_prime_power::Bool = true, + hamming_bound_isinteger::Bool = true, + hamming_bound_is_lower_bound::Bool = true, +) df = filter_df( df; - q_not_one=q_not_one, - d_not_one=d_not_one, - d_not_two=d_not_two, - d_not_n=d_not_n, - d_leq_n=d_leq_n, - not_hamming_perfect=not_hamming_perfect, - d_not_even=d_not_even, - q_prime_power=q_prime_power, - q_not_prime_power=q_not_prime_power, - hamming_bound_isinteger=hamming_bound_isinteger, - hamming_bound_is_lower_bound=hamming_bound_is_lower_bound - ) - + q_not_one = q_not_one, + d_not_one = d_not_one, + d_not_two = d_not_two, + d_not_n = d_not_n, + d_leq_n = d_leq_n, + not_hamming_perfect = not_hamming_perfect, + d_not_even = d_not_even, + q_prime_power = q_prime_power, + q_not_prime_power = q_not_prime_power, + hamming_bound_isinteger = hamming_bound_isinteger, + hamming_bound_is_lower_bound = hamming_bound_is_lower_bound, + ) + # sort last! df = sort(df, :hamming_bound) # df = sort(df, :size_of_perfect_code) - - println("The filtered list has $(format(nrow(df), commas = true)) observations matching the filtered criteria. Now showing the 10 smallest values.") - # return select(df, Not([:smallest_bound]))[1:nrows, :] - return df[1:nrows, :] + + println("The filtered list has $(format(nrow(df), commas = true)) observations matching the filtered criteria. Now showing the 10 smallest values.",) + # return select(df, Not([:smallest_bound]))[1:nrows, :] + return df[1:nrows, :] # return select(df, Not([:singleton_bound]))[1:nrows, :] end - - - - ############################################################################################################## datafile = "/Users/jakeireland/projects/CodingTheory.jl/other/research/bounds/bound_integers_100000_very_loose_q_neq_1.csv" # <===== @@ -153,19 +147,19 @@ datafile = "/Users/jakeireland/projects/CodingTheory.jl/other/research/bounds/bo df = view_df( load_df(datafile); - nrows=10, - q_not_one=true, - d_not_one=true, # distance equal to one is trivially the list of all words with q and n - d_not_two=true, # distance equal to two is trivial as the words' coordinates sum to zero modulo q - d_not_n=true, # distance equal to two is trivially the repetition code - d_leq_n=true, - not_hamming_perfect=true, - d_not_even=false, # cannot get a perfect code if d is even, because of how hamming bound is calculated - q_prime_power=false, # inverse of the below - q_not_prime_power=true, # we are interested in alphabet size not prime power - hamming_bound_isinteger=true, - hamming_bound_is_lower_bound=true # cannnot obtain perfect code if singleton bound is less than hamming bound - ) + nrows = 10, + q_not_one = true, + d_not_one = true, # distance equal to one is trivially the list of all words with q and n + d_not_two = true, # distance equal to two is trivial as the words' coordinates sum to zero modulo q + d_not_n = true, # distance equal to two is trivially the repetition code + d_leq_n = true, + not_hamming_perfect = true, + d_not_even = false, # cannot get a perfect code if d is even, because of how hamming bound is calculated + q_prime_power = false, # inverse of the below + q_not_prime_power = true, # we are interested in alphabet size not prime power + hamming_bound_isinteger = true, + hamming_bound_is_lower_bound = true, # cannnot obtain perfect code if singleton bound is less than hamming bound +) println(df) @@ -198,10 +192,6 @@ println(df) # println(size(D)) # println(sort(D, :size_of_perfect_code)) - - - - ## original code: # df = round.(BigInt, DataFrame(CSV.read(datafile))) # df = @where(df, .! isprimepower.(:q)) diff --git a/examples/turtle.jl b/examples/turtle.jl index 48811b0..7e8d96a 100644 --- a/examples/turtle.jl +++ b/examples/turtle.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - using CSV, DataFrames, Luxor, Colors, Cairo include(joinpath(dirname(dirname(@__FILE__)), "src", "messages.jl")) @@ -15,17 +9,17 @@ function integer_search(stop_at::Integer)::Array{Array{Number, 1}} A = [] upper_bound = 10 increment_bound = 10 - processed = [] - i = 0 - + processed = [] + i = 0 + while true for q in 1:upper_bound, n in 1:upper_bound, d in 1:upper_bound - # skip configurations that have already been processed - # arelessthan(upper_bound - increment_bound, q, n, d) && continue - (q, n, d) ∈ processed && continue + # skip configurations that have already been processed + # arelessthan(upper_bound - increment_bound, q, n, d) && continue + (q, n, d) ∈ processed && continue hb = hamming_bound(q, n, d, no_round) - push!(processed, (q, n, d)) - # i += 1; println("$i:\t$q, $n, $d") + push!(processed, (q, n, d)) + # i += 1; println("$i:\t$q, $n, $d") push!(A, [q, n, d, hb]) isequal(length(A), stop_at) && return A end @@ -35,60 +29,76 @@ end function make_integer_df(stop_at::Integer) A = integer_search(stop_at) - D = DataFrame(q = Number[], n = Number[], d = Number[], hamming_bound = Number[]) + D = DataFrame(; q = Number[], n = Number[], d = Number[], hamming_bound = Number[]) for i in A push!(D, i) end - - return D + + return D end -function draw_turtle(stop_at::Integer; specific::Bool=false) - rigidness = "non_rigid" - if specific - rigidness = "rigid" - end - - df = make_integer_df(stop_at) - - golden_angle = 137.5 - angle = golden_angle - # angle = 42 - - # Drawing(600, 400, joinpath(dirname(dirname(@__FILE__)), "other", "turtle_drawing.pdf")) - Drawing(10000, 10000, joinpath(homedir(), "projects", "CodingTheory.jl", "other", "$(rigidness)", "turtle_drawing_$(stop_at)_$(angle).pdf")) - origin() - # background("midnightblue") +function draw_turtle(stop_at::Integer; specific::Bool = false) + rigidness = "non_rigid" + if specific + rigidness = "rigid" + end + + df = make_integer_df(stop_at) + + golden_angle = 137.5 + angle = golden_angle + # angle = 42 + + # Drawing(600, 400, joinpath(dirname(dirname(@__FILE__)), "other", "turtle_drawing.pdf")) + Drawing( + 10000, + 10000, + joinpath( + homedir(), + "projects", + "CodingTheory.jl", + "other", + "$(rigidness)", + "turtle_drawing_$(stop_at)_$(angle).pdf", + ), + ) + origin() + # background("midnightblue") + + 🐢 = Turtle() + Pencolor(🐢, "black") + Penwidth(🐢, 1.5) + n = 5 + + # distance shouldn't be larger than the block length; filter trivial distances of one; we filter out combinations when all are equal; filter out codes that are generalised hamming codes; filter out even distances + for row in eachrow(df) + if specific + if row.d < row.n && + !isone(row.d) && + !allequal(row.q, row.n, row.d) && + !iszero(mod(row.d, 2)) && + !isone(row.hamming_bound) && + !ishammingbound(BigInt(row.q), BigInt(row.n), BigInt(row.d)) + continue + end + else + if row.d < row.n && !isone(row.d) + continue + end + end + + if isinteger(row.hamming_bound) + Turn(🐢, angle) + # HueShift(🐢) + end + Forward(🐢, n) + # n += 0.75 + end + + # fontsize(20) + # Message(🐢, "finished") + Luxor.finish() - 🐢 = Turtle() - Pencolor(🐢, "black") - Penwidth(🐢, 1.5) - n = 5 - - # distance shouldn't be larger than the block length; filter trivial distances of one; we filter out combinations when all are equal; filter out codes that are generalised hamming codes; filter out even distances - for row in eachrow(df) - if specific - if row.d < row.n && ! isone(row.d) && ! allequal(row.q, row.n, row.d) && ! iszero(mod(row.d, 2)) && ! isone(row.hamming_bound) && ! ishammingbound(BigInt(row.q), BigInt(row.n), BigInt(row.d)) - continue - end - else - if row.d < row.n && ! isone(row.d) - continue - end - end - - if isinteger(row.hamming_bound) - Turn(🐢, angle) - # HueShift(🐢) - end - Forward(🐢, n) - # n += 0.75 - end - - # fontsize(20) - # Message(🐢, "finished") - Luxor.finish() - - return nothing + return nothing end diff --git a/justfile b/justfile index 58f2295..b24d2ec 100644 --- a/justfile +++ b/justfile @@ -7,14 +7,14 @@ # as the root project directory that Julia uses. # # [1]: https://just.systems/man/en/functions.html#justfile-and-justfile-directory -project_dir := justfile_dir() + "/" -test_dir := project_dir + "test/" -test_file := test_dir + "runtests.jl" -docs_dir := project_dir + "docs/" -docs_mk_file := docs_dir + "make.jl" -dev_dir := project_dir + "dev/" -bench_dir := project_dir + "perf/" -bench_file := bench_dir + "runbenchmarks.jl" +project_dir := justfile_dir() / "/" +test_dir := project_dir / "test/" +test_file := test_dir / "runtests.jl" +docs_dir := project_dir / "docs/" +docs_mk_file := docs_dir / "make.jl" +dev_dir := project_dir / "dev/" +bench_dir := project_dir / "perf/" +bench_file := bench_dir / "runbenchmarks.jl" standard_instantiate_code := """ import Pkg Pkg.instantiate() @@ -43,10 +43,11 @@ docs: (instantiate-dev docs_dir) bench: (instantiate-dev bench_dir) julia --project={{bench_dir}} {{bench_file}} -# check formatting +# Check formatting with blue style [group: 'ci'] fmt: - julia --project={{dev_dir}} -e 'using JuliaFormatter; format(".")' + # https://github.com/invenia/BlueStyle + julia --project={{dev_dir}} -e 'using JuliaFormatter; format("{{project_dir}}", style=BlueStyle())' # Instantiate main project instantiate: diff --git a/perf/runbenchmarks.jl b/perf/runbenchmarks.jl index 75b1dad..3cb861b 100644 --- a/perf/runbenchmarks.jl +++ b/perf/runbenchmarks.jl @@ -5,8 +5,19 @@ function BM() hamming_distance("ABC", "DBC") == 1 hamming_distance("ABC", "DEF") == 3 - hamming_ball([[1, 0, 1], [0, 1, 1], [1, 0, 0]], [1, 0, 0], 2) == deepsym([[1, 0, 1], [1, 0, 0]]) - hamming_ball(list_span([2, 2, 2], [1, 2, 0], [0, 2, 1], 3), [1, 0, 0], 3) == deepsym([[0, 0, 0], [0, 2, 1], [0, 1, 2], [1, 2, 0], [1, 1, 1], [1, 0, 2], [2, 1, 0], [2, 0, 1], [2, 2, 2]]) + hamming_ball([[1, 0, 1], [0, 1, 1], [1, 0, 0]], [1, 0, 0], 2) == + deepsym([[1, 0, 1], [1, 0, 0]]) + hamming_ball(list_span([2, 2, 2], [1, 2, 0], [0, 2, 1], 3), [1, 0, 0], 3) == deepsym([ + [0, 0, 0], + [0, 2, 1], + [0, 1, 2], + [1, 2, 0], + [1, 1, 1], + [1, 0, 2], + [2, 1, 0], + [2, 0, 1], + [2, 2, 2], + ]) rate(3, 5, 4) ≈ 0.3662433802 @@ -14,8 +25,11 @@ function BM() t_error_correcting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false t_error_detecting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false - find_error_detection_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) == 1 - find_error_correction_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) == 0 + find_error_detection_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) == + 1 + find_error_correction_max( + [[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2 + ) == 0 find_error_correction_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 0 find_error_detection_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 1 @@ -37,35 +51,150 @@ function BM() mod(rem(a, b), 2) == Polynomial([1]) rref([1 0 1 1 1; 1 1 1 0 1; 0 1 1 1 1], 2) == [1 0 0 1 0; 0 1 0 1 0; 0 0 1 0 1] - rref([1 0 1 0 1 0; 0 1 0 0 1 0; 1 1 1 1 1 1], 2) == [1 0 1 0 1 0; 0 1 0 0 1 0; 0 0 0 1 1 1] - rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=false) == [1 0 3 0 2 2; 0 1 2 0 1 1; 0 0 0 1 0 4] - rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=true) == [1 0 0 3 2 2; 0 1 0 2 1 1; 0 0 1 0 0 4] - rref([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] - rref([0 0 0 0 0; 1 0 1 0 1; 0 1 0 1 1; 1 1 1 1 0], 2) == [1 0 1 0 1; 0 1 0 1 1; 0 0 0 0 0; 0 0 0 0 0] + rref([1 0 1 0 1 0; 0 1 0 0 1 0; 1 1 1 1 1 1], 2) == + [1 0 1 0 1 0; 0 1 0 0 1 0; 0 0 0 1 1 1] + rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5; colswap = false) == + [1 0 3 0 2 2; 0 1 2 0 1 1; 0 0 0 1 0 4] + rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5; colswap = true) == + [1 0 0 3 2 2; 0 1 0 2 1 1; 0 0 1 0 0 4] + rref([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] + rref([0 0 0 0 0; 1 0 1 0 1; 0 1 0 1 1; 1 1 1 1 0], 2) == + [1 0 1 0 1; 0 1 0 1 1; 0 0 0 0 0; 0 0 0 0 0] rref([1 1 1 0; 1 1 0 1; 0 0 1 1], 2) == [1 1 0 1; 0 0 1 1; 0 0 0 0] - multiplication_table(2, 3) == Polynomial[Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]); Polynomial([0]) Polynomial([1]) Polynomial([2]) Polynomial([0, 1]) Polynomial([1, 1]) Polynomial([2, 1]) Polynomial([0, 2]) Polynomial([1, 2]) Polynomial([2, 2]); Polynomial([0]) Polynomial([2]) Polynomial([1]) Polynomial([0, 2]) Polynomial([2, 2]) Polynomial([1, 2]) Polynomial([0, 1]) Polynomial([2, 1]) Polynomial([1, 1]); Polynomial([0]) Polynomial([0, 1]) Polynomial([0, 2]) Polynomial([0, 0, 1]) Polynomial([0, 1, 1]) Polynomial([0, 2, 1]) Polynomial([0, 0, 2]) Polynomial([0, 1, 2]) Polynomial([0, 2, 2]); Polynomial([0]) Polynomial([1, 1]) Polynomial([2, 2]) Polynomial([0, 1, 1]) Polynomial([1, 2, 1]) Polynomial([2, 0, 1]) Polynomial([0, 2, 2]) Polynomial([1, 0, 2]) Polynomial([2, 1, 2]); Polynomial([0]) Polynomial([2, 1]) Polynomial([1, 2]) Polynomial([0, 2, 1]) Polynomial([2, 0, 1]) Polynomial([1, 1, 1]) Polynomial([0, 1, 2]) Polynomial([2, 2, 2]) Polynomial([1, 0, 2]); Polynomial([0]) Polynomial([0, 2]) Polynomial([0, 1]) Polynomial([0, 0, 2]) Polynomial([0, 2, 2]) Polynomial([0, 1, 2]) Polynomial([0, 0, 1]) Polynomial([0, 2, 1]) Polynomial([0, 1, 1]); Polynomial([0]) Polynomial([1, 2]) Polynomial([2, 1]) Polynomial([0, 1, 2]) Polynomial([1, 0, 2]) Polynomial([2, 2, 2]) Polynomial([0, 2, 1]) Polynomial([1, 1, 1]) Polynomial([2, 0, 1]); Polynomial([0]) Polynomial([2, 2]) Polynomial([1, 1]) Polynomial([0, 2, 2]) Polynomial([2, 1, 2]) Polynomial([1, 0, 2]) Polynomial([0, 1, 1]) Polynomial([2, 0, 1]) Polynomial([1, 2, 1])] + multiplication_table(2, 3) == Polynomial[ + Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) + Polynomial([0]) Polynomial([1]) Polynomial([2]) Polynomial([0, 1]) Polynomial([1, 1]) Polynomial([2, 1]) Polynomial([0, 2]) Polynomial([1, 2]) Polynomial([2, 2]) + Polynomial([0]) Polynomial([2]) Polynomial([1]) Polynomial([0, 2]) Polynomial([2, 2]) Polynomial([1, 2]) Polynomial([0, 1]) Polynomial([2, 1]) Polynomial([1, 1]) + Polynomial([0]) Polynomial([0, 1]) Polynomial([0, 2]) Polynomial([0, 0, 1]) Polynomial([0, 1, 1]) Polynomial([0, 2, 1]) Polynomial([0, 0, 2]) Polynomial([0, 1, 2]) Polynomial([0, 2, 2]) + Polynomial([0]) Polynomial([1, 1]) Polynomial([2, 2]) Polynomial([0, 1, 1]) Polynomial([1, 2, 1]) Polynomial([2, 0, 1]) Polynomial([0, 2, 2]) Polynomial([1, 0, 2]) Polynomial([2, 1, 2]) + Polynomial([0]) Polynomial([2, 1]) Polynomial([1, 2]) Polynomial([0, 2, 1]) Polynomial([2, 0, 1]) Polynomial([1, 1, 1]) Polynomial([0, 1, 2]) Polynomial([2, 2, 2]) Polynomial([1, 0, 2]) + Polynomial([0]) Polynomial([0, 2]) Polynomial([0, 1]) Polynomial([0, 0, 2]) Polynomial([0, 2, 2]) Polynomial([0, 1, 2]) Polynomial([0, 0, 1]) Polynomial([0, 2, 1]) Polynomial([0, 1, 1]) + Polynomial([0]) Polynomial([1, 2]) Polynomial([2, 1]) Polynomial([0, 1, 2]) Polynomial([1, 0, 2]) Polynomial([2, 2, 2]) Polynomial([0, 2, 1]) Polynomial([1, 1, 1]) Polynomial([2, 0, 1]) + Polynomial([0]) Polynomial([2, 2]) Polynomial([1, 1]) Polynomial([0, 2, 2]) Polynomial([2, 1, 2]) Polynomial([1, 0, 2]) Polynomial([0, 1, 1]) Polynomial([2, 0, 1]) Polynomial([1, 2, 1]) + ] - list_span([2, 1, 1], [1, 1, 1], 3) == [[0, 0, 0], [1, 1, 1], [2, 2, 2], [2, 1, 1], [0, 2, 2], [1, 0, 0], [1, 2, 2], [2, 0, 0], [0, 1, 1]] + list_span([2, 1, 1], [1, 1, 1], 3) == [ + [0, 0, 0], + [1, 1, 1], + [2, 2, 2], + [2, 1, 1], + [0, 2, 2], + [1, 0, 0], + [1, 2, 2], + [2, 0, 0], + [0, 1, 1], + ] - islinear([[0,0,0],[1,1,1],[1,0,1],[1,1,0]], 2) == false - islinear([[0,0,0],[1,1,1],[1,0,1],[0,1,0]], 2) == true + islinear([[0, 0, 0], [1, 1, 1], [1, 0, 1], [1, 1, 0]], 2) == false + islinear([[0, 0, 0], [1, 1, 1], [1, 0, 1], [0, 1, 0]], 2) == true - code_distance([[0,0,0,0,0],[1,0,1,0,1],[0,1,0,1,0],[1,1,1,1,1]]) == 2 - code_distance([[0,0,0,0,0],[1,1,1,0,0],[0,0,0,1,1],[1,1,1,1,1],[1,0,0,1,1],[0,1,1,0,0]]) == 1 + code_distance([[0, 0, 0, 0, 0], [1, 0, 1, 0, 1], [0, 1, 0, 1, 0], [1, 1, 1, 1, 1]]) == 2 + code_distance([ + [0, 0, 0, 0, 0], + [1, 1, 1, 0, 0], + [0, 0, 0, 1, 1], + [1, 1, 1, 1, 1], + [1, 0, 0, 1, 1], + [0, 1, 1, 0, 0], + ]) == 1 Alphabet("123") == deepsym([1, 2, 3]) Alphabet([1, 2, 3]) == deepsym([1, 2, 3]) Alphabet(["1", "2", "3"]) == deepsym([1, 2, 3]) # TODO: write test for CodeUniverse struct [i for i in CodeUniverseIterator(["a", "b", "c"], 3)] == get_all_words(["a", "b", "c"], 3) # implicitly tests CodeUniverseIterator - collect(CodeUniverseIterator(["a", "b", "c"], 4)) == Tuple[(:a, :a, :a, :a), (:b, :a, :a, :a), (:c, :a, :a, :a), (:a, :b, :a, :a), (:b, :b, :a, :a), (:c, :b, :a, :a), (:a, :c, :a, :a), (:b, :c, :a, :a), (:c, :c, :a, :a), (:a, :a, :b, :a), (:b, :a, :b, :a), (:c, :a, :b, :a), (:a, :b, :b, :a), (:b, :b, :b, :a), (:c, :b, :b, :a), (:a, :c, :b, :a), (:b, :c, :b, :a), (:c, :c, :b, :a), (:a, :a, :c, :a), (:b, :a, :c, :a), (:c, :a, :c, :a), (:a, :b, :c, :a), (:b, :b, :c, :a), (:c, :b, :c, :a), (:a, :c, :c, :a), (:b, :c, :c, :a), (:c, :c, :c, :a), (:a, :a, :a, :b), (:b, :a, :a, :b), (:c, :a, :a, :b), (:a, :b, :a, :b), (:b, :b, :a, :b), (:c, :b, :a, :b), (:a, :c, :a, :b), (:b, :c, :a, :b), (:c, :c, :a, :b), (:a, :a, :b, :b), (:b, :a, :b, :b), (:c, :a, :b, :b), (:a, :b, :b, :b), (:b, :b, :b, :b), (:c, :b, :b, :b), (:a, :c, :b, :b), (:b, :c, :b, :b), (:c, :c, :b, :b), (:a, :a, :c, :b), (:b, :a, :c, :b), (:c, :a, :c, :b), (:a, :b, :c, :b), (:b, :b, :c, :b), (:c, :b, :c, :b), (:a, :c, :c, :b), (:b, :c, :c, :b), (:c, :c, :c, :b), (:a, :a, :a, :c), (:b, :a, :a, :c), (:c, :a, :a, :c), (:a, :b, :a, :c), (:b, :b, :a, :c), (:c, :b, :a, :c), (:a, :c, :a, :c), (:b, :c, :a, :c), (:c, :c, :a, :c), (:a, :a, :b, :c), (:b, :a, :b, :c), (:c, :a, :b, :c), (:a, :b, :b, :c), (:b, :b, :b, :c), (:c, :b, :b, :c), (:a, :c, :b, :c), (:b, :c, :b, :c), (:c, :c, :b, :c), (:a, :a, :c, :c), (:b, :a, :c, :c), (:c, :a, :c, :c), (:a, :b, :c, :c), (:b, :b, :c, :c), (:c, :b, :c, :c), (:a, :c, :c, :c), (:b, :c, :c, :c), (:c, :c, :c, :c)] + collect(CodeUniverseIterator(["a", "b", "c"], 4)) == Tuple[ + (:a, :a, :a, :a), + (:b, :a, :a, :a), + (:c, :a, :a, :a), + (:a, :b, :a, :a), + (:b, :b, :a, :a), + (:c, :b, :a, :a), + (:a, :c, :a, :a), + (:b, :c, :a, :a), + (:c, :c, :a, :a), + (:a, :a, :b, :a), + (:b, :a, :b, :a), + (:c, :a, :b, :a), + (:a, :b, :b, :a), + (:b, :b, :b, :a), + (:c, :b, :b, :a), + (:a, :c, :b, :a), + (:b, :c, :b, :a), + (:c, :c, :b, :a), + (:a, :a, :c, :a), + (:b, :a, :c, :a), + (:c, :a, :c, :a), + (:a, :b, :c, :a), + (:b, :b, :c, :a), + (:c, :b, :c, :a), + (:a, :c, :c, :a), + (:b, :c, :c, :a), + (:c, :c, :c, :a), + (:a, :a, :a, :b), + (:b, :a, :a, :b), + (:c, :a, :a, :b), + (:a, :b, :a, :b), + (:b, :b, :a, :b), + (:c, :b, :a, :b), + (:a, :c, :a, :b), + (:b, :c, :a, :b), + (:c, :c, :a, :b), + (:a, :a, :b, :b), + (:b, :a, :b, :b), + (:c, :a, :b, :b), + (:a, :b, :b, :b), + (:b, :b, :b, :b), + (:c, :b, :b, :b), + (:a, :c, :b, :b), + (:b, :c, :b, :b), + (:c, :c, :b, :b), + (:a, :a, :c, :b), + (:b, :a, :c, :b), + (:c, :a, :c, :b), + (:a, :b, :c, :b), + (:b, :b, :c, :b), + (:c, :b, :c, :b), + (:a, :c, :c, :b), + (:b, :c, :c, :b), + (:c, :c, :c, :b), + (:a, :a, :a, :c), + (:b, :a, :a, :c), + (:c, :a, :a, :c), + (:a, :b, :a, :c), + (:b, :b, :a, :c), + (:c, :b, :a, :c), + (:a, :c, :a, :c), + (:b, :c, :a, :c), + (:c, :c, :a, :c), + (:a, :a, :b, :c), + (:b, :a, :b, :c), + (:c, :a, :b, :c), + (:a, :b, :b, :c), + (:b, :b, :b, :c), + (:c, :b, :b, :c), + (:a, :c, :b, :c), + (:b, :c, :b, :c), + (:c, :c, :b, :c), + (:a, :a, :c, :c), + (:b, :a, :c, :c), + (:c, :a, :c, :c), + (:a, :b, :c, :c), + (:b, :b, :c, :c), + (:c, :b, :c, :c), + (:a, :c, :c, :c), + (:b, :c, :c, :c), + (:c, :c, :c, :c), + ] - sphere_covering_bound(5,7,3) == 215 - sphere_packing_bound(5,7,3) == 2693 + sphere_covering_bound(5, 7, 3) == 215 + sphere_packing_bound(5, 7, 3) == 2693 - construct_ham_matrix(3,2) == [0 0 0 1 1 1 1; 0 1 1 0 0 1 1; 1 0 1 0 1 0 1] - construct_ham_matrix(3,3) == [0 0 0 0 0 0 0 0 1 1 1 1 1; 0 0 1 1 1 2 2 2 0 0 0 1 1; 1 2 0 1 2 0 1 2 0 1 2 0 1] + construct_ham_matrix(3, 2) == [0 0 0 1 1 1 1; 0 1 1 0 0 1 1; 1 0 1 0 1 0 1] + construct_ham_matrix(3, 3) == + [0 0 0 0 0 0 0 0 1 1 1 1 1; 0 0 1 1 1 2 2 2 0 0 0 1 1; 1 2 0 1 2 0 1 2 0 1 2 0 1] isperfect(11, 6, 5, 3) == true isperfect(23, 12, 7, 2) == true @@ -81,16 +210,49 @@ function BM() length(get_codewords_greedy(5, 5, 3)) == 74 randq, randn = rand(1:8, 2) length(get_all_words(randq, randn)) == big(randq)^randn - get_codewords([1 0 1 0; 0 1 1 1], 2) == [[0, 0, 0, 0], [1, 0, 1, 0], [0, 1, 1, 1], [1, 1, 0, 1]] - get_codewords([1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1], 2) == [[0, 0, 0, 0, 0, 0], [1, 0, 0, 1, 1, 0], [0, 1, 0, 1, 0, 1], [1, 1, 0, 0, 1, 1], [0, 0, 1, 0, 1, 1], [1, 0, 1, 1, 0, 1], [0, 1, 1, 1, 1, 0], [1, 1, 1, 0, 0, 0]] - - syndrome([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == [0 0 0] - parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] - normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] - equivalent_code([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] - parity_check(normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3), 3) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] - isincode([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == true - isincode([1, 0, 2, 2, 1, 2, 1], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == false + get_codewords([1 0 1 0; 0 1 1 1], 2) == + [[0, 0, 0, 0], [1, 0, 1, 0], [0, 1, 1, 1], [1, 1, 0, 1]] + get_codewords([1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1], 2) == [ + [0, 0, 0, 0, 0, 0], + [1, 0, 0, 1, 1, 0], + [0, 1, 0, 1, 0, 1], + [1, 1, 0, 0, 1, 1], + [0, 0, 1, 0, 1, 1], + [1, 0, 1, 1, 0, 1], + [0, 1, 1, 1, 1, 0], + [1, 1, 1, 0, 0, 0], + ] + + syndrome( + [0, 2, 1, 2, 0, 1, 0], + transpose(parity_check( + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3 + )), + 3, + ) == [0 0 0] + parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3) == + [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] + normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] + equivalent_code([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] + parity_check( + normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3), 3 + ) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] + isincode( + [0, 2, 1, 2, 0, 1, 0], + transpose(parity_check( + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3 + )), + 3, + ) == true + return isincode( + [1, 0, 2, 2, 1, 2, 1], + transpose(parity_check( + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3 + )), + 3, + ) == false end # end runtests @btime BM() diff --git a/src/CodingTheory.jl b/src/CodingTheory.jl index ea90beb..293c957 100755 --- a/src/CodingTheory.jl +++ b/src/CodingTheory.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - module CodingTheory ENV["HECKE_PRINT_BANNER"] = "false" @@ -22,29 +16,73 @@ export sizeof_perfect_code, sizeof_all_words, has_identity, has_identity_on_left include("utils.jl") # Abstract types -export FinitePolynomial, AbstractCode, NonStaticAbstractWord, AbstractWord, Word, - Codewords, Alphabet, Messages, CodeUniverse, CodeUniverseIterator, - UniverseParameters -export no_round, getindex, setindex!, firstindex, lastindex, size, length, rand, - gensym, genalphabet, eltype, isword, isabstractword +export FinitePolynomial, + AbstractCode, + NonStaticAbstractWord, + AbstractWord, + Word, + Codewords, + Alphabet, + Messages, + CodeUniverse, + CodeUniverseIterator, + UniverseParameters +export no_round, + getindex, + setindex!, + firstindex, + lastindex, + size, + length, + rand, + gensym, + genalphabet, + eltype, + isword, + isabstractword # RREF export rref, rref! # Bounds, Messages, Distance, and Primes -export rate, sphere_covering_bound, sphere_packing_bound, hamming_bound, - singleton_bound, gilbert_varshamov_bound, elias_bassalygo_bound, - plotkin_bound, construct_ham_matrix, isperfect, ishammingperfect, - isgolayperfect +export rate, + sphere_covering_bound, + sphere_packing_bound, + hamming_bound, + singleton_bound, + gilbert_varshamov_bound, + elias_bassalygo_bound, + plotkin_bound, + construct_ham_matrix, + isperfect, + ishammingperfect, + isgolayperfect export get_codewords_greedy, get_codewords_random, get_all_words, get_codewords -export hamming_distance, hamming_ball, code_distance, t_error_detecting, - t_error_correcting, find_error_detection_max, find_error_correction_max +export hamming_distance, + hamming_ball, + code_distance, + t_error_detecting, + t_error_correcting, + find_error_detection_max, + find_error_correction_max export isprimepower # Algebra -export FinitePolynomial, list_polys, multiplication_table, list_span, islinear, - isirreducible, normal_form!, normal_form, equivalent_code!, equivalent_code, - generator!, generator, parity_check, syndrome, isincode +export FinitePolynomial, + list_polys, + multiplication_table, + list_span, + islinear, + isirreducible, + normal_form!, + normal_form, + equivalent_code!, + equivalent_code, + generator!, + generator, + parity_check, + syndrome, + isincode # Powers export isprimepower, isperfectpower diff --git a/src/abstract_types.jl b/src/abstract_types.jl index 4068808..cfb77c4 100644 --- a/src/abstract_types.jl +++ b/src/abstract_types.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - """ ```julia abstract type FiniteField end @@ -23,7 +17,8 @@ abstract type AbstractCode end NonStaticAbstractWord{N, T} = Union{NTuple{N, T}, AbstractVector{T}, AbstractString} where {N, T} ``` """ -NonStaticAbstractWord{N, T} = Union{NTuple{N, T}, AbstractVector{T}, AbstractString} where {N, T} +NonStaticAbstractWord{N, T} = + Union{NTuple{N, T}, AbstractVector{T}, AbstractString} where {N, T} """ ```julia @@ -38,14 +33,13 @@ A Word is an `StaticArrays.MVector` which is efficient (like tuple) yet mutable. """ mutable struct Word{N, T} w::MVector{N, T} - + Word(w::NTuple{N, T}) where {N, T} = new{N, T}(MVector{N, T}(w)) - Word(w::AbstractVector{T}) where {T} = - (len = length(w); new{len, T}(MVector{len, T}(w))) - Word(w::AbstractString) = - MVector{length(w), eltype(w)}(collect(w)) - Word(i::T...) where {T} = - (len = length(i); new{len, T}(MVector{len, T}(i))) + function Word(w::AbstractVector{T}) where {T} + return (len = length(w); new{len, T}(MVector{len, T}(w))) + end + Word(w::AbstractString) = MVector{length(w), eltype(w)}(collect(w)) + Word(i::T...) where {T} = (len = length(i); new{len, T}(MVector{len, T}(i))) end # Indexing Interface @@ -63,7 +57,8 @@ Base.length(::Word{N, T}) where {N, T} = N AbstractWord{N} = Union{NonStaticAbstractWord{N, T}, MVector{T}} where {T} ``` """ -AbstractWord{N, T} = Union{Word{N, T}, NonStaticAbstractWord{N, T}, MVector{N, T}} where {N, T} +AbstractWord{N, T} = + Union{Word{N, T}, NonStaticAbstractWord{N, T}, MVector{N, T}} where {N, T} isword(w) = w isa Word isword(i...) = isword(i) @@ -77,10 +72,9 @@ Codewords{N} <: AbstractCode Simply a wrapper type for a vector of abstract words of length `N`. """ -Codewords{N} = Vector{AbstractWord{N, T}} where T +Codewords{N} = Vector{AbstractWord{N, T}} where {T} # Codewords{N} = Vector{Word{N, T}} where T - """ ```julia struct Alphabet <: AbstractCode @@ -95,7 +89,7 @@ Alphabet(Σ::AbstractArray) ``` A constructor method for the struct Alphabet. Takes an array of letters in the alphabet, and attempts to parse them as 64-bit Ints. - + --- ```julia @@ -105,11 +99,13 @@ Alphabet{N}(Σ::AbstractString) A constructor method for the struct Alphabet. Takes in a symbols and splits it into constituent characters. Those symbols are the letters in the alphabet. `N` is the number of letters in the alphabet """ -struct Alphabet{N} <: AbstractVector{Symbol} where N +struct Alphabet{N} <: AbstractVector{Symbol} where {N} Σ::AbstractVector{Symbol} - Alphabet(Σ::Union{Vector{T}, String}) where {T} = - (Σ_unique = Set(Σ); new{length(Σ_unique)}(ensure_symbolic(Σ_unique))) + function Alphabet(Σ::Union{Vector{T}, String}) where {T} + return (Σ_unique = Set(Σ); + new{length(Σ_unique)}(ensure_symbolic(Σ_unique))) + end end # end struct # Indexing Interface @@ -159,13 +155,14 @@ UniverseParameters(q::Int, n::Int) An inner constructor function on the structure `UniverseParameters`. """ struct UniverseParameters <: AbstractCode - Σ::Alphabet{N} where N + Σ::Alphabet{N} where {N} q::Int # size of alphabet n::Int # block length - + UniverseParameters(Σ::Alphabet{N}, n::Int) where {N} = new(Σ, N, n) UniverseParameters(Σ::Alphabet{N}, q::Int, n::Int) where {N} = new(Σ, q, n) - UniverseParameters(Σ::AbstractVector{T}, n::Int) where {T} = (Σ = Alphabet(Σ); new(Σ, length(Σ), n)) + UniverseParameters(Σ::AbstractVector{T}, n::Int) where {T} = (Σ = Alphabet(Σ); + new(Σ, length(Σ), n)) UniverseParameters(Σ::AbstractVector{T}, q::Int, n::Int) where {T} = new(Σ, q, n) UniverseParameters(q::Int, n::Int) = new(genalphabet(q), q, n) end @@ -187,7 +184,8 @@ Given universe parameters 𝒰 and a code C, return a tuple including - A random letter in the alphabet; and - A random index in the block length. """ -Base.rand(𝒰::UniverseParameters, C::AbstractArray) = tuple(Word(rand(C)), rand(𝒰.Σ), rand(1:𝒰.n)) +Base.rand(𝒰::UniverseParameters, C::AbstractArray) = + tuple(Word(rand(C)), rand(𝒰.Σ), rand(1:(𝒰.n))) """ ```julia @@ -219,41 +217,46 @@ struct CodeUniverseIterator <: AbstractCode function CodeUniverseIterator(𝒰::UniverseParameters) Σ = 𝒰.Σ - return Iterators.product(Vector{eltype(Σ)}[Σ for _ in 1:𝒰.n]...) + return Iterators.product(Vector{eltype(Σ)}[Σ for _ in 1:(𝒰.n)]...) end - CodeUniverseIterator(Σ::Union{Alphabet{N}, AbstractVector{T}}, q::Int, n::Int) where {N, T} = - CodeUniverseIterator(UniverseParameters(Σ, q, n)) - CodeUniverseIterator(Σ::Union{Alphabet{N}, AbstractVector{T}}, n::Int) where {N, T} = - CodeUniverseIterator(UniverseParameters(Σ, n)) - CodeUniverseIterator(q::Int, n::Int) = - CodeUniverseIterator(UniverseParameters(q, n)) + function CodeUniverseIterator( + Σ::Union{Alphabet{N}, AbstractVector{T}}, q::Int, n::Int + ) where {N, T} + return CodeUniverseIterator(UniverseParameters(Σ, q, n)) + end + function CodeUniverseIterator( + Σ::Union{Alphabet{N}, AbstractVector{T}}, n::Int + ) where {N, T} + return CodeUniverseIterator(UniverseParameters(Σ, n)) + end + CodeUniverseIterator(q::Int, n::Int) = CodeUniverseIterator(UniverseParameters(q, n)) end """ ```julia struct CodeUniverse <: AbstractCode ``` - + Defines a structure for the messages in the code. Parameters are the abstract array of messages `𝒰`, and the length of the messages `n`. ```julia CodeUniverse(𝒰::AbstractArray, Σ::Alphabet) ``` - + An inner constructor function on the structure `CodeUniverse`. """ struct CodeUniverse <: AbstractCode 𝒰::Codewords{M} where {M} - Σ::Alphabet{N} where N + Σ::Alphabet{N} where {N} q::Int # alphabet size n::Int # block length - + function CodeUniverse(𝒰::Codewords{M}, Σ::Alphabet{N}) where {N, M} message_length_error = "We have fixed the block length of each message. Please ensure all messages are of equal length." _allequal_length_(𝒰) || throw(error(message_length_error)) - - new(𝒰, Σ, length(Σ), length(first(𝒰))) + + return new(𝒰, Σ, length(Σ), length(first(𝒰))) end # end constructor function end diff --git a/src/algebra.jl b/src/algebra.jl index 3a2c0ef..95910ba 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - """ ```julia struct FinitePolynomial <: FiniteField @@ -20,13 +14,13 @@ FinitePolynomial(p::AbstractPolynomial, n::Int) A constructor method for `FinitePolynomial`. Takes in a polynomial `p` and a number `n`, and constructs a polynomial under modulo n. """ struct FinitePolynomial <: FiniteField - p::AbstractPolynomial - n::Int - - function FinitePolynomial(p::AbstractPolynomial, n::Int) - p = Polynomial(mod.(p.coeffs, n)) - new(p, n) - end + p::AbstractPolynomial + n::Int + + function FinitePolynomial(p::AbstractPolynomial, n::Int) + p = Polynomial(mod.(p.coeffs, n)) + return new(p, n) + end end """ @@ -76,7 +70,8 @@ Returns: - `Array`: An array of polynomials of degree less than n, under modulo m. """ function list_polys(n::Int, m::Int) - return collect(Polynomial(collect(t)) for t in Iterators.product([0:(m-1) for i in 1:n]...)) + return collect(Polynomial(collect(t)) for + t in Iterators.product([0:(m - 1) for i in 1:n]...)) end """ @@ -112,21 +107,21 @@ julia> julia> multiplication_table(2, 3) # multiplication table of all polynomia ``` """ function multiplication_table(degree::Int, modulo::Int) - polys = list_polys(degree, modulo) - number_of_polys = length(polys) - poly_matrix = Matrix{Polynomial}(undef, number_of_polys, number_of_polys) - - for i in 1:number_of_polys, j in 1:number_of_polys - poly_matrix[i,j] = mod(polys[i] * polys[j], modulo) - end + polys = list_polys(degree, modulo) + number_of_polys = length(polys) + poly_matrix = Matrix{Polynomial}(undef, number_of_polys, number_of_polys) + + for i in 1:number_of_polys, j in 1:number_of_polys + poly_matrix[i, j] = mod(polys[i] * polys[j], modulo) + end - return poly_matrix + return poly_matrix end -function __list_span_inner(modulo::Int, u̲::Vector{T}...) where T <: Number +function __list_span_inner(modulo::Int, u̲::Vector{T}...) where {T <: Number} n_vec = length(u̲) span = Vector{T}[zero(u̲[1])] - + for n in 1:n_vec new_span = copy(span) for v in span, λ in 1:(modulo - 1) @@ -137,8 +132,8 @@ function __list_span_inner(modulo::Int, u̲::Vector{T}...) where T <: Number end span = new_span end - - return span + + return span end """ @@ -175,7 +170,7 @@ julia> list_span([2, 1, 1], [1, 1, 1], 3) # list the span of two vectors modulo [0, 1, 1] ``` """ -list_span(args...) = __list_span_inner(last(args), args[1:end-1]...) +list_span(args...) = __list_span_inner(last(args), args[1:(end - 1)]...) """ ```julia @@ -203,36 +198,37 @@ false ``` """ function islinear(C::Vector, modulo::Int; verbose::Bool = false) - allequal_length(C) || return false # not all codes are of the same length - block_length = length(C[1]) - 𝟎 = fill(0, block_length) - - if 𝟎 ∉ C - verbose && println("The zero vector 0̲ is not in C.\n") - return false # the zero vector is not in the code - end - - for c̲ ∈ C - for λ in 0:modulo-1 - if mod.(λ*c̲, modulo) ∉ C - verbose && println(λ, " ⋅ ", c̲, " = ", mod.(λ*c̲, modulo), " ∉ C\n") - return false # this code isn't closed under scalar multiplication - end - end - - for c̲ ∈ C, c̲′ ∈ C - if c̲ ≠ c̲′ - if mod.(c̲ + c̲′, modulo) ∉ C - verbose && println(c̲, " + ", c̲′, " = ", mod.(c̲ + c̲′, modulo), " ∉ C\n") - return false # this code isn't closed under addition - end - end - end - end - - verbose && println() - - return true + allequal_length(C) || return false # not all codes are of the same length + block_length = length(C[1]) + 𝟎 = fill(0, block_length) + + if 𝟎 ∉ C + verbose && println("The zero vector 0̲ is not in C.\n") + return false # the zero vector is not in the code + end + + for c̲ in C + for λ in 0:(modulo - 1) + if mod.(λ * c̲, modulo) ∉ C + verbose && println(λ, " ⋅ ", c̲, " = ", mod.(λ * c̲, modulo), " ∉ C\n") + return false # this code isn't closed under scalar multiplication + end + end + + for c̲ in C, c̲′ in C + if c̲ ≠ c̲′ + if mod.(c̲ + c̲′, modulo) ∉ C + verbose && + println(c̲, " + ", c̲′, " = ", mod.(c̲ + c̲′, modulo), " ∉ C\n") + return false # this code isn't closed under addition + end + end + end + end + + verbose && println() + + return true end """ @@ -259,15 +255,15 @@ true ``` """ function isirreducible(f::Polynomial, modulo::Int) - deg = length(f.coeffs) - 1 - f = mod(f, deg) - polys = list_polys(deg, modulo) + deg = length(f.coeffs) - 1 + f = mod(f, deg) + polys = list_polys(deg, modulo) - for a in polys, b in polys - isequal(f, mod(a*b, modulo)) && return false - end - - return true + for a in polys, b in polys + isequal(f, mod(a * b, modulo)) && return false + end + + return true end """ @@ -299,7 +295,7 @@ julia> normal_form!([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2] ``` """ function normal_form!(M::AbstractArray{Int}, n::Int) - return rref!(M, n, colswap = false) + return rref!(M, n; colswap = false) end """ @@ -331,7 +327,7 @@ julia> normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], ``` """ function normal_form(M::AbstractArray{Int}, n::Int) - normal_form!(copy(M), n) + return normal_form!(copy(M), n) end """ @@ -363,7 +359,7 @@ julia> equivalent_code!([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 ``` """ function equivalent_code!(M::AbstractArray{Int}, n::Int) - return rref!(M, n, colswap=true) + return rref!(M, n; colswap = true) end """ @@ -394,7 +390,7 @@ julia> equivalent_code([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 ``` """ function equivalent_code(M::AbstractArray{Int}, n::Int) - equivalent_code!(copy(M), n) + return equivalent_code!(copy(M), n) end """ @@ -413,7 +409,7 @@ Returns: - `Matrix{Int}`: A generating matrix. """ function generator!(M::AbstractArray{Int}, n::Int; colswap::Bool = false) - return ifelse(colswap, equivalent_code!(M, n), normal_form!(M, n)) + return ifelse(colswap, equivalent_code!(M, n), normal_form!(M, n)) end """ @@ -432,7 +428,7 @@ Returns: - `Matrix{Int}`: A generating matrix. """ function generator(M::AbstractArray{Int}, n::Int; colswap::Bool = true) - generator!(copy(M), n, colswap=colswap) + return generator!(copy(M), n; colswap = colswap) end """ @@ -450,10 +446,11 @@ Returns: - `Matrix{Int}`: A parity check matrix. """ function parity_check(M::AbstractArray{Int}, n::Int) - has_identity_on_left(M) || throw(error("This matrix is not in normal form. Use normal_form or equivalent_code.")) - - minus_Dᵀ = mod.(-one(Int) .* transpose(M[:, (size(M, 1) + 1):end]), n) - return H = Int[minus_Dᵀ I(size(minus_Dᵀ, 1))] + has_identity_on_left(M) || + throw(error("This matrix is not in normal form. Use normal_form or equivalent_code.")) + + minus_Dᵀ = mod.(-one(Int) .* transpose(M[:, (size(M, 1) + 1):end]), n) + return H = Int[minus_Dᵀ I(size(minus_Dᵀ, 1))] end """ @@ -485,7 +482,7 @@ julia> syndrome([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 ``` """ function syndrome(v̲::Vector, Hᵀ::AbstractArray{Int}, n::Int) - return mod.(v̲' * Hᵀ, n) + return mod.(v̲' * Hᵀ, n) end """ @@ -512,5 +509,5 @@ true ``` """ function isincode(v̲::Vector, Hᵀ::AbstractArray{Int}, n::Int) - return iszero(syndrome(v̲, Hᵀ, n)) + return iszero(syndrome(v̲, Hᵀ, n)) end diff --git a/src/bounds.jl b/src/bounds.jl index 91a4cc7..c7df3b4 100644 --- a/src/bounds.jl +++ b/src/bounds.jl @@ -2,18 +2,20 @@ ```julia rate(q::Integer, M::Integer, n::Integer) -> Real ``` - + Calculate the rate of a code. That is, how efficient the code is. Parameters: + - `q::Integer`: the number of symbols in the code. - `M::Integer`: the size/number of elements in the code. - `n::Integer`: The word length. Returns: + - `Real`: Rate of the code. ---- +* * * ### Examples @@ -24,146 +26,164 @@ julia> rate(3, 5, 4) # the rate of the code which has 3 symbols, 5 words in the """ rate(q::T, M::T, n::T) where {T <: Integer} = log(q, M) / n -__spheres(q::T, n::T, r::T) where {T <: Integer} = sum(Integer[((big(q) - 1)^i) * binomial(big(n), big(i)) for i in 0:r]) -__sphere_bound(round_func::Function, q::T, n::T, d::T) where {T <: Integer} = round_func((big(q)^n) / __spheres(q, n, d)) +function __spheres(q::T, n::T, r::T) where {T <: Integer} + return sum(Integer[((big(q) - 1)^i) * binomial(big(n), big(i)) for i in 0:r]) +end +function __sphere_bound(round_func::Function, q::T, n::T, d::T) where {T <: Integer} + return round_func((big(q)^n) / __spheres(q, n, d)) +end """ ```julia sphere_covering_bound(q::Integer, n::Integer, d::Integer) -> Integer ``` - + Computes the sphere covering bound of a ``[n, d]_q``-code. Parameters: + - `q::Integer`: the number of symbols in the code. - `n::Integer`: the word length. - `d::Integer`: the distance of the code. - + Returns: + - `Integer`: the sphere covering bound. ---- +* * * ### Examples ```julia -julia> sphere_covering_bound(5,7,3) +julia> sphere_covering_bound(5, 7, 3) 215 ``` """ -sphere_covering_bound(q::T, n::T, d::T) where {T <: Integer} = __sphere_bound(ceil, q, n, d - 1) +sphere_covering_bound(q::T, n::T, d::T) where {T <: Integer} = + __sphere_bound(ceil, q, n, d - 1) """ ```julia sphere_packing_bound(q::Integer, n::Integer, d::Integer) -> Integer sphere_packing_bound(q::Integer, n::Integer, d::Integer, ::Rounding) -> Real ``` - + Computes the sphere packing bound of a ``[n, d]_q``-code. The sphere packing bound is also known as the hamming bound. You can use `hamming_bound` to compute the same thing. Parameters: + - `q::Integer`: the number of symbols in the code. - `n::Integer`: the word length. - `d::Integer`: the distance of the code. - - `::Rounding`: use the argument `no_round` in this position to preserve the rounding of the code — which usually by default rounds down. - + - `::Rounding`: use the argument `no_round` in this position to preserve the rounding of the code — which usually by default rounds down. + Returns: + - `Integer`: the sphere packing bound. ---- +* * * ### Examples ```julia -julia> sphere_packing_bound(5,7,3) +julia> sphere_packing_bound(5, 7, 3) 2693 ``` """ -sphere_packing_bound(q::T, n::T, d::T) where T <: Integer = - __sphere_bound(a -> floor(T, a), q, n, floor(T, (d - 1) / 2)) -sphere_packing_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer = - __sphere_bound(identity, q, n, floor(T, (d - 1) / 2)) -hamming_bound(q::T, n::T, d::T) where T <: Integer = - sphere_packing_bound(q, n, d) -hamming_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer = - sphere_packing_bound(q, n, d, no_round) +sphere_packing_bound(q::T, n::T, d::T) where {T <: Integer} = + __sphere_bound(a -> floor(T, a), q, n, floor(T, (d - 1) / 2)) +function sphere_packing_bound(q::T, n::T, d::T, ::Rounding) where {T <: Integer} + return __sphere_bound(identity, q, n, floor(T, (d - 1) / 2)) +end +hamming_bound(q::T, n::T, d::T) where {T <: Integer} = sphere_packing_bound(q, n, d) +function hamming_bound(q::T, n::T, d::T, ::Rounding) where {T <: Integer} + return sphere_packing_bound(q, n, d, no_round) +end """ ```julia sphere_packing_bound(q::Integer, n::Integer, d::Integer) -> Real ``` - + Computes the Singleton bound of a ``[n, d]_q``-code. Parameters: - `q::Integer`: the number of symbols in the code. - `n::Integer`: the word length. - `d::Integer`: the distance of the code. - + Returns: - `Real`: the Singleton bound. Can round down, as it is an equivalent to the Hamming bound in that it is an upper bound. """ # promote() # _T = typeof(T) -singleton_bound(q::T, n::T, d::T) where T <: Integer = - floor(T, float(big(q))^(big(n) - big(d) + 1)) -singleton_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer = - float(big(q))^(big(n) - big(d) + 1) - - -gilbert_varshamov_bound(q::T, n::T, d::T) where T <: Integer = - __sphere_bound(a -> floor(T, a), q, n, d - 1) -gilbert_varshamov_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer = - __sphere_bound(identity, q, n, d - 1) - -function __plotkin_bound_core(round_func::Function, q::T, n::T, d::T) where T <: Integer - if ! isequal(q, 2) - throw(error("The Plotkin bound only works for the binary code.")) - end - - if iseven(d) && 2d > n - return round_func((d) / (2d + 1 - n)) - elseif isodd(d) && 2d + 1 > n - return round_func((d + 1) / (2d + 1 - n)) - elseif iseven(d) - return T(4d) - # return A_2(2d, d) ≤ 4d - elseif isodd(d) - return T(4d + 4) - # return A_2(2d + 1, d) ≤ 4d + 4 - end +function singleton_bound(q::T, n::T, d::T) where {T <: Integer} + return floor(T, float(big(q))^(big(n) - big(d) + 1)) +end +function singleton_bound(q::T, n::T, d::T, ::Rounding) where {T <: Integer} + return float(big(q))^(big(n) - big(d) + 1) end -plotkin_bound(q::T, n::T, d::T) where T <: Integer = - __plotkin_bound_core(a -> floor(T, a), q, n, d) -plotkin_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer = - __plotkin_bound_core(identity, q, n, d, no_round) - -elias_bassalygo_bound(q::T, n::T, d::T) where T <: Integer = -elias_bassalygo_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer = - -function __johnson_bound_core(round_func::Function, q::T, n::T, d::T) where T <: Integer - if isinteger((d - 1) / 2) # is odd - t = T((d - 1) / 2) - __sphere_bound(round_func, q, n, t) # if d = 2t + 1 - elseif isinteger(d / 2) - t = T(d / 2) - __sphere_bound(round_func, q, n, t) - end +function gilbert_varshamov_bound(q::T, n::T, d::T) where {T <: Integer} + return __sphere_bound(a -> floor(T, a), q, n, d - 1) +end +function gilbert_varshamov_bound(q::T, n::T, d::T, ::Rounding) where {T <: Integer} + return __sphere_bound(identity, q, n, d - 1) +end + +function __plotkin_bound_core(round_func::Function, q::T, n::T, d::T) where {T <: Integer} + if !isequal(q, 2) + throw(error("The Plotkin bound only works for the binary code.")) + end + + if iseven(d) && 2d > n + return round_func((d) / (2d + 1 - n)) + elseif isodd(d) && 2d + 1 > n + return round_func((d + 1) / (2d + 1 - n)) + elseif iseven(d) + return T(4d) + # return A_2(2d, d) ≤ 4d + elseif isodd(d) + return T(4d + 4) + # return A_2(2d + 1, d) ≤ 4d + 4 + end +end + +function plotkin_bound(q::T, n::T, d::T) where {T <: Integer} + return __plotkin_bound_core(a -> floor(T, a), q, n, d) +end +function plotkin_bound(q::T, n::T, d::T, ::Rounding) where {T <: Integer} + return __plotkin_bound_core(identity, q, n, d, no_round) +end + +function elias_bassalygo_bound(q::T, n::T, d::T) where {T <: Integer} + return function elias_bassalygo_bound(q::T, n::T, d::T, ::Rounding) where {T <: Integer} + return function __johnson_bound_core( + round_func::Function, q::T, n::T, d::T + ) where {T <: Integer} + if isinteger((d - 1) / 2) # is odd + t = T((d - 1) / 2) + __sphere_bound(round_func, q, n, t) # if d = 2t + 1 + elseif isinteger(d / 2) + t = T(d / 2) + __sphere_bound(round_func, q, n, t) + end + end + end end @doc raw""" ```julia construct_ham_matrix(r::Int, q::Int) -> Matrix ``` - + Construct a Hamming parity-check matrix. Parameters: - `r::Int`: number of rows of a parity check matrix. - `q:::Int`: The size of the alphabet of the code. - + Returns: - `Matrix`: The Hamming matrix, denoted as ``\text{Ham}(r, q)`` @@ -182,11 +202,11 @@ julia> construct_ham_matrix(3,2) function construct_ham_matrix(r::Int, q::Int) ncols = Int(floor((q^r - 1) / (q - 1))) M = Matrix{Int}(undef, r, ncols) - + for i in 1:ncols - M[:, i] = reverse(digits(parse(Int, string(i, base = q)), pad = r), dims = 1) + M[:, i] = reverse(digits(parse(Int, string(i; base = q)); pad = r); dims = 1) end - + return M end @@ -194,19 +214,21 @@ end ```julia isperfect(n::Int, k::Int, d::Int, q::Int) -> Bool ``` - + Checks if a code is perfect. That is, checks if the number of words in the code is exactly the "Hamming bound", or the "Sphere Packing Bound". - + Parameters: + - `q:::Int`: The size of the alphabet of the code. - `n::Int`: The length of the words in the code (block length). - `d::Int`: The distance of the code. - `k::Int`: The dimension of the code. - + Returns: + - `Bool`: true or false ---- +* * * ### Examples @@ -215,26 +237,29 @@ julia> isperfect(11, 6, 5, 3) true ``` """ -function isperfect(n::T, k::T, d::T, q::T) where T <: Int - isprimepower(q) || throw(error("Cannot check if the code is perfect with q not a prime power.")) +function isperfect(n::T, k::T, d::T, q::T) where {T <: Int} + isprimepower(q) || + throw(error("Cannot check if the code is perfect with q not a prime power.")) M = q^k - + isequal(sphere_packing_bound(q, n, d), M) && return true - return false + return false end """ ```julia ishammingbound(r::Int, q::Int) -> Bool ``` - + Checks if the code is a perfect code that is of the form of a generalised Hamming code. - + Parameters: + - `r::Int`: number of rows of a parity check matrix. - `q::Int`: The size of the alphabet of the code. - + Returns: + - `Bool`: true or false """ function ishammingperfect(r::Int, q::Int) @@ -246,11 +271,11 @@ function ishammingperfect(r::Int, q::Int) # r is dim of dueal code; dim of code itself is block length minus r # println(n) # println((q^r - 1) / (q - 1)) - + isequal(n, (q^r - 1) / (q - 1)) && \ - isequal(d, 3) && \ - isequal(M, q^(((q^r - 1) / (q - 1)) - r)) && \ - return true + isequal(d, 3) && \ + isequal(M, q^(((q^r - 1) / (q - 1)) - r)) && \ + return true return false end @@ -259,19 +284,21 @@ end ishammingperfect(n::Int, k::Int, d::Int, q::Int) -> Bool ishammingperfect(q::Int, n::Int, d::Int) -> Bool ``` - + Checks if the code is a perfect code that is of the form of a generalised Hamming code. - + Parameters: + - `q:::Int`: The size of the alphabet of the code. - `n::Int`: The length of the words in the code (block length). - `d::Int`: The distance of the code. - `k::Int`: The dimension of the code. - + Returns: + - `Bool`: true or false ---- +* * * ### Examples @@ -280,47 +307,51 @@ julia> isgolayperfect(11, 6, 5, 3) # this is one of golay's perfect codes true ``` """ -function ishammingperfect(n::T, k::T, d::T, q::T) where T <: Int +function ishammingperfect(n::T, k::T, d::T, q::T) where {T <: Int} isprimepower(q) || return false - + M = q^k r = log(ℯ, ((n * log(ℯ, 1)) / (log(ℯ, 2))) + 1) / log(ℯ, 2) - - if isequal(n, (q^(r - 1)) / (q - 1)) && isequal(d, 3) && isequal(M, q^(((q^r - 1) / (q - 1)) - r)) + + if isequal(n, (q^(r - 1)) / (q - 1)) && + isequal(d, 3) && + isequal(M, q^(((q^r - 1) / (q - 1)) - r)) return true end - + return false end function ishammingperfect(q::Int, n::Int, d::Int) - isprimepower(q) || return false # we are working in finite fields, so q must be a prime power - d ≠ 3 && return false - - r = 1 - while ((q^r - 1) / (q - 1)) < n - r = r + 1 - end - - return ifelse(isequal(((q^r - 1) / (q - 1)), n), true, false) + isprimepower(q) || return false # we are working in finite fields, so q must be a prime power + d ≠ 3 && return false + + r = 1 + while ((q^r - 1) / (q - 1)) < n + r = r + 1 + end + + return ifelse(isequal(((q^r - 1) / (q - 1)), n), true, false) end """ ```julia isgolayperfect(n::Int, k::Int, d::Int, q::Int) -> Bool ``` - + Golay found two perfect codes. `isgolayperfect` checks if a code of block length n, distance d, alphabet size q, and dimension k, is a perfect code as described by Golay. Parameters: + - `n::Int`: The block length of words in the code (e.g., word "abc" has block length 3). - `k::Int`: The dimension of the code. - `d::Int`: The distance of the code (i.e., the minimum distance between codewords in the code). - `q::Int`: An Int that is a prime power. The modulus of the finite field. - + Returns: + - `Bool`: true or false. ---- +* * * ### Examples @@ -329,8 +360,8 @@ julia> isgolayperfect(11, 6, 5, 3) # this is one of golay's perfect codes true ``` """ -function isgolayperfect(n::T, k::T, d::T, q::T) where T <: Int - isprimepower(q) || false # we are working in finite fields, so q must be a prime power +function isgolayperfect(n::T, k::T, d::T, q::T) where {T <: Int} + isprimepower(q) || false # we are working in finite fields, so q must be a prime power M = q^k (isequal(q, 2) && isequal(n, 23) && isequal(d, 7) && isequal(M, 2^12)) && return true (isequal(q, 3) && isequal(n, 11) && isequal(d, 5) && isequal(M, 3^6)) && return true diff --git a/src/distance.jl b/src/distance.jl index 7989a71..98f5e5a 100755 --- a/src/distance.jl +++ b/src/distance.jl @@ -1,25 +1,20 @@ - -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - """ ```julia hamming_distance(w₁, w₂) -> Int ``` The Hamming distance of two words is the number of changes that need to be made to each letter in the word for the words to be the same. This does not work for words of unequal length. - + Parameters: + - `w₁`: A word. - `w₂`: Another word. - + Returns: + - `Int`: the number of changes needing to be made to one word for it to be identical to the other. ---- +* * * ### Examples @@ -28,27 +23,29 @@ julia> hamming_distance("ABC", "BBC") # computes the hamming distance 1 ``` """ -function hamming_distance(w₁::T, w₂::T) where T <: AbstractWord - if ! isequal(length(w₁), length(w₂)) +function hamming_distance(w₁::T, w₂::T) where {T <: AbstractWord} + if !isequal(length(w₁), length(w₂)) throw(error("Cannot compute Hamming Distance on strings of unequal length.")) end - + distance = 0 - - for (s₁, s₂) in zip(w₁, w₂) - if s₁ ≠ s₂ - distance += 1 - end - end - + + for (s₁, s₂) in zip(w₁, w₂) + if s₁ ≠ s₂ + distance += 1 + end + end + return distance end -function __hamming_space(relation::Function, Σⁿ::AbstractArray{T}, w::AbstractArray, e::Int) where T <: AbstractWord - e < 0 && throw(error("e (the ball/sphere \"radius\") must be a non-negative number.")) - Σⁿ, w = deepsym(Σⁿ), ensure_symbolic(w) - - return eltype(Σⁿ)[v for v in deepsym(Σⁿ) if relation(hamming_distance(w, v), e)] +function __hamming_space( + relation::Function, Σⁿ::AbstractArray{T}, w::AbstractArray, e::Int +) where {T <: AbstractWord} + e < 0 && throw(error("e (the ball/sphere \"radius\") must be a non-negative number.")) + Σⁿ, w = deepsym(Σⁿ), ensure_symbolic(w) + + return eltype(Σⁿ)[v for v in deepsym(Σⁿ) if relation(hamming_distance(w, v), e)] end """ @@ -59,14 +56,16 @@ hamming_ball(Σⁿ::AbstractArray, w::Vector, e::Int) -> Vector{Vector} Get the codewords of radius ``e`` of a ball centered at word ``w``. That is, all words whose distance from w is less than or equal to the radius. Parameters: + - `Σⁿ::AbstractArray`: An array of words in the code. - `w::Vector`: A word. - `e::Int`: The radius of the ball. - + Returns: + - AbstractArray: The list of words in Σⁿ whose distance from w is less than or equal to e. Returns an array of array of symbols. ---- +* * * ### Examples @@ -78,7 +77,7 @@ julia> hamming_ball([[1, 0, 1], [0, 1, 1], [1, 0, 0]], [1, 0, 0], 2) # given a l ``` """ hamming_ball(Σⁿ::AbstractArray{T}, w::Vector{S}, e::Int) where {T <: AbstractWord, S} = - __hamming_space(≤, Σⁿ, w, e) + __hamming_space(≤, Σⁿ, w, e) """ ```julia @@ -88,44 +87,48 @@ hamming_sphere(Σⁿ::AbstractArray, w::Vector, e::Int) -> Vector{Vector} Get the codewords of radius e of a sohere centered at word w. That is, all words whose distance from w is exactly equal to to the radius. Parameters: + - `Σⁿ::AbstractArray`: An array of words in the code. - `w::Vector`: A word. - `e::Int`: The radius of the ball. - + Returns: + - `AbstractArray`: The list of words in Σⁿ whose distance from w is exactly equal to e. Returns an array of array of symbols. """ hamming_sphere(Σⁿ::AbstractArray{T}, w::Vector{S}, e::Int) where {T <: AbstractWord, S} = - __hamming_space(isequal, Σⁿ, w, e) + __hamming_space(isequal, Σⁿ, w, e) """ ```julia code_distance(C::AbstractArray) -> Int ``` - + Finds the distance of the code. That is, given a code C, finds the minimum distance between any two words in the code, which are not the same. (Find the minimum hamming distance in the code for all unique letters). Parameters: + - `C::AbstractArray`: An array of words in the code. - + Return: + - `Int`: the distance of the code. """ -function code_distance(C::AbstractArray{T}) where T <: AbstractWord - isempty(C) && return nothing - C = deepsym(C) - min_distance = nothing - - for c in C, c′ in C - if c ≠ c′ - ham_dist = hamming_distance(c, c′) - if isnothing(min_distance) || ham_dist < min_distance - min_distance = ham_dist - end - end - end - - return min_distance +function code_distance(C::AbstractArray{T}) where {T <: AbstractWord} + isempty(C) && return nothing + C = deepsym(C) + min_distance = nothing + + for c in C, c′ in C + if c ≠ c′ + ham_dist = hamming_distance(c, c′) + if isnothing(min_distance) || ham_dist < min_distance + min_distance = ham_dist + end + end + end + + return min_distance end """ @@ -136,15 +139,17 @@ code_distance!(C::AbstractArray{T}, w::T) -> Int A wrapper to get the code distance after pushing a word to the code. *This directly changes the matrix M. Use `code_distance` for a non-mutating version of this function.* Parameters: + - `C::AbstractArray`: An array of words in the code. - `w`: A word to be appended to C. - + Returns: + - `Int`: The code distance after adding w to C. """ -function code_distance!(C::AbstractArray{T}, w::T) where T <: AbstractWord - push!(C, w) - return code_distance(C) +function code_distance!(C::AbstractArray{T}, w::T) where {T <: AbstractWord} + push!(C, w) + return code_distance(C) end """ @@ -155,40 +160,44 @@ code_distance(C::AbstractArray{T}, w::T) -> Int A wrapper to get the code distance after pushing a word to the code. Parameters: + - `C::AbstractArray`: An array of words in the code. - `w`: A word to be appended to C. - + Returns: + - `Int`: The code distance after adding w to C. ---- +* * * ### Examples ```julia -julia> code_distance([[0,0,0,0,0],[1,0,1,0,1],[0,1,0,1,0],[1,1,1,1,1]]) # gets the minimum distance between two vectors in an array of vectors +julia> code_distance([[0, 0, 0, 0, 0], [1, 0, 1, 0, 1], [0, 1, 0, 1, 0], [1, 1, 1, 1, 1]]) # gets the minimum distance between two vectors in an array of vectors 2 ``` """ -function code_distance(C::AbstractArray{T}, w::T) where T <: AbstractWord - code_distance!(copy(C), w) +function code_distance(C::AbstractArray{T}, w::T) where {T <: AbstractWord} + return code_distance!(copy(C), w) end """ ```julia t_error_detecting(C::AbstractArray{T}, t::Int) -> Bool ``` - + Check if a given code C can detect t many errors. - + Parameters: + - `C::AbstractArray`: An array of words in the code. - `t::Int`: The number of errors you want to check that the code can detect. - + Returns: + - `Bool`: Yes, C can detect t errors, or no it cannot (true of false). ---- +* * * ### Examples @@ -197,74 +206,87 @@ julia> t_error_detecting([[1, 0, 1], [0, 1, 1], [1, 0, 0]], 3) false ``` """ -function t_error_detecting(C::AbstractArray{T}, t::Int) where T <: Union{AbstractArray{Int}, AbstractWord} - code_distance(C) ≥ t + 1 && return true - return false +function t_error_detecting( + C::AbstractArray{T}, t::Int +) where {T <: Union{AbstractArray{Int}, AbstractWord}} + code_distance(C) ≥ t + 1 && return true + return false end """ ```julia t_error_correcting(C::AbstractArray{T}, t::Int) -> Bool ``` - + Check if a given code C can correct t many errors. - + Parameters: + - `C::AbstractArray`: An array of words in the code. - `t::Int`: The number of errors you want to check that the code can correct. - + Returns: + - Bool: Yes, C can correct t errors, or no it cannot (true of false). """ -function t_error_correcting(C::AbstractArray{T}, t::Int) where T <: Union{AbstractArray{Int}, AbstractWord} - code_distance(C) ≥ 2*t + 1 && return true - return false +function t_error_correcting( + C::AbstractArray{T}, t::Int +) where {T <: Union{AbstractArray{Int}, AbstractWord}} + code_distance(C) ≥ 2 * t + 1 && return true + return false end - """ ```julia find_error_detection_max(C::AbstractArray{T}, modulo::Int) -> Int ``` - + Finds the greatest number t such that C is error t error detecting. - + Parameters: + - `C::AbstractArray`: An array of words in the code. - `moldulo::Int`: The modulus of the finite field. The upper bound of t. - + Returns: + - `Int`: The maximum number t such that the code is t error detecting. ---- +* * * ```julia julia> find_error_detection_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) 1 ``` """ -function find_error_detection_max(C::AbstractArray{T}, modulo::Int) where T <: Union{AbstractArray{Int}, AbstractWord} - for t in (modulo - 1):-1:0 - t_error_detecting(C, t) && return t - end +function find_error_detection_max( + C::AbstractArray{T}, modulo::Int +) where {T <: Union{AbstractArray{Int}, AbstractWord}} + for t in (modulo - 1):-1:0 + t_error_detecting(C, t) && return t + end end """ ```julia find_error_correction_max(C::AbstractArray{T}, modulo::Int) -> Int ``` - + Finds the greatest number t such that C is error t error correcting. - + Parameters: + - `C::AbstractArray`: An array of words in the code. - `moldulo::Int`: The modulus of the finite field. The upper bound of t. - + Returns: + - `Int`: The maximum number t such that the code is t error correcting. """ -function find_error_correction_max(C::AbstractArray{T}, modulo::Int) where T <: Union{AbstractArray{Int}, AbstractWord} - for t in modulo-1:-1:0 - t_error_correcting(C, t) && return t - end +function find_error_correction_max( + C::AbstractArray{T}, modulo::Int +) where {T <: Union{AbstractArray{Int}, AbstractWord}} + for t in (modulo - 1):-1:0 + t_error_correcting(C, t) && return t + end end diff --git a/src/levenshtein.jl b/src/levenshtein.jl index 850296c..787ab08 100755 --- a/src/levenshtein.jl +++ b/src/levenshtein.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - # Adapted from Levenshtein.jl """ @@ -15,17 +9,18 @@ levenshtein( target::AbstractString, deletion_cost::R, insertion_cost::S, - substitution_cost::T) -> Integer + substitution_cost::T, +) -> Integer levenshtein!( source::AbstractString, target::AbstractString, deletion_cost::R, insertion_cost::S, - substitution_cost::T, - costs::Matrix = Array{promote_type(R, S, T)}(undef, 2, length(target) + 1) + substitution_cost::T; + costs::Matrix = Array{promote_type(R, S, T)}(undef, 2, length(target) + 1), ) -> Integer ``` - + Computes the Levenshtein distance. *These methods are adapted from Levenshtein.jl, by Roger Tu.* @@ -43,9 +38,8 @@ function levenshtein( target::AbstractString, deletion_cost::R, insertion_cost::S, - substitution_cost::T + substitution_cost::T, ) where {R <: Real, S <: Real, T <: Real} - return levenshtein!(target, source, insertion_cost, deletion_cost, substitution_cost) end @@ -55,14 +49,15 @@ function levenshtein!( deletion_cost::R, insertion_cost::S, substitution_cost::T, - costs::Matrix = Array{promote_type(R, S, T)}(undef, 2, length(target) + 1) + costs::Matrix = Array{promote_type(R, S, T)}(undef, 2, length(target) + 1), ) where {R <: Real, S <: Real, T <: Real} + cost_type = promote_type(R, S, T) - cost_type = promote_type(R,S,T) - if length(source) < length(target) # Space complexity of function = O(length(target)) - return levenshtein!(target, source, insertion_cost, deletion_cost, substitution_cost, costs) + return levenshtein!( + target, source, insertion_cost, deletion_cost, substitution_cost, costs + ) else if iszero(length(target)) return length(source) * deletion_cost @@ -70,25 +65,25 @@ function levenshtein!( old_cost_index = 1 new_cost_index = 2 costs[old_cost_index, 1] = 0 - + for i in 1:length(target) costs[old_cost_index, i + 1] = i * insertion_cost end i = 0 - + for r in source i += 1 # Delete i characters from source to get empty target costs[new_cost_index, 1] = i * deletion_cost j = 0 - + for c in target j += 1 deletion = costs[old_cost_index, j + 1] + deletion_cost insertion = costs[new_cost_index, j] + insertion_cost substitution = costs[old_cost_index, j] - + if r ≠ c substitution += substitution_cost end @@ -100,7 +95,7 @@ function levenshtein!( end new_cost_index = old_cost_index - + return costs[new_cost_index, length(target) + 1] end end diff --git a/src/messages.jl b/src/messages.jl index 9840db4..3596983 100755 --- a/src/messages.jl +++ b/src/messages.jl @@ -1,9 +1,3 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - include("distance.jl") include("powers.jl") @@ -14,17 +8,17 @@ push_if_allowed!(C::AbstractArray{T}, w::T, d::Int) Takes in an array and a word. As long as the word does not mean that the distance is smaller than d, we add w to the array. If we are successful in doing this, return true. Otherwise, return false. *This is a mutating function.* """ -function push_if_allowed!(C::AbstractArray{T}, w::T, d::Int) where T <: AbstractWord - isempty(C) && (push!(C, w); return true) - - for c in C - if hamming_distance(c, w) < d - return false - end - end - - push!(C, w) - return true +function push_if_allowed!(C::AbstractArray{T}, w::T, d::Int) where {T <: AbstractWord} + isempty(C) && (push!(C, w); return true) + + for c in C + if hamming_distance(c, w) < d + return false + end + end + + push!(C, w) + return true end """ @@ -34,17 +28,19 @@ push_if_allowed!(C::AbstractArray{T}, C′::AbstractArray{T}, w::T, d::Int) Takes in two arrays, A and B. If w is allowed in C given distance d, push to C′. If we are successful in doing this, return true. Otherwise, return false. *This is a mutating function.* """ -function push_if_allowed!(C::AbstractVector{T}, C′::AbstractVector{T}, w::T, d::Int) where T <: AbstractWord - isempty(C) && (push!(C′, w); return true) - - for c in C - if hamming_distance(c, w) < d - return false - end - end - - push!(C′, w) - return true +function push_if_allowed!( + C::AbstractVector{T}, C′::AbstractVector{T}, w::T, d::Int +) where {T <: AbstractWord} + isempty(C) && (push!(C′, w); return true) + + for c in C + if hamming_distance(c, w) < d + return false + end + end + + push!(C′, w) + return true end """ @@ -54,37 +50,42 @@ replace_if_allowed!(C::AbstractArray, d::Int, w, w′) -> Bool Takes in an array and a word. As long as the word does not mean that the distance is smaller than d, we replace a with b in the array. Replaces and returns true if allowed; otherwise returns false. *This is a mutating function.* """ -function replace_if_allowed!(C::AbstractVector{T}, d::Int, ws::Pair) where T <: AbstractWord - w, w′ = ws - - for c in C - if hamming_distance(c, w′) < d - return false - end - end - - replace!(C, ws) - return true +function replace_if_allowed!( + C::AbstractVector{T}, d::Int, ws::Pair +) where {T <: AbstractWord} + w, w′ = ws + + for c in C + if hamming_distance(c, w′) < d + return false + end + end + + replace!(C, ws) + return true end -function replace_if_allowed!(C::AbstractVector{T}, d::Int, w::T, w′::T) where T <: AbstractWord - return replace_if_allowed!(C, d, Pair(w, w′)) +function replace_if_allowed!( + C::AbstractVector{T}, d::Int, w::T, w′::T +) where {T <: AbstractWord} + return replace_if_allowed!(C, d, Pair(w, w′)) end -_mutate_codeword(w::Word{N, T}, i::Int, a::T) where {T, N} = - setindex!(w, a, i) - +_mutate_codeword(w::Word{N, T}, i::Int, a::T) where {T, N} = setindex!(w, a, i) + """ ```julia -mutate_codeword(w::NonStaticAbstractWord{N, T}, n::Int, i::Int, a::T) where {T, N} -> MVector{N, T} +mutate_codeword(w::NonStaticAbstractWord{N, T}, n::Int, i::Int, a::T) where {T, N} -> + MVector{N, T} mutate_codeword(w::Word{N, T}, n::Int, i::Int, a::T) where {T, N} -> MVector{N, T} ``` Mutates the word w, which is an `MVector` of length N, changing its iᵗʰ index to a. """ mutate_codeword(w::NonStaticAbstractWord{N, T}, i::Int, a::T) where {T, N} = - _mutate_codeword(Word{N, T}(w), i, a) -mutate_codeword(w::Word{N, T}, n::Int, i::Int, a::T) where {T, N} = - _mutate_codeword(w, i, a) + _mutate_codeword(Word{N, T}(w), i, a) +function mutate_codeword(w::Word{N, T}, n::Int, i::Int, a::T) where {T, N} + return _mutate_codeword(w, i, a) +end """ ```julia @@ -94,20 +95,22 @@ get_all_words(Σ::AbstractArray, q::Int, n::Int) -> Codewords{M} get_all_words(Σ::AbstractArray, n::Int) -> Codewords{M} get_all_words(q::Int, n::Int) -> Codewords{M} ``` - + Get the universe of *all* codewords of a given alphabet. The alphabet will be uniquely generated if none is given. - + Parameters: + - `Σ::AbstractArray`: The alphabet allowed. - `q::Int`: The size of the alphabet. - `n::Int`: The (fixed) length of the words in the code. - `d::Int`: The minimum distance between words in the code. - `𝒰::AbstractArray`: The universe of all codewords of q many letters of block length n. - + Returns: + - `Codewords{M}`: An array of codewords, each of length `M`. Each codewords is a tuple, and each character in said word is a symbol. ---- +* * * ### Examples @@ -119,13 +122,10 @@ julia> get_all_words(2, 2) # all words of block length 2 using 2 unique symbols ``` """ get_all_words(Σ::Alphabet{N}, q::Int, n::Int) where {N} = - collect(CodeUniverseIterator(UniverseParameters(Σ, q, n))) -get_all_words(Σ::Alphabet{N}, n::Int) where {N} = - get_all_words(Σ, length(Set(Σ)), n) # if alphabet is given, then q is the length of that alphabet -get_all_words(Σ::AbstractArray, n::Int) = - get_all_words(Alphabet(Σ), length(Set(Σ)), n) # if alphabet is given, then q is the length of that alphabet -get_all_words(q::Int, n::Int) = - get_all_words(genalphabet(q), q, n) # generate symbols if no alphabet is given + collect(CodeUniverseIterator(UniverseParameters(Σ, q, n))) +get_all_words(Σ::Alphabet{N}, n::Int) where {N} = get_all_words(Σ, length(Set(Σ)), n) # if alphabet is given, then q is the length of that alphabet +get_all_words(Σ::AbstractArray, n::Int) = get_all_words(Alphabet(Σ), length(Set(Σ)), n) # if alphabet is given, then q is the length of that alphabet +get_all_words(q::Int, n::Int) = get_all_words(genalphabet(q), q, n) # generate symbols if no alphabet is given """ ```julia @@ -135,63 +135,83 @@ get_codewords_greedy(Σ::Alphabet{N}, n::Int, d::Int) -> Codewords{M} get_codewords_greedy(q::Int, n::Int, d::Int) -> Codewords{M} get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int) -> Codewords{M} get_codewords_greedy(Σ::AbstractArray, n::Int, d::Int) -> Codewords{M} -get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray) -> Codewords{M} +get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray) -> + Codewords{M} ``` - + Search through the universe of all codewords and find a code of block length n and distance d, using the alphabet Σ. The alphabet will be uniquely generated if none is given. This uses a greedy algorithm, simply iterating through all words (see above) and choosing them if they fit in the code. In some cases the greedy algorithm is the best, but in others it is very much not. - + Parameters: + - `𝒰::UniverseParameters`: The parameters of the universe of all codewords of q many letters of block length n. - `Σ::AbstractArray`: The alphabet allowed. - `q::Int`: The size of the alphabet. - `n::Int`: The (fixed) length of the words in the code. - `d::Int`: The minimum distance between words in the code. - + Returns: + - `Codewords{M}`: An array of codewords, each of length `M`. Each codewords is a tuple, and each character in said word is a symbol. """ function get_codewords_greedy(𝒰::UniverseParameters, d::Int) - C = eltype(𝒰)[] - - for wᵢ in CodeUniverseIterator(𝒰) - push_if_allowed!(C, wᵢ, d) - end - - return C + C = eltype(𝒰)[] + + for wᵢ in CodeUniverseIterator(𝒰) + push_if_allowed!(C, wᵢ, d) + end + + return C +end +function get_codewords_greedy(Σ::Alphabet{N}, q::Int, n::Int, d::Int) where {N} + return get_codewords_greedy(UniverseParameters(Σ, q, n), d) +end +function get_codewords_greedy(Σ::Alphabet{N}, n::Int, d::Int) where {N} + return get_codewords_greedy(UniverseParameters(Σ, n), d) +end +function get_codewords_greedy(q::Int, n::Int, d::Int) + return get_codewords_greedy(UniverseParameters(q, n), d) +end +function get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int) + return get_codewords_greedy(Alphabet(Σ), q, n, d) +end +function get_codewords_greedy(Σ::AbstractArray, n::Int, d::Int) + return get_codewords_greedy(UniverseParameters(Σ, n), d) +end +function get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray) + return get_codewords_greedy(Alphabet(Σ), q, n, d, 𝒰) +end + +function argmaxminima(A::AbstractArray; dims::Int) + return getindex(argmin(A; dims = dims), argmax(argmin(A; dims = dims))) +end +function maxminima(A::AbstractArray; dims::Int) + return getindex(minimum(A; dims = dims), maximum(minimum(A; dims = dims))) +end +function argminmaxima(A::AbstractArray; dims::Int) + return getindex(argmax(A; dims = dims), argmin(argmax(A; dims = dims))) +end +function minmaxima(A::AbstractArray; dims::Int) + return getindex(maximum(A; dims = dims), minimum(maximum(A; dims = dims))) end -get_codewords_greedy(Σ::Alphabet{N}, q::Int, n::Int, d::Int) where {N} = - get_codewords_greedy(UniverseParameters(Σ, q, n), d) -get_codewords_greedy(Σ::Alphabet{N}, n::Int, d::Int) where {N} = - get_codewords_greedy(UniverseParameters(Σ, n), d) -get_codewords_greedy(q::Int, n::Int, d::Int) = - get_codewords_greedy(UniverseParameters(q, n), d) -get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int) = - get_codewords_greedy(Alphabet(Σ), q, n, d) -get_codewords_greedy(Σ::AbstractArray, n::Int, d::Int) = - get_codewords_greedy(UniverseParameters(Σ, n), d) -get_codewords_greedy(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray) = - get_codewords_greedy(Alphabet(Σ), q, n, d, 𝒰) - - -argmaxminima(A::AbstractArray; dims::Int) = getindex(argmin(A, dims = dims), argmax(argmin(A, dims = dims))) -maxminima(A::AbstractArray; dims::Int) = getindex(minimum(A, dims = dims), maximum(minimum(A, dims = dims))) -argminmaxima(A::AbstractArray; dims::Int) = getindex(argmax(A, dims = dims), argmin(argmax(A, dims = dims))) -minmaxima(A::AbstractArray; dims::Int) = getindex(maximum(A, dims = dims), minimum(maximum(A, dims = dims))) """ ```julia get_codewords_random(𝒰::UniverseParameters, d::Int; m::Int = 1000) -> Codewords{M} -get_codewords_random(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=1000) -> Codewords{M} -get_codewords_random(Σ::Alphabet{N}, n::Int, d::Int; m::Int=1000) -> Codewords{M} -get_codewords_random(q::Int, n::Int, d::Int; m::Int=1000) -> Codewords{M} -get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=1000) -> Codewords{M} -get_codewords_random(Σ::AbstractArray, n::Int, d::Int; m::Int=1000) -> Codewords{M} -get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=1000) -> Codewords{M} +get_codewords_random(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int = 1000) -> Codewords{M} +get_codewords_random(Σ::Alphabet{N}, n::Int, d::Int; m::Int = 1000) -> Codewords{M} +get_codewords_random(q::Int, n::Int, d::Int; m::Int = 1000) -> Codewords{M} +get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int = 1000) -> + Codewords{M} +get_codewords_random(Σ::AbstractArray, n::Int, d::Int; m::Int = 1000) -> Codewords{M} +get_codewords_random( + Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int = 1000 +) -> Codewords{M} ``` Search through the universe of all codewords at random and find a code of block length n and distance d, using the alphabet Σ. The alphabet will be uniquely generated if none is given. This is a cleverer algorithm than the greedy algorithm. Increasing the `m` keyword argument arbitrarily _should_ produce a maximal code, as for each codeword it chooses, it collects a list of `m` many random words, and chooses the best one from that intermediate list. Parameters: + - `Σ::AbstractArray`: The alphabet allowed. - `q::Int`: The size of the alphabet. - `n::Int`: The (fixed) length of the words in the code. @@ -199,44 +219,54 @@ Parameters: - `𝒰::AbstractArray`: The universe of all codewords of q many letters of block length n. Returns: - - `Codewords{M}`: An array of codewords, each of length `M`. Each codewords is a tuple, and each character in said word is a symbol. + + - `Codewords{M}`: An array of codewords, each of length `M`. Each codewords is a tuple, and each character in said word is a symbol. # get a random word in the code start """ function get_codewords_random(𝒰::UniverseParameters, d::Int; m::Int = 1000) - C = eltype(𝒰)[] - - starting_word = rand(𝒰) # get a random word in the code start - push!(C, starting_word) - - for _ in 1:length(𝒰) - C′ = eltype(𝒰)[] - # while length(C′) < m - for _ in 1:m - wᵣ = rand(𝒰) - push_if_allowed!(C, C′, wᵣ, d) # if allowed in C, push to C′ - # wᵣ ∉ C′ && push_if_allowed!(C, C′, wᵣ, d) - end - isempty(C′) && break - # [push_if_allowed!(C, C′, w, d) for _ in 1:m] - distances = Int[hamming_distance(wᵢ, wⱼ) for wᵢ in C, wⱼ in C′] - best_word = getindex(C′, getindex(argmaxminima(distances, dims = 1), 2)) - push!(C, best_word) - end - - return C + C = eltype(𝒰)[] + + starting_word = rand(𝒰) # get a random word in the code start + push!(C, starting_word) + + for _ in 1:length(𝒰) + C′ = eltype(𝒰)[] + # while length(C′) < m + for _ in 1:m + wᵣ = rand(𝒰) + push_if_allowed!(C, C′, wᵣ, d) # if allowed in C, push to C′ + # wᵣ ∉ C′ && push_if_allowed!(C, C′, wᵣ, d) + end + isempty(C′) && break + # [push_if_allowed!(C, C′, w, d) for _ in 1:m] + distances = Int[hamming_distance(wᵢ, wⱼ) for wᵢ in C, wⱼ in C′] + best_word = getindex(C′, getindex(argmaxminima(distances; dims = 1), 2)) + push!(C, best_word) + end + + return C +end +function get_codewords_random( + Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int = 1000 +) where {N} + return get_codewords_random(UniverseParameters(Σ, q, n), d; m = m) +end +function get_codewords_random(Σ::Alphabet{N}, n::Int, d::Int; m::Int = 1000) where {N} + return get_codewords_random(UniverseParameters(Σ, n), d; m = m) +end +function get_codewords_random(q::Int, n::Int, d::Int; m::Int = 1000) + return get_codewords_random(UniverseParameters(q, n), d; m = m) +end +function get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int = 1000) + return get_codewords_random(Alphabet(Σ), q, n, d; m = m) +end +function get_codewords_random(Σ::AbstractArray, n::Int, d::Int; m::Int = 1000) + return get_codewords_random(UniverseParameters(Σ, n), d; m = m) +end +function get_codewords_random( + Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int = 1000 +) + return get_codewords_random(Alphabet(Σ), q, n, d, 𝒰; m = m) end -get_codewords_random(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=1000) where {N} = - get_codewords_random(UniverseParameters(Σ, q, n), d, m=m) -get_codewords_random(Σ::Alphabet{N}, n::Int, d::Int; m::Int=1000) where {N} = - get_codewords_random(UniverseParameters(Σ, n), d, m=m) -get_codewords_random(q::Int, n::Int, d::Int; m::Int=1000) = - get_codewords_random(UniverseParameters(q, n), d, m=m) -get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=1000) = - get_codewords_random(Alphabet(Σ), q, n, d, m=m) -get_codewords_random(Σ::AbstractArray, n::Int, d::Int; m::Int=1000) = - get_codewords_random(UniverseParameters(Σ, n), d, m=m) -get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=1000) = - get_codewords_random(Alphabet(Σ), q, n, d, 𝒰, m=m) - # using Mmap # function get_codewords_random_mmap(mmappath::AbstractString, 𝒰::UniverseParameters, d::Int) @@ -270,18 +300,32 @@ get_codewords_random(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractAr """ ```julia -get_codewords(𝒰::UniverseParameters, d::Int; m::Int=10, m_random::Int = 1000) -> Codewords{M} -get_codewords(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M} -get_codewords(q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M} -get_codewords(q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M} -get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M} -get_codewords(Σ::AbstractArray, n::Int, d::Int; m::Int=10, m_random::Int=1000) -> Codewords{M} -get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=10, m_random::Int=1000) -> Codewords{M} +get_codewords(𝒰::UniverseParameters, d::Int; m::Int = 10, m_random::Int = 1000) -> + Codewords{M} +get_codewords(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int = 10, m_random::Int = 1000) -> + Codewords{M} +get_codewords(q::Int, n::Int, d::Int; m::Int = 10, m_random::Int = 1000) -> Codewords{M} +get_codewords(q::Int, n::Int, d::Int; m::Int = 10, m_random::Int = 1000) -> Codewords{M} +get_codewords( + Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int = 10, m_random::Int = 1000 +) -> Codewords{M} +get_codewords(Σ::AbstractArray, n::Int, d::Int; m::Int = 10, m_random::Int = 1000) -> + Codewords{M} +get_codewords( + Σ::AbstractArray, + q::Int, + n::Int, + d::Int, + 𝒰::AbstractArray; + m::Int = 10, + m_random::Int = 1000, +) -> Codewords{M} ``` Use function `get_codewords_random` `m` many times (with `get_codewords_random(..., m = m_random)`), and `get_codewords_greedy`. Return the code with the greatest number of words. The alphabet will be uniquely generated if none is given. You can omit Σ and 𝒰. You can omit q if Σ is given. - + Parameters: + - `Σ::AbstractArray`: The alphabet allowed. - `q::Int`: The size of the alphabet. - `n::Int`: The (fixed) length of the words in the code. @@ -289,15 +333,16 @@ Parameters: - `𝒰::AbstractArray`: The universe of all codewords of q many letters of block length n. - `m::Int` (kwarg): Try a random code m many times. - `m_random::Int` (kwarg): The number of possible words `get_codewords_random` chooses from for _each_ word it selects. - + Returns: + - `Codewords{M}`: An array of codewords, each of length `M`. Each codewords is a tuple, and each character in said word is a symbol. !!! note *If you are looking for a _maximal_ code, this is likely the function you need.* Increasing `m` and `m_random` arbitrarily should ensure a maximal code—_however_, that computing power/time in not always possible, as it requires a lot of RAM to store certain codes in memeory. Efforts are being made to make this process better by using memory-mapped filed instead of storing codewords in RAM, but this will make it much slower as well. Help with this would be much appreciated. ---- +* * * ### Examples @@ -315,41 +360,60 @@ julia> get_codewords(["a", "b", "c"], 3, 2) # get codewords of block length 3 wi (:b, :b, :b) ``` """ -function get_codewords(𝒰::UniverseParameters, d::Int; m::Int=10, m_random::Int = 1000) - code_size = 0 - C = eltype(𝒰)[] - - - for _ in 1:m - random_code = get_codewords_random(𝒰, d; m = m_random) - random_size = length(random_code) - if random_size > code_size - code_size = random_size - C = random_code - end - end - - greedy_code = get_codewords_greedy(𝒰, d) - greedy_size = length(greedy_code) - if greedy_size > code_size - code_size = greedy_size - C = greedy_code - end - - return C +function get_codewords(𝒰::UniverseParameters, d::Int; m::Int = 10, m_random::Int = 1000) + code_size = 0 + C = eltype(𝒰)[] + + for _ in 1:m + random_code = get_codewords_random(𝒰, d; m = m_random) + random_size = length(random_code) + if random_size > code_size + code_size = random_size + C = random_code + end + end + + greedy_code = get_codewords_greedy(𝒰, d) + greedy_size = length(greedy_code) + if greedy_size > code_size + code_size = greedy_size + C = greedy_code + end + + return C +end +function get_codewords( + Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int = 10, m_random::Int = 1000 +) where {N} + return get_codewords(UniverseParameters(Σ, q, n), d; m = m, m_random = m_random) +end +function get_codewords( + Σ::Alphabet{N}, n::Int, d::Int; m::Int = 10, m_random::Int = 1000 +) where {N} + return get_codewords(UniverseParameters(Σ, n), d; m = m, m_random = m_random) +end +function get_codewords(q::Int, n::Int, d::Int; m::Int = 10, m_random::Int = 1000) + return get_codewords(UniverseParameters(q, n), d; m = m, m_random = m_random) +end +function get_codewords( + Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int = 10, m_random::Int = 1000 +) + return get_codewords(Alphabet(Σ), q, n, d; m = m, m_random = m_random) +end +function get_codewords(Σ::AbstractArray, n::Int, d::Int; m::Int = 10, m_random::Int = 1000) + return get_codewords(UniverseParameters(Σ, n), d; m = m, m_random = m_random) +end +function get_codewords( + Σ::AbstractArray, + q::Int, + n::Int, + d::Int, + 𝒰::AbstractArray; + m::Int = 10, + m_random::Int = 1000, +) + return get_codewords(Alphabet(Σ), q, n, d, 𝒰; m = m, m_random = m_random) end -get_codewords(Σ::Alphabet{N}, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) where {N} = - get_codewords(UniverseParameters(Σ, q, n), d, m=m, m_random=m_random) -get_codewords(Σ::Alphabet{N}, n::Int, d::Int; m::Int=10, m_random::Int=1000) where {N} = - get_codewords(UniverseParameters(Σ, n), d, m=m, m_random=m_random) -get_codewords(q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) = - get_codewords(UniverseParameters(q, n), d, m=m, m_random=m_random) -get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int; m::Int=10, m_random::Int=1000) = - get_codewords(Alphabet(Σ), q, n, d, m=m, m_random=m_random) -get_codewords(Σ::AbstractArray, n::Int, d::Int; m::Int=10, m_random::Int=1000) = - get_codewords(UniverseParameters(Σ, n), d, m=m, m_random=m_random) -get_codewords(Σ::AbstractArray, q::Int, n::Int, d::Int, 𝒰::AbstractArray; m::Int=10, m_random::Int=1000) = - get_codewords(Alphabet(Σ), q, n, d, 𝒰, m=m, m_random=m_random) """ ```julia @@ -357,33 +421,35 @@ get_codewords(G::AbstractArray, m::Int) -> Codewords{M} ``` Get codewords of a code from the _generating matrix_ under a finite field of modulo `m`. Precisely, computes all linear combinations of the rows of the generating matrix. - + Parameters: + - `G::AbstractArray`: A matrix of Ints which generates the code. - `m::Int`: The bounds of the finite field (i.e., the molulus you wish to work in). - + Returns: + - `Codewords{M}`: An array of codewords, each of length `M`. Each codewords is a tuple, and each character in said word is a symbol. """ function get_codewords(G::AbstractArray, m::Int) - codewords = Vector() - rows = Vector(undef, size(G, 2)) - - for i in 1:size(G, 1) - rows[i] = [G[i, j] for j in 1:size(G, 2)] - end - - for c in Base.Iterators.product([0:m-1 for _ in 1:size(G, 1)]...) - word = Ref(c[1]) .* rows[1] - - for i in 2:size(G, 1) - word = mod.(word .+ (Ref(c[i]) .* rows[i]), m) - end - - push!(codewords, word) - end - - return codewords + codewords = Vector() + rows = Vector(undef, size(G, 2)) + + for i in 1:size(G, 1) + rows[i] = [G[i, j] for j in 1:size(G, 2)] + end + + for c in Base.Iterators.product([0:(m - 1) for _ in 1:size(G, 1)]...) + word = Ref(c[1]) .* rows[1] + + for i in 2:size(G, 1) + word = mod.(word .+ (Ref(c[i]) .* rows[i]), m) + end + + push!(codewords, word) + end + + return codewords end # function obtain_maximal_code(𝒰::UniverseParameters, d::Int) diff --git a/src/powers.jl b/src/powers.jl index 8a9408a..924077c 100644 --- a/src/powers.jl +++ b/src/powers.jl @@ -82,7 +82,7 @@ end function Ω(n) n == fmpz(0) && return 0 __isprime(n) && return fmpz(1) - + return sum(e for (__, e) ∈ __factors(n)) end diff --git a/src/rref.jl b/src/rref.jl index 80918fb..8e25505 100644 --- a/src/rref.jl +++ b/src/rref.jl @@ -1,66 +1,70 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - #= Gauss-Jordan elimination over finite fields, modulo n =# -function swaprows!(A::Matrix{T}, i::Int, j::Int) where T +function swaprows!(A::Matrix{T}, i::Int, j::Int) where {T} for k in axes(A, 2) A[i, k], A[j, k] = A[j, k], A[i, k] end - + return A end -function swapcols!(A::Matrix{T}, i::Int, j::Int) where T +function swapcols!(A::Matrix{T}, i::Int, j::Int) where {T} for k in axes(A, 1) A[k, i], A[k, j] = A[k, j], A[k, i] end - + return A end """ ```julia -rref!(A::Matrix{Int}, n::Int; colswap::Bool=false, verbose::Bool=false, vverbose::Bool=false) -> Matrix{Int} +rref!( + A::Matrix{Int}, + n::Int; + colswap::Bool = false, + verbose::Bool = false, + vverbose::Bool = false, +) -> Matrix{Int} ``` - + Performs Gauss-Jordan Elimination on a matrix `A`. *This directly changes the matrix A. Use `rref` for a non-mutating version of this function.* Parameters: + - `A::Matrix{Int}`: A matrix of Ints you wish to perform Gauss-Jordan elimiation on. - `n::Int`: The modulus of the finite field you are working under. - `colswap::Bool` (kwarg): Whether or not you allow for column swapping. - `verbose::Bool` (kwarg): Print the row operations. - `vverbose::Bool` (kwarg): Print the intermediate matrices of the algorithm. - + Returns: + - `Matrix{Int}`: a matrix in row echelon form. ---- +* * * ### Examples ```julia -julia> rref!([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=true) # gauss-jordan elimitation modulo 5 with column swapping +julia> rref!([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5; colswap = true) # gauss-jordan elimitation modulo 5 with column swapping 3×6 Array{Int64,2}: 1 0 0 3 2 2 0 1 0 2 1 1 0 0 1 0 0 4 +``` # Rule 1: Swap zero rows if out of order and ensure leading ones cascade down diagonally. ``` """ -function rref!(A::Matrix{Int}, +function rref!( + A::Matrix{Int}, n::Int; colswap::Bool = false, verbose::Bool = false, - vverbose::Bool = false) - + vverbose::Bool = false, +) i = j = 1 - + while i ≤ size(A, 1) && j ≤ size(A, 2) # Rule 1: Swap zero rows if out of order and ensure leading ones cascade down diagonally. s = findfirst(!iszero, A[i:end, :]) @@ -69,46 +73,49 @@ function rref!(A::Matrix{Int}, if verbose println("r$(i) ⟷ r$(i + s[1] - 1)") if vverbose - displaymatrix(A); println() + displaymatrix(A) + println() end end - + # Rule 2: Normalize each row - s = findfirst(!iszero, A[i,:]) + s = findfirst(!iszero, A[i, :]) isnothing(s) && break α = invmod(A[i, s], n) - A[i,:] .= mod.(A[i,:] .* α, n) + A[i, :] .= mod.(A[i, :] .* α, n) if verbose - if ! iszero(α) + if !iszero(α) isone(α) || println("r$(i) ⟵ r$(i) × $(α)") - if vverbose && ! isone(α) - displaymatrix(A); println() + if vverbose && !isone(α) + displaymatrix(A) + println() end end end - + # Rule 3: Subtract it from the others - s = findfirst(!iszero, A[i,:]) + s = findfirst(!iszero, A[i, :]) isnothing(s) && break for k in axes(A, 1) if i ≠ k β = A[k, s] A[k, :] .= mod.(A[k, :] .- β .* A[i, :], n) if verbose - if ! iszero(β) + if !iszero(β) if isone(β) println("r$(k) ⟵ r$(k) - r$(i)") else println("r$(k) ⟵ r$(k) - $(β)r$(i)") end if vverbose - displaymatrix(A); println() + displaymatrix(A) + println() end end end end end - + # Rule 4: Swap columns if needed (and allowed) if colswap if iszero(A[i, j]) @@ -117,17 +124,18 @@ function rref!(A::Matrix{Int}, if verbose println("c$(j) ⟷ c$(s)") if vverbose - displaymatrix(A); println() + displaymatrix(A) + println() end end end end - + # increment row and column counter respectively i += 1 j += 1 end - + if verbose println() end @@ -137,36 +145,45 @@ end """ ```julia -rref(A::Matrix{Int}, n::Int; colswap::Bool=false, verbose::Bool=false, vverbose::Bool=false) -> Matrix{Int} +rref( + A::Matrix{Int}, + n::Int; + colswap::Bool = false, + verbose::Bool = false, + vverbose::Bool = false, +) -> Matrix{Int} ``` - + Performs Gauss-Jordan Elimination on a matrix `A`. Parameters: + - `A::Matrix{Int}`: A matrix of Ints you wish to perform Gauss-Jordan elimiation on. - `n::Int`: The modulus of the finite field you are working under. - `colswap::Bool` (kwarg): Whether or not you allow for column swapping. - `verbose::Bool` (kwarg): Print the row operations. - `vverbose::Bool` (kwarg): Print the intermediate matrices of the algorithm. - + Returns: + - `Matrix{Int}`: a matrix in row echelon form. ---- +* * * ### Examples ```julia -julia> rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=true) # gauss-jordan elimitation modulo 5 with column swapping +julia> rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5; colswap = true) # gauss-jordan elimitation modulo 5 with column swapping 3×6 Array{Int64,2}: 1 0 0 3 2 2 0 1 0 2 1 1 0 0 1 0 0 4 ``` """ -rref(A::Matrix{Int}, +rref( + A::Matrix{Int}, n::Int; colswap::Bool = false, verbose::Bool = false, - vverbose::Bool = false - ) = rref!(copy(A), n::Int; colswap = colswap, verbose = verbose, vverbose = vverbose) + vverbose::Bool = false, +) = rref!(copy(A), n::Int; colswap = colswap, verbose = verbose, vverbose = vverbose) diff --git a/src/utils.jl b/src/utils.jl index 9cbbf04..f38e9cb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,21 +1,18 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - """ ```julia displaymatrix(M::AbstractArray) ``` - + Displays a matrix `M` in a compact form from the terminal. """ function displaymatrix(M::AbstractArray) - return show(IOContext(stdout, :limit => true, :compact => true, :short => true), "text/plain", M); print("\n") + return show( + IOContext(stdout, :limit => true, :compact => true, :short => true), "text/plain", M + ) + return print("\n") end -_Iterable{T} = Union{AbstractArray{T}, NTuple{N, T}} where N +_Iterable{T} = Union{AbstractArray{T}, NTuple{N, T}} where {N} """ ```julia @@ -25,13 +22,13 @@ allequal_length(a, b...) -> Bool Check that all elements in a list are of equal length. """ -@inline function allequal_length(A::_Iterable{T}) where T +@inline function allequal_length(A::_Iterable{T}) where {T} length(A) < 2 && return true - + @inbounds for i in 2:length(A) isequal(length(first(A)), length(A[i])) || return false end - + return true end @inline allequal_length(a::T...) where {T} = allequal_length(a) @@ -44,19 +41,19 @@ allequal(a, b...) -> Bool Check that all elements in a list are equal to each other. """ -@inline function allequal(A::_Iterable{T}) where T +@inline function allequal(A::_Iterable{T}) where {T} length(A) < 2 && return true - + @inbounds for i in 2:length(A) first(A) ≠ A[i] && return false end - + return true end @inline allequal(a::T...) where {T} = allequal(a) @inline function allequal2(A) - return !A=A[1:1]==A∪A + return !A = A[1:1] == A ∪ A end """ @@ -67,14 +64,14 @@ aredistinct(a, b...) -> Bool Check that all elements in a list are distinct from every other element in the list. """ -@inline function aredistinct(A::_Iterable{T}) where T +@inline function aredistinct(A::_Iterable{T}) where {T} length(A) < 2 && return true - - while ! iszero(length(A)) + + while !iszero(length(A)) a = pop!(A) a ∉ A || return false end - + return true end @inline aredistinct(a::T...) where {T} = aredistinct(a) @@ -87,11 +84,11 @@ arelessthan(x::Number, a, b...) -> Bool Check that all elements in a list are less than a given x. """ -@inline function arelessthan(x::Number, A::_Iterable{T}) where T +@inline function arelessthan(x::Number, A::_Iterable{T}) where {T} @inbounds for a in A a < x || return false end - + return true end @inline arelessthan(x::Number, a::Number...) = arelessthan(x, a) @@ -104,11 +101,11 @@ areequalto(x::Number, a, b...) -> Bool Check that all elements in a list are equal to a given x. """ -@inline function areequalto(x, A::_Iterable{T}) where T +@inline function areequalto(x, A::_Iterable{T}) where {T} @inbounds for a in A a != x || return false end - + return true end @inline areequalto(x, a::Number...) = areequalto(x, a) @@ -132,7 +129,7 @@ Returns the type of the inner-most element in a nested array structure. """ deepeltype(a) = deepeltype(typeof(a)) deepeltype(::Type{T}) where {T <: _Iterable{R}} where {R} = deepeltype(eltype(T)) -deepeltype(::Type{T}) where T = T +deepeltype(::Type{T}) where {T} = T """ ```julia @@ -156,26 +153,32 @@ ensure_symbolic(Σ) = ensure_symbolic!(copy(Σ)) ```julia _has_identity(M::Matrix, n::Union{T, Base.OneTo{T}, AbstractRange{T}}) where {T <: Integer} ``` + Inner function on `has_identity` function. Will return a tuple of: + - Whether or not the matrix has an identity; - Where that identity matrix starts; and - How big the identity matrix is. """ -function _has_identity(M::Matrix, n::Union{T, Base.OneTo{T}, AbstractRange{T}}) where {T <: Integer} - for ID_size in n # iterate through possible identity sizes; try identity matrixes from the maximum size down - n isa Integer || (isone(ID_size) && continue) # skip I(1) unless that is the integer specifically given - max_starting_row, max_starting_col = size(M) .- (ID_size - 1) - for i in CartesianIndices(M) # iterate through all possible starting positions - m = M[i] - (isone(m) || iszero(m)) || continue # skip starting positions that are invalid - _row, _col = Tuple(i) - if _row ≤ max_starting_row && _col ≤ max_starting_col # ensure the matrix at position i has enough space for an identity - isequal(M[_row:(_row + ID_size - 1), _col:(_col + ID_size - 1)], I(ID_size)) && return true, (_row, _col), ID_size - end - end - end - - return false, ntuple(_ -> nothing, ndims(M)), nothing +function _has_identity( + M::Matrix, n::Union{T, Base.OneTo{T}, AbstractRange{T}} +) where {T <: Integer} + for ID_size in n # iterate through possible identity sizes; try identity matrixes from the maximum size down + n isa Integer || (isone(ID_size) && continue) # skip I(1) unless that is the integer specifically given + max_starting_row, max_starting_col = size(M) .- (ID_size - 1) + for i in CartesianIndices(M) # iterate through all possible starting positions + m = M[i] + (isone(m) || iszero(m)) || continue # skip starting positions that are invalid + _row, _col = Tuple(i) + if _row ≤ max_starting_row && _col ≤ max_starting_col # ensure the matrix at position i has enough space for an identity + isequal( + M[_row:(_row + ID_size - 1), _col:(_col + ID_size - 1)], I(ID_size) + ) && return true, (_row, _col), ID_size + end + end + end + + return false, ntuple(_ -> nothing, ndims(M)), nothing end """ @@ -185,10 +188,10 @@ has_identity(M::Matrix, n::Integer) -> Bool ``` Checks if a matrix has an identity in it. If given a number `n`, the function will specifically check if it has an identity matrix _of size n_ in `M`. - + See `CodingTheory._has_identity` for more information. ---- +* * * ### Examples @@ -211,7 +214,14 @@ julia> C = [-96 -66 20 1 0 0; -65 59 -82 0 1 0; -16 87 -113 0 0 1] -65 59 -82 0 1 0 -16 87 -113 0 0 1 -julia> D = [78 -99 125 -123 -111 -71 17; -115 78 40 -88 81 -40 78; -99 126 -54 1 0 0 24; -55 88 42 0 1 0 -8; 119 55 2 0 0 1 -92; -40 -21 -89 -79 59 -44 9] +julia> D = [ + 78 -99 125 -123 -111 -71 17 + -115 78 40 -88 81 -40 78 + -99 126 -54 1 0 0 24 + -55 88 42 0 1 0 -8 + 119 55 2 0 0 1 -92 + -40 -21 -89 -79 59 -44 9 + ] 6×7 Array{Int64,2}: 78 -99 125 -123 -111 -71 17 -115 78 40 -88 81 -40 78 @@ -239,10 +249,8 @@ julia> has_identity(D, 4) false ``` """ -has_identity(M::Matrix) = - _has_identity(M, reverse(minimum(axes(M))))[1] -has_identity(M::Matrix, n::T) where {T <: Integer} = - _has_identity(M, n)[1] +has_identity(M::Matrix) = _has_identity(M, reverse(minimum(axes(M))))[1] +has_identity(M::Matrix, n::T) where {T <: Integer} = _has_identity(M, n)[1] """ ```julia @@ -251,7 +259,7 @@ has_identity_on_left(M::Matrix) -> Bool Checks if the left-hand side of a matrix contains an identity matrix. ---- +* * * ### Examples @@ -286,9 +294,11 @@ sizeof_perfect_code(q::Number, n::Number, d::Number) -> Number Calculates the number of gigabytes required to store a perfect code of parameters q, n, and d. """ function sizeof_perfect_code(q::Int, n::Int, d::Int) - return (sizeof(ntuple(_ -> gensym(), n)) * hamming_bound(big.([q, n, d])...)) / (2^30) + return (sizeof(ntuple(_ -> gensym(), n)) * hamming_bound(big.([q, n, d])...)) / (2^30) +end +function sizeof_perfect_code(q::Number, n::Number, d::Number) + return sizeof_perfect_code(round.(BigInt, [q, n, d])...) end -sizeof_perfect_code(q::Number, n::Number, d::Number) = sizeof_perfect_code(round.(BigInt, [q, n, d])...) """ ```julia @@ -298,6 +308,6 @@ sizeof_all_words(q::Number, n::Number) -> Number Calculates the number of gigabytes required to store all unique words of length n from an alphabet of size q. """ function sizeof_all_words(q::Int, n::Int) - return (sizeof(ntuple(_ -> gensym(), n)) * big(q)^n) / (2^30) + return (sizeof(ntuple(_ -> gensym(), n)) * big(q)^n) / (2^30) end sizeof_all_words(q::Number, n::Number) = sizeof_all_words(round.(BigInt, (q, n))...) diff --git a/src/wip/primepowers.jl b/src/wip/primepowers.jl index 33753c7..0a0f34e 100644 --- a/src/wip/primepowers.jl +++ b/src/wip/primepowers.jl @@ -4,18 +4,18 @@ using Primes: primes, factor, isprime ### Helpers -nroot(x::Number, n::Number) = x^(1/n) +nroot(x::Number, n::Number) = x^(1 / n) function factors(n::Integer) f = [one(n)] - for (p,e) in factor(n) - f = reduce(vcat, [f*p^j for j in 1:e], init=f) + for (p, e) in factor(n) + f = reduce(vcat, [f * p^j for j in 1:e]; init = f) end return ifelse(isone(length(f)), [one(n), n], sort!(f)) end - + """ ord(q::Integer, n::Integer) -> Integer @@ -27,7 +27,7 @@ function ord(q::Integer, n::Integer) for (p, e) in factor(q) m = p^e - t = div(q, p) * (p-1) + t = div(q, p) * (p - 1) for f in factors(t) if isone(powermod(n, f, q)) res = lcm(res, f) @@ -39,8 +39,6 @@ function ord(q::Integer, n::Integer) return res end - - ### Main code #= \begin{enumerate} \item For each prime power $q$ such that $2^q\leq n$, write down a positive integer $r_q$ such that if $n$ is a $q$\th power, then $n=r^q_q$.\item Find a finite coprime set $P$ of integers larger than 1 such that each of $n,r_2,r_3,r_4, r_5,r_7,\ldots$ is a product of powers of elements of $P$. (In this paper, ``coprime'' means ``pairwise coprime.'')\item Factor $n$ as $\Pi_{p\in P}p^{n_p}$, and compute $k=\text{gcd}n_p:p\in P$. \end{enumerate}=# @@ -59,7 +57,7 @@ From \texttt{Powers.pdf}, section 21:\begin{quotation} \textbf{Lemma 21.1.} \emph{If $2i\equiv 2j \mod 2^{b+1}$, and $b\geq 1$, then $i^2\equiv j^2\mod w^{b+1}.$}\par \textbf{2-adoc approximate powers.} Fix positive integers $k$ and $b$. For any integer $m$ define $\text{pow}_{2,b}(m,k)=m^k\mod 2^b$. See section 6 for methods of computing $\text{pow}_{2, b}$ without many multiplications. As I compute $\text{pow}_{2,b}(m, k)$, I keep track of an ``overflow bit'' to figure out whether $m^k\mod 2^b=m^k$. \end{quotation} - + From \texttt{Powers.pdf}, section 6:\begin{quotation} Let $r$ be a positive floating-point number, and let $k$ be positive integers. Then $\text{pow}_{b}(r, k)$, the \textbf{$\mathbf{b}$-bit approximate $\mathbf{k}th$ power of $\mathbf{r}$} is a floating-point approximation to $r^k$. In this section, I show how to compute $\text{pow}_b(r,k)$ in $M$-time at most $P(k)M(b), where $P(k)\leq 2\lfloor\text{lg} k\rfloor$.\par Define $P(k)$ for $k\geq 1$ as follows: $P(1)=0;P(2k)=P(k)+1;P(2k+1)=P(2k)+1$.\par @@ -73,11 +71,10 @@ From \texttt{Powers.pdf}, section 6:\begin{quotation} \end{quotation} From \texttt{Powers.pdf}, section 22:\begin{quotation} - + \end{quotation} =# - lg(x) = log2(x) #= @@ -87,52 +84,52 @@ Design an algorithm that divides a nonnegative n-place integer $(u_1u_2\ldots u_ =# """ - divrem(b::Integer, u::Integer, v::Real) + divrem(b::Integer, u::Integer, v::Real) This algorithm divides a non-negative n-place integer (u1 u2 ... un) by v, where v is a single-precision number (i.e., 0 < v < b), producing the quotient (w1 w2 ... wn) and remainder r. This algorithm was written by Donald Knuth, in The Art of Computer Programming, semi-numerical algorithms, 2nd edition, exercise 4.3.1, 16. """ function Base.divrem(b::Integer, u::Integer, v::Real) - r, j = 0, 1 - n = ndigits(u) #; base = b - u̲ = round.(reverse(digits(u)), base = b) - w̲ = Vector{Int}(undef, n) - - while j ≤ n - w̲[j] = floor(Int, (r*b + u̲[j]) / v) - r = mod(r*b + u̲[j], v) - j += 1 - end - - # w = sum(w̲[i] * 10 ^ (length(w̲) - i) for i in eachindex(w̲)) - w = parse(Int, join(string.(w̲))) - # w = foldr((i, j) -> 10 * j + i, reverse(w̲)) - # - # n = 0 - # reduce(w̲) do i, j - # 10 * j + i - # # global n = n + ndigits(w̲[i]) - # end - # n = 0 - # w = sum(eachindex(w̲)) do i - # w̲[i] * 10 ^ (length(w̲) - i + n) - # global n = n + ndigits(w̲[i]) - # end - # n = 0 - - return w, r + r, j = 0, 1 + n = ndigits(u) #; base = b + u̲ = round.(reverse(digits(u)), base = b) + w̲ = Vector{Int}(undef, n) + + while j ≤ n + w̲[j] = floor(Int, (r * b + u̲[j]) / v) + r = mod(r * b + u̲[j], v) + j += 1 + end + + # w = sum(w̲[i] * 10 ^ (length(w̲) - i) for i in eachindex(w̲)) + w = parse(Int, join(string.(w̲))) + # w = foldr((i, j) -> 10 * j + i, reverse(w̲)) + # + # n = 0 + # reduce(w̲) do i, j + # 10 * j + i + # # global n = n + ndigits(w̲[i]) + # end + # n = 0 + # w = sum(eachindex(w̲)) do i + # w̲[i] * 10 ^ (length(w̲) - i + n) + # global n = n + ndigits(w̲[i]) + # end + # n = 0 + + return w, r end """ - div(b::Integer, u::Integer, v::Real) + div(b::Integer, u::Integer, v::Real) This algorithm divides a non-negative n-place integer (u1 u2 ... un) by v, where v is a single-precision number (i.e., 0 < v < b), producing the quotient (w1 w2 ... wn) and ommiting the remained. To return the quotient and the remainder, see `divrem`. """ Base.div(b::Integer, u::Integer, v::Real) = first(divrem(b, u, v)) """ - rem(b::Integer, u::Integer, v::Real) + rem(b::Integer, u::Integer, v::Real) This algorithm divides a non-negative n-place integer (u1 u2 ... un) by v, where v is a single-precision number (i.e., 0 < v < b), producing a remainder r. To return the quotient and the remainder, see `divrem`. """ @@ -143,13 +140,13 @@ In this section, I define truncation to $b$ bits, written $\text{trunc}_b$, and =# """ - trunc(b::Integer, n::Real) - + trunc(b::Integer, n::Real) + Truncates a number `n` to `b` bits. """ function Base.trunc(b::Integer, n::Real) - return n & (1 << b - 1) - # return div(b, n, 1) + return n & (1 << b - 1) + # return div(b, n, 1) end #= @@ -162,55 +159,55 @@ Given a positive floating-point number r and two positive integers b, k, to prim =# """ - pow(b::Integer, r::Real, k::Integer) - + pow(b::Integer, r::Real, k::Integer) + Given a positive floating-point number r and two positive integers b and k, pow(b, r, k) computes the b-bit approximate kth power of r. pow_b(r, k) ≤ r^k < pow_b(r, k)(1 + 2^(1 - b))^(2k - 1). """ function pow(b::Integer, r::Real, k::Integer) - if r ≤ 0 || b ≤ 0 || k ≤ 0 - throw(error("r, b, and k must be positive numbers to compute pow_b(r, k).")) - end - isone(k) && return trunc(b, r) - iszero(mod(k, 2)) && return trunc(b, pow(b, r, k ÷ 2)^2) # use \div because k is even - return trunc(b, pow(b, r, k - 1) * trunc(b, r)) + if r ≤ 0 || b ≤ 0 || k ≤ 0 + throw(error("r, b, and k must be positive numbers to compute pow_b(r, k).")) + end + isone(k) && return trunc(b, r) + iszero(mod(k, 2)) && return trunc(b, pow(b, r, k ÷ 2)^2) # use \div because k is even + return trunc(b, pow(b, r, k - 1) * trunc(b, r)) end """ - pow2(b::Integer, r::Real, k::Integer) - + pow2(b::Integer, r::Real, k::Integer) + Given a positive floating-point number r and two positive integers b and k, pow2(b, r, k) computes the 2-adic b-bit approximate kth power of r. """ function pow2(b::Integer, m::Integer, k::Integer) - if k ≤ 0 || b ≤ 0 - throw(error("k, and b must be positive integers to find a 2-adic approximate power.")) - end - - return mod(m^k, 2^b) + if k ≤ 0 || b ≤ 0 + throw(error("k, and b must be positive integers to find a 2-adic approximate power.")) + end + + return mod(m^k, 2^b) end """ - tentative_check_kth_roots(n::Integer, x::Integer, k::Integer) + tentative_check_kth_roots(n::Integer, x::Integer, k::Integer) A straight forward algorithm for checking whether x^k = n. """ function tentative_check_kth_roots(n::Integer, x::Integer, k::Integer) - if n ≤ 0 || x ≤ 0 || k ≤ 0 - throw(error("n, x, and k must be positive integers to compute check whether n=x^k.")) - end - - f = floor(Integer, lg(2*n)) - - isone(x) && return ifelse(isone(n), true, false) # return ifelse(isone(n), 0, 2) - - b = 1 - while true - # r = pow(x, b, k) # 2, b - r = pow2(b, x, k) - isequal(mod(n, 2^b), r) || return false # 2 - b ≥ f && return ifelse(isequal(r, x^k), true, false) # ifelse(isequal(r, x^k), 0, 2) - - b = min(2*b, f) - end + if n ≤ 0 || x ≤ 0 || k ≤ 0 + throw(error("n, x, and k must be positive integers to compute check whether n=x^k.")) + end + + f = floor(Integer, lg(2 * n)) + + isone(x) && return ifelse(isone(n), true, false) # return ifelse(isone(n), 0, 2) + + b = 1 + while true + # r = pow(x, b, k) # 2, b + r = pow2(b, x, k) + isequal(mod(n, 2^b), r) || return false # 2 + b ≥ f && return ifelse(isequal(r, x^k), true, false) # ifelse(isequal(r, x^k), 0, 2) + + b = min(2 * b, f) + end end #= @@ -230,60 +227,61 @@ Write a MIX program that multiplies $(u_1u_2\ldots u_n)_b$ by $v$, where $v$ is LDA U,1 N DEC1 1 N MUL V N J1P 2B N SLC 5 N STX W 1 - + The running time is 23N + K + 5 cycles, and K is rougly (1/2)N. =# """ - mul2(b::Integer, m::Integer, k::Integer) + mul2(b::Integer, m::Integer, k::Integer) The notation mul(r, k) is used in the Bernstein paper to denote k * r. The notation mul_{2, b}(m, k) represents 2-adic multiplication. """ function mul2(b::Integer, m::Integer, k::Integer) - if k ≤ 0 - throw(error("Cannot compute mul(m, k) when k is negative using the 2-adic method.")) - end - - return mod(k*m, 2^b) + if k ≤ 0 + throw(error("Cannot compute mul(m, k) when k is negative using the 2-adic method.")) + end + + return mod(k * m, 2^b) end """ - div2(b::Integer, m::Integer, k::Integer) + div2(b::Integer, m::Integer, k::Integer) + div_b is defined as a floating-point approximation to r/k, so that r/k div_b(r, k) is between 1 and 1 + 2^(1 - b). div2 is the 2-adic approximation of this. """ function div2(b::Integer, m::Integer, k::Integer) - if k ≤ 0 - throw(error("Cannot compute mul(m, k) when k is negative using the 2-adic method.")) - end - - for m in 0:2^b - if isequal(m, mod(div(b, m, k), 2^b)) - return div(b, m, k) - end - end + if k ≤ 0 + throw(error("Cannot compute mul(m, k) when k is negative using the 2-adic method.")) + end + + for m in 0:(2^b) + if isequal(m, mod(div(b, m, k), 2^b)) + return div(b, m, k) + end + end end """ - odd_nroot2(b::Integer, y::Integer, k::Integer) + odd_nroot2(b::Integer, y::Integer, k::Integer) Fix an odd integer y and a positive odd integer k. odd_nroot2 finds an approximate negative kth root of y by Newton's method. """ function odd_nroot2(b::Integer, y::Integer, k::Integer) - if k ≤ 0 || b ≤ 0 - throw(error("Cannot use this method for negative number k and b.")) - end - - if iseven(y) || iseven(k) - throw(error("Cannot use this method for even y and k.")) - end - - b′ = ceil(Integer, b / 2) - - isone(b) && return 1 - z = odd_nroot2(b′, y, k) - r₂ = mul2(b, z, k + 1) - r₃ = mod(y * pow2(b, z, k + 1), 2^b) - return r₄ = div2(b, r₂ - r₃, k) + if k ≤ 0 || b ≤ 0 + throw(error("Cannot use this method for negative number k and b.")) + end + + if iseven(y) || iseven(k) + throw(error("Cannot use this method for even y and k.")) + end + + b′ = ceil(Integer, b / 2) + + isone(b) && return 1 + z = odd_nroot2(b′, y, k) + r₂ = mul2(b, z, k + 1) + r₃ = mod(y * pow2(b, z, k + 1), 2^b) + return r₄ = div2(b, r₂ - r₃, k) end #= @@ -292,28 +290,28 @@ end =# """ - sqrt_nroot2(b::Integer, y::Integer) + sqrt_nroot2(b::Integer, y::Integer) Fix odd integer y. For each b ≥ 1, define and construct nroot_(2, b)(y, 2). """ function sqrt_nroot2(b::Integer, y::Integer) - # if b ≤ 0 - # throw(error("Cannot use this method for negative number b.")) - # end - # - # if iseven(y) - # throw(error("Cannot use this method for even y.")) - # end - - b′ = ceil(Integer, (b + 1) / 2) - - isone(b) && return ifelse(isone(mod(y, 4)), 1, 0) - isequal(b, 2) && return ifelse(isone(mod(y, 8)), 1, 0) - iszero(z) && return 0 - - r₂ = mul2(b + 1, z, 3) - r₃ = mod(y * pow2(b + 1, z, 3), 2^(b + 1)) - return r₄ = mod(r₂ - r₃, 2^b) + # if b ≤ 0 + # throw(error("Cannot use this method for negative number b.")) + # end + # + # if iseven(y) + # throw(error("Cannot use this method for even y.")) + # end + + b′ = ceil(Integer, (b + 1) / 2) + + isone(b) && return ifelse(isone(mod(y, 4)), 1, 0) + isequal(b, 2) && return ifelse(isone(mod(y, 8)), 1, 0) + iszero(z) && return 0 + + r₂ = mul2(b + 1, z, 3) + r₃ = mod(y * pow2(b + 1, z, 3), 2^(b + 1)) + return r₄ = mod(r₂ - r₃, 2^b) end #= @@ -322,31 +320,31 @@ end =# """ - perfect_power_decomp_odd(n::Integer, k::Integer, y::Integer) + perfect_power_decomp_odd(n::Integer, k::Integer, y::Integer) Given a positive odd integer n, an integer k ≥ 2 such that either k = 2 or k is odd, and an odd integer y, `perfect_power_decomp_odd` checks if n is a kth power. """ function perfect_power_decomp_odd(n::Integer, k::Integer, y::Integer) - if iseven(n) && !isequal(k, 2) - throw(error("Cannot use this algorithm on even n.")) - end - - if k < 2 - throw(error("k must be ≥ 2.")) - else - if iseven(k) && !isequal(k, 2) - throw(error("k must be 2 or odd.")) - end - end - - f = floor(Integer, lg(2*n)) - b = ceil(Integer, f / k) - - r = odd_nroot2(b, y, k) - (isequal(k, 2) && iszero(r)) && return 0 - tentative_check_kth_roots(n, r, k) && return r - (isequal(k, 2) && tentative_check_kth_roots(n, 2^b - r, k)) && return 2^b - r - return 0 + if iseven(n) && !isequal(k, 2) + throw(error("Cannot use this algorithm on even n.")) + end + + if k < 2 + throw(error("k must be ≥ 2.")) + else + if iseven(k) && !isequal(k, 2) + throw(error("k must be 2 or odd.")) + end + end + + f = floor(Integer, lg(2 * n)) + b = ceil(Integer, f / k) + + r = odd_nroot2(b, y, k) + (isequal(k, 2) && iszero(r)) && return 0 + tentative_check_kth_roots(n, r, k) && return r + (isequal(k, 2) && tentative_check_kth_roots(n, 2^b - r, k)) && return 2^b - r + return 0 end #= @@ -355,21 +353,21 @@ end =# """ - odd_decomp_perf_pow(n::Integer) + odd_decomp_perf_pow(n::Integer) Given an odd integer n ≥ 2, `odd_decomp_perf_pow` decomposes n as a perfect power if possible. """ function perfectpower(n::Integer) # isprime(n) && return n, 1 - # if iseven(n) && !isequal(n, 2) - # throw(error("This method must have odd n.")) - # end - - if n < 2 - throw(error("n must be ≥ 2.")) - end - + # if iseven(n) && !isequal(n, 2) + # throw(error("This method must have odd n.")) + # end + + if n < 2 + throw(error("n must be ≥ 2.")) + end + isprime(n) && return n, 1 # sqrt_val = sqrt(n) @@ -383,16 +381,16 @@ function perfectpower(n::Integer) end end - f = floor(Integer, lg(2*n)) - - y = odd_nroot2(ceil(Integer, f / 2) + 1, n, 1) + f = floor(Integer, lg(2 * n)) + + y = odd_nroot2(ceil(Integer, f / 2) + 1, n, 1) + + for p in primes(f - 1) + x = perfect_power_decomp_odd(n, p, y) + x > 0 && return x, p + end - for p in primes(f - 1) - x = perfect_power_decomp_odd(n, p, y) - x > 0 && return x, p - end - - return n, 1 + return n, 1 end #= @@ -411,6 +409,6 @@ As usual, fix $n\geq 2$. In this section I discuss several tricks based on comp function lower_exp_bound(n::Integer, x::Integer, k::Integer) # if isodd(n) - # continue + # continue # end end diff --git a/src/wip/primes.jl b/src/wip/primes.jl index b1d237d..2727298 100644 --- a/src/wip/primes.jl +++ b/src/wip/primes.jl @@ -1,19 +1,13 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - # TAKE ONE using Nemo, Primes function _isprime(n) - return Nemo.isprime(Nemo.fmpz(n)) + return Nemo.isprime(Nemo.fmpz(n)) end function _factors(n) - return n == 0 ? [] : Nemo.factor(Nemo.fmpz(n)) + return n == 0 ? [] : Nemo.factor(Nemo.fmpz(n)) end # function _factors(n) # first(i) for i in Primes.factor(n) @@ -23,32 +17,32 @@ function primedivisors(n) n == 0 && return [] _isprime(n) && return [Nemo.fmpz(n)] f = _factors(n) - return sort!([p for (p, e) ∈ f]) + return sort!([p for (p, e) in f]) end function ω(n) - nprimedivisors = 0 - if n == 0 - nprimedivisors = 0 - end - if _isprime(n) - nprimedivisors = length(Nemo.fmpz(n)) - end + nprimedivisors = 0 + if n == 0 + nprimedivisors = 0 + end + if _isprime(n) + nprimedivisors = length(Nemo.fmpz(n)) + end - return length(_factors(n)) + return length(_factors(n)) - # return fmpz(length(primedivisors(n))) + # return fmpz(length(primedivisors(n))) end function Ω(n) n == fmpz(0) && return 0 _isprime(n) && return fmpz(1) - - return sum(e for (_, e) ∈ _factors(n)) + + return sum(e for (_, e) in _factors(n)) end function isperfectpower(n) - return ω(n) == 1 && Ω(n) != 1 + return ω(n) == 1 && Ω(n) != 1 end #--------------------------- @@ -57,47 +51,47 @@ end You find the smallest pair (a, n) (small with respect to n) such that a^n is your target number. Then you need to check if a is prime. =# function isprimepower(n::Integer) - # return any(iszero(mod(n, p)) for p in primes(ceil(Int, sqrt(n)))) - - n ≤ 1 && return false - isprime(n) && return true - ub = isqrt(n) + 1 - - for a in primes(ub) - for m in 2:ub # m > 1 - res = a^m - res > n && break - isequal(res, n) && return true - end - end - - return false + # return any(iszero(mod(n, p)) for p in primes(ceil(Int, sqrt(n)))) + + n ≤ 1 && return false + isprime(n) && return true + ub = isqrt(n) + 1 + + for a in primes(ub) + for m in 2:ub # m > 1 + res = a^m + res > n && break + isequal(res, n) && return true + end + end + + return false end function primepower_brute_force(n::Integer) - n ≤ 1 && return nothing - isprime(n) && return n, 1 - - for a in primes(ceil(Int, sqrt(n))) - for m in 2:ceil(Int, sqrt(n)) # m > 1 - a^m > n && break - isequal(a^m, n) && return a, m - end - end - - return nothing + n ≤ 1 && return nothing + isprime(n) && return n, 1 + + for a in primes(ceil(Int, sqrt(n))) + for m in 2:ceil(Int, sqrt(n)) # m > 1 + a^m > n && break + isequal(a^m, n) && return a, m + end + end + + return nothing end function perfectpower_brute_force(n::Integer) - n ≤ 1 && return nothing - - for a in 1:n - for m in 2:n - isequal(a^m, n) && return m, a - end - end - - return nothing + n ≤ 1 && return nothing + + for a in 1:n + for m in 2:n + isequal(a^m, n) && return m, a + end + end + + return nothing end #========================================================================# @@ -105,8 +99,8 @@ end ## The following function are obtained from DJB struct DJBFloat - a::Signed - n::Unsigned + a::Signed + n::Unsigned end lg(x) = log2(x) @@ -118,34 +112,34 @@ Design an algorithm that divides a nonnegative n-place integer $(u_1u_2\ldots u_ =# function Base.divrem(u::Integer, v::Real, b::Integer) - r, j = 0, 1 - n = ndigits(u) #; base = b - u̲ = round.(reverse(digits(u)), base = b) - w̲ = Vector{Int}(undef, n) - - while j ≤ n - w̲[j] = floor(Int, (r*b + u̲[j]) / v) - r = mod(r*b + u̲[j], v) - j += 1 - end - - # w = sum(w̲[i] * 10 ^ (length(w̲) - i) for i in eachindex(w̲)) - w = parse(Int, join(string.(w̲))) - # w = foldr((i, j) -> 10 * j + i, reverse(w̲)) - # - # n = 0 - # reduce(w̲) do i, j - # 10 * j + i - # # global n = n + ndigits(w̲[i]) - # end - # n = 0 - # w = sum(eachindex(w̲)) do i - # w̲[i] * 10 ^ (length(w̲) - i + n) - # global n = n + ndigits(w̲[i]) - # end - # n = 0 - - return w, r + r, j = 0, 1 + n = ndigits(u) #; base = b + u̲ = round.(reverse(digits(u)), base = b) + w̲ = Vector{Int}(undef, n) + + while j ≤ n + w̲[j] = floor(Int, (r * b + u̲[j]) / v) + r = mod(r * b + u̲[j], v) + j += 1 + end + + # w = sum(w̲[i] * 10 ^ (length(w̲) - i) for i in eachindex(w̲)) + w = parse(Int, join(string.(w̲))) + # w = foldr((i, j) -> 10 * j + i, reverse(w̲)) + # + # n = 0 + # reduce(w̲) do i, j + # 10 * j + i + # # global n = n + ndigits(w̲[i]) + # end + # n = 0 + # w = sum(eachindex(w̲)) do i + # w̲[i] * 10 ^ (length(w̲) - i + n) + # global n = n + ndigits(w̲[i]) + # end + # n = 0 + + return w, r end Base.div(u::Integer, v::Real, b::Integer) = first(divrem(u, v, b)) @@ -161,10 +155,9 @@ In this section, I define truncation to $b$ bits, written $\text{trunc}_b$, and =# function Base.trunc(n::Real, b::Integer) return n & (1 << b - 1) - # return div(n, 1, b) + # return div(n, 1, b) end - #= Write a MIX program that multiplies $(u_1u_2\ldots u_n)_b$ by $v$, where $v$ is a single-precision number (i.e., $0 \leq v < b$), producing the answer $(w_0w_1\ldots w_n)_b$. How much running time is required? @@ -175,13 +168,13 @@ Write a MIX program that multiplies $(u_1u_2\ldots u_n)_b$ by $v$, where $v$ is LDA U,1 N DEC1 1 N MUL V N J1P 2B N SLC 5 N STX W 1 - + The running time is 23N + K + 5 cycles, and K is rougly (1/2)N. =# function mul(u::Integer, v::Real, b::Integer) - n = ndigits(u) #; base = b - u̲ = round.(reverse(digits(u)), base = b) - w̲ = Vector{Int}(undef, n) + n = ndigits(u) #; base = b + u̲ = round.(reverse(digits(u)), base = b) + return w̲ = Vector{Int}(undef, n) end #= @@ -193,13 +186,13 @@ Given a positive floating-point number r and two positive integers b, k, to prim \end{aligned}\end{equation} =# function __pow(r::Real, b::Integer, k::Integer) - if r ≤ 0 || b ≤ 0 || k ≤ 0 - throw(error("r, b, and k must be positive numbers to compute pow_b(r, k).")) - end - println("pow: $r, $b, $k") - isone(k) && return trunc(r, b) - iszero(mod(k, 2)) && return trunc(__pow(r, b, k ÷ 2)^2, b) # use \div because k is even - return trunc(__pow(r, b, k - 1) * trunc(r, b), b) + if r ≤ 0 || b ≤ 0 || k ≤ 0 + throw(error("r, b, and k must be positive numbers to compute pow_b(r, k).")) + end + println("pow: $r, $b, $k") + isone(k) && return trunc(r, b) + iszero(mod(k, 2)) && return trunc(__pow(r, b, k ÷ 2)^2, b) # use \div because k is even + return trunc(__pow(r, b, k - 1) * trunc(r, b), b) end #= @@ -207,113 +200,115 @@ I am trying to find a root $z$ of $z^ky-1$. To compute nroot_b(y, k) for 1 \le b \le \ceil{lg8k}: In advance, find the exponent g satisfying 2^{g- 1} < y < 2^g, and set a = \floor{-g / k}, so that 2^a \le y^{-1/k} < 2^{a+1}. Also set B = \ceil{lg(66(2k + 1))} =# function __nroot(y::Real, k::Real, b::Real) - # find the exponent g satisfying 2^(g - 1) < y < 2^g - g = nothing - i = 0 - while true - if 2^(float(i) - 1) < y && y ≤ float(2)^i - g = i - break - end - i += 1 - end - - a = floor(Int, -g / k) - B = ceil(Int, lg(66 * (2*k + 1))) - z = 2^float(a) + 2^(float(a) - 1) - println("nroot: $z") - j = 1 - - while ! isequal(j, b) - println("nroot: $y, $k, $j, $B") - z = __nroot(y, k, j) - r = trunc(__pow(z, k, B) * trunc(y, B), B) - z = ifelse(r ≤ 993//1024, z + 2^(a - j + 1), z) - z = ifelse(r > 1, z - 2^(a - j - 1), z) - j = j + 1 - end - - return z + # find the exponent g satisfying 2^(g - 1) < y < 2^g + g = nothing + i = 0 + while true + if 2^(float(i) - 1) < y && y ≤ float(2)^i + g = i + break + end + i += 1 + end + + a = floor(Int, -g / k) + B = ceil(Int, lg(66 * (2 * k + 1))) + z = 2^float(a) + 2^(float(a) - 1) + println("nroot: $z") + j = 1 + + while !isequal(j, b) + println("nroot: $y, $k, $j, $B") + z = __nroot(y, k, j) + r = trunc(__pow(z, k, B) * trunc(y, B), B) + z = ifelse(r ≤ 993//1024, z + 2^(a - j + 1), z) + z = ifelse(r > 1, z - 2^(a - j - 1), z) + j = j + 1 + end + + return z end #= To compute nroot_b(y, k) for b \ge 4 + \ceil{lg 8k} + 1: In advance set b' = \ceil{lg 2k} + \ceil{(b - \ceil{lg 2k}) / 2} and B = 2b' + 4 - \ceil{lg k}. Note that b; < b. =# function __newton_nroot(y::Real, k::Real, b::Real) - b′ = ceil(Int, (b - ceil(Int, lg(2*k))) / 2) - B = 2*b′ + 4 - ceil(Int, lg(l)) - z = ifelse(b′ ≤ ceil(Int, lg(8*k)), __nroot(__nroot(y, k, b′)), __newton_nroot(y, k, b′)) # if else, then b′ ≥ ceil(Int, lg(8*k)) + 1 - r₂ = trunc(z, B) * (k + 1) - r₃ = trunc(__pow(z, k + 1, B) * trunc(y, B), B) - - return r₄ = div(r₂ - r₃, k, B) + b′ = ceil(Int, (b - ceil(Int, lg(2 * k))) / 2) + B = 2 * b′ + 4 - ceil(Int, lg(l)) + z = ifelse( + b′ ≤ ceil(Int, lg(8 * k)), __nroot(__nroot(y, k, b′)), __newton_nroot(y, k, b′) + ) # if else, then b′ ≥ ceil(Int, lg(8*k)) + 1 + r₂ = trunc(z, B) * (k + 1) + r₃ = trunc(__pow(z, k + 1, B) * trunc(y, B), B) + + return r₄ = div(r₂ - r₃, k, B) end #= Given positive integers n, x, k, to compute the sign of n - x^k: In advance, set f = \ceil{lg 2n} =# function __sign(n::Integer, x::Integer, k::Integer) - if n ≤ 0 || x ≤ 0 || k ≤ 0 - throw(error("n, x, and k must be positive numbers to compute the sign of n - x^k using this algorithm.")) - end - - f = ceil(Int, lg(2n)) - b = 1 - while true - r = __pow(x, k, b + ceil(Int, lg(8k))) - - n < r && return -1 - r*(1 + 2^(-b)) ≤ n && return 1 - b ≥ f && return 0 - - b = min(2*b, f) - end + if n ≤ 0 || x ≤ 0 || k ≤ 0 + throw(error("n, x, and k must be positive numbers to compute the sign of n - x^k using this algorithm.",),) + end + + f = ceil(Int, lg(2n)) + b = 1 + while true + r = __pow(x, k, b + ceil(Int, lg(8k))) + + n < r && return -1 + r * (1 + 2^(-b)) ≤ n && return 1 + b ≥ f && return 0 + + b = min(2 * b, f) + end end #= Given integers n \ge 2 and k \ge 2, and a positive floating-point number y, to see if n is a kth power: in advance, set f = \floor{lg 2n} and b = 3 + \ceil{f / k} =# function __iskthpower(n::Integer, k::Integer, y::Real) - if n < 2 || k < 2 || y ≤ 0 - thow(error("n, k, and y must be positive (and n and k must be more than one) to see if n is a kth power using this algorithm.")) - end - - f = floor(Int, lg(2*n)) - b = 3 + ceil(Int, f / k) - r = __nroot(y, k, b) - - x = nothing - for i in 1:r - if abs(r - x) ≤ 5//8 - x = i - end - end - - (iszero(x) || abs(r - x) ≥ 1//4) && return 0 - - sign = __sign(n, x, k) # why do we do this step?! - isequal(n, x^k) && return x - - return 0 + if n < 2 || k < 2 || y ≤ 0 + thow(error("n, k, and y must be positive (and n and k must be more than one) to see if n is a kth power using this algorithm.",),) + end + + f = floor(Int, lg(2 * n)) + b = 3 + ceil(Int, f / k) + r = __nroot(y, k, b) + + x = nothing + for i in 1:r + if abs(r - x) ≤ 5//8 + x = i + end + end + + (iszero(x) || abs(r - x) ≥ 1//4) && return 0 + + sign = __sign(n, x, k) # why do we do this step?! + isequal(n, x^k) && return x + + return 0 end #= Given an integer n \ge 2, to decompose n as a perfect power if possible: in advance, set f = \floor{lg 2n} =# function __find_prime_power(n::Integer) - if n < 2 - throw(error("This algorithm to compute the prime power if possible only works for integers more than one.")) - end - - # isequal(n, 2) && return n, 1 - - f = floor(Int, lg(2*n)) - y = __nroot(n, 1, 3 + ceil(Int, f / 2)) - - for p in primes(f) - x = __iskthpower(n, p, y) - x > 0 && return x, p - end - - return x, 1 + if n < 2 + throw(error("This algorithm to compute the prime power if possible only works for integers more than one.",),) + end + + # isequal(n, 2) && return n, 1 + + f = floor(Int, lg(2 * n)) + y = __nroot(n, 1, 3 + ceil(Int, f / 2)) + + for p in primes(f) + x = __iskthpower(n, p, y) + x > 0 && return x, p + end + + return x, 1 end diff --git a/test/runtests.jl b/test/runtests.jl index c59c106..dd4c767 100755 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,162 +1,362 @@ -#!/usr/bin/env bash - #= - exec julia --project="$(realpath $(dirname $(realpath $(dirname $0))))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \ - "${BASH_SOURCE[0]}" "$@" - =# - -include(joinpath(dirname(dirname(@__FILE__)), "src", "CodingTheory.jl")) - -using .CodingTheory -import .CodingTheory.deepsym, .CodingTheory.displaymatrix +using CodingTheory +import CodingTheory.deepsym, CodingTheory.displaymatrix using Test using Polynomials - function aresamearrays(A::AbstractVector{T}, B::AbstractVector{R}) where {T, R} - A === B && return true - length(A) ≠ length(B) && return false - - for (a, b) in zip(sort(A), sort(B)) - if a ≠ b - return false - end - end - - return true + A === B && return true + length(A) ≠ length(B) && return false + + for (a, b) in zip(sort(A), sort(B)) + if a ≠ b + return false + end + end + + return true end """ Given two arrays A and B, ensures they have the same elements """ function aresamearrays(A::AbstractArray{T, N}, B::AbstractArray{R, M}) where {T, R, N, M} - A === B && return true - length(A) ≠ length(B) && return false - - for (a, b) in zip(sort(reshape(A, :)), sort(reshape(B, :))) - # for (a, b) in zip(sort(A, dims = 1), sort(B, dims = 1)) - if a ≠ b - return false - end - end - - return true + A === B && return true + length(A) ≠ length(B) && return false + + for (a, b) in zip(sort(reshape(A, :)), sort(reshape(B, :))) + # for (a, b) in zip(sort(A, dims = 1), sort(B, dims = 1)) + if a ≠ b + return false + end + end + + return true end @time @testset "CodingTheory.jl" begin - @test has_identity([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0], 1) == true - @test has_identity([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true - @test has_identity([1 0 1 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true - @test has_identity([1 0 1 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0], 3) == false - @test has_identity([-96 -66 20 1 0 0; -65 59 -82 0 1 0; -16 87 -113 0 0 1]) == true - @test has_identity([78 -99 125 -123 -111 -71 17; -115 78 40 -88 81 -40 78; -99 126 -54 1 0 0 24; -55 88 42 0 1 0 -8; 119 55 2 0 0 1 -92; -40 -21 -89 -79 59 -44 9]) == true - @test has_identity([78 -99 125 -123 -111 -71 17; -115 78 40 -88 81 -40 78; -99 126 -54 1 0 0 24; -55 88 42 0 1 0 -8; 119 55 2 0 0 1 -92; -40 -21 -89 -79 59 -44 9], 4) == false - @test has_identity_on_left([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true - @test has_identity_on_left([-96 -66 20 1 0 0; -65 59 -82 0 1 0; -16 87 -113 0 0 1]) == false - - @test hamming_distance("ABC", "DBC") == 1 - @test hamming_distance("ABC", "DEF") == 3 - - @test aresamearrays(hamming_ball([[1, 0, 1], [0, 1, 1], [1, 0, 0]], [1, 0, 0], 2), deepsym([[1, 0, 1], [1, 0, 0]])) - @test aresamearrays(hamming_ball(list_span([2, 2, 2], [1, 2, 0], [0, 2, 1], 3), [1, 0, 0], 3), deepsym([[0, 0, 0], [0, 2, 1], [0, 1, 2], [1, 2, 0], [1, 1, 1], [1, 0, 2], [2, 1, 0], [2, 0, 1], [2, 2, 2]])) - - @test rate(3, 5, 4) ≈ 0.3662433802 - - @test t_error_correcting([[0, 0, 0, 0], [2, 2, 2, 2]], 1) == true - @test t_error_correcting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false - @test t_error_detecting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false - - @test find_error_detection_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) == 1 - @test find_error_correction_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) == 0 - @test find_error_correction_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 0 - @test find_error_detection_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 1 - - @test Polynomial([1, 2, 3, 4, 5, 6, 7, 8, 9], 3) == Polynomial([1, 2, 0, 1, 2, 0, 1, 2]) - - @test isirreducible(Polynomial([1, 1, 0, 0, 1]), 2) == true - @test isirreducible(Polynomial([1, 1, 1, 0, 1, 1]), 2) == true - @test isirreducible(Polynomial([4, 4, 1]), 2) == false - @test isirreducible(Polynomial([1, 0, 1]), 2) == false - @test isirreducible(Polynomial([-2, 0, 1]), 2) == false - @test isirreducible(Polynomial([1, 1]), 2) == false - - p = Polynomial([1, 1, 2, 0, 1, 2, 1]) - q = Polynomial([2, 1, 1]) - a = Polynomial([0, 1, 0, 1, 0, 0, 1]) - b = Polynomial([1, 0, 1]) - - @test mod(rem(p, q), 3) == Polynomial([1, 1]) - @test mod(rem(a, b), 2) == Polynomial([1]) - + @test has_identity([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0], 1) == true + @test has_identity([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true + @test has_identity([1 0 1 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true + @test has_identity([1 0 1 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0], 3) == false + @test has_identity([-96 -66 20 1 0 0; -65 59 -82 0 1 0; -16 87 -113 0 0 1]) == true + @test has_identity([ + 78 -99 125 -123 -111 -71 17 + -115 78 40 -88 81 -40 78 + -99 126 -54 1 0 0 24 + -55 88 42 0 1 0 -8 + 119 55 2 0 0 1 -92 + -40 -21 -89 -79 59 -44 9 + ],) == true + @test has_identity( + [ + 78 -99 125 -123 -111 -71 17 + -115 78 40 -88 81 -40 78 + -99 126 -54 1 0 0 24 + -55 88 42 0 1 0 -8 + 119 55 2 0 0 1 -92 + -40 -21 -89 -79 59 -44 9 + ], + 4, + ) == false + @test has_identity_on_left([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true + @test has_identity_on_left([-96 -66 20 1 0 0; -65 59 -82 0 1 0; -16 87 -113 0 0 1]) == + false + + @test hamming_distance("ABC", "DBC") == 1 + @test hamming_distance("ABC", "DEF") == 3 + + @test aresamearrays( + hamming_ball([[1, 0, 1], [0, 1, 1], [1, 0, 0]], [1, 0, 0], 2), + deepsym([[1, 0, 1], [1, 0, 0]]), + ) + @test aresamearrays( + hamming_ball(list_span([2, 2, 2], [1, 2, 0], [0, 2, 1], 3), [1, 0, 0], 3), + deepsym([ + [0, 0, 0], + [0, 2, 1], + [0, 1, 2], + [1, 2, 0], + [1, 1, 1], + [1, 0, 2], + [2, 1, 0], + [2, 0, 1], + [2, 2, 2], + ]), + ) + + @test rate(3, 5, 4) ≈ 0.3662433802 + + @test t_error_correcting([[0, 0, 0, 0], [2, 2, 2, 2]], 1) == true + @test t_error_correcting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false + @test t_error_detecting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false + + @test find_error_detection_max( + [[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2 + ) == 1 + @test find_error_correction_max( + [[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2 + ) == 0 + @test find_error_correction_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 0 + @test find_error_detection_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 1 + + @test Polynomial([1, 2, 3, 4, 5, 6, 7, 8, 9], 3) == Polynomial([1, 2, 0, 1, 2, 0, 1, 2]) + + @test isirreducible(Polynomial([1, 1, 0, 0, 1]), 2) == true + @test isirreducible(Polynomial([1, 1, 1, 0, 1, 1]), 2) == true + @test isirreducible(Polynomial([4, 4, 1]), 2) == false + @test isirreducible(Polynomial([1, 0, 1]), 2) == false + @test isirreducible(Polynomial([-2, 0, 1]), 2) == false + @test isirreducible(Polynomial([1, 1]), 2) == false + + p = Polynomial([1, 1, 2, 0, 1, 2, 1]) + q = Polynomial([2, 1, 1]) + a = Polynomial([0, 1, 0, 1, 0, 0, 1]) + b = Polynomial([1, 0, 1]) + + @test mod(rem(p, q), 3) == Polynomial([1, 1]) + @test mod(rem(a, b), 2) == Polynomial([1]) + @test isperfectpower(36) == true - @test isprimepower(36) == false - @test isperfectpower(9) == true - @test isprimepower(9) == true - @test isperfectpower(5) == false - @test isprimepower(2) == true - @test isprimepower(18) == false - - @test rref([1 0 1 1 1; 1 1 1 0 1; 0 1 1 1 1], 2) == [1 0 0 1 0; 0 1 0 1 0; 0 0 1 0 1] - @test rref([1 0 1 0 1 0; 0 1 0 0 1 0; 1 1 1 1 1 1], 2) == [1 0 1 0 1 0; 0 1 0 0 1 0; 0 0 0 1 1 1] - @test rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=false) == [1 0 3 0 2 2; 0 1 2 0 1 1; 0 0 0 1 0 4] - @test rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=true) == [1 0 0 3 2 2; 0 1 0 2 1 1; 0 0 1 0 0 4] - @test rref([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] - @test rref([0 0 0 0 0; 1 0 1 0 1; 0 1 0 1 1; 1 1 1 1 0], 2) == [1 0 1 0 1; 0 1 0 1 1; 0 0 0 0 0; 0 0 0 0 0] - @test rref([1 1 1 0; 1 1 0 1; 0 0 1 1], 2) == [1 1 0 1; 0 0 1 1; 0 0 0 0] - - @test multiplication_table(2, 3) == Polynomial[Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]); Polynomial([0]) Polynomial([1]) Polynomial([2]) Polynomial([0, 1]) Polynomial([1, 1]) Polynomial([2, 1]) Polynomial([0, 2]) Polynomial([1, 2]) Polynomial([2, 2]); Polynomial([0]) Polynomial([2]) Polynomial([1]) Polynomial([0, 2]) Polynomial([2, 2]) Polynomial([1, 2]) Polynomial([0, 1]) Polynomial([2, 1]) Polynomial([1, 1]); Polynomial([0]) Polynomial([0, 1]) Polynomial([0, 2]) Polynomial([0, 0, 1]) Polynomial([0, 1, 1]) Polynomial([0, 2, 1]) Polynomial([0, 0, 2]) Polynomial([0, 1, 2]) Polynomial([0, 2, 2]); Polynomial([0]) Polynomial([1, 1]) Polynomial([2, 2]) Polynomial([0, 1, 1]) Polynomial([1, 2, 1]) Polynomial([2, 0, 1]) Polynomial([0, 2, 2]) Polynomial([1, 0, 2]) Polynomial([2, 1, 2]); Polynomial([0]) Polynomial([2, 1]) Polynomial([1, 2]) Polynomial([0, 2, 1]) Polynomial([2, 0, 1]) Polynomial([1, 1, 1]) Polynomial([0, 1, 2]) Polynomial([2, 2, 2]) Polynomial([1, 0, 2]); Polynomial([0]) Polynomial([0, 2]) Polynomial([0, 1]) Polynomial([0, 0, 2]) Polynomial([0, 2, 2]) Polynomial([0, 1, 2]) Polynomial([0, 0, 1]) Polynomial([0, 2, 1]) Polynomial([0, 1, 1]); Polynomial([0]) Polynomial([1, 2]) Polynomial([2, 1]) Polynomial([0, 1, 2]) Polynomial([1, 0, 2]) Polynomial([2, 2, 2]) Polynomial([0, 2, 1]) Polynomial([1, 1, 1]) Polynomial([2, 0, 1]); Polynomial([0]) Polynomial([2, 2]) Polynomial([1, 1]) Polynomial([0, 2, 2]) Polynomial([2, 1, 2]) Polynomial([1, 0, 2]) Polynomial([0, 1, 1]) Polynomial([2, 0, 1]) Polynomial([1, 2, 1])] - - @test aresamearrays(list_span([2, 1, 1], [1, 1, 1], 3), [[0, 0, 0], [1, 1, 1], [2, 2, 2], [2, 1, 1], [0, 2, 2], [1, 0, 0], [1, 2, 2], [2, 0, 0], [0, 1, 1]]) - - @test islinear([[0,0,0],[1,1,1],[1,0,1],[1,1,0]], 2) == false - @test islinear([[0,0,0],[1,1,1],[1,0,1],[0,1,0]], 2) == true - - @test code_distance([[0,0,0,0,0],[1,0,1,0,1],[0,1,0,1,0],[1,1,1,1,1]]) == 2 - @test code_distance([[0,0,0,0,0],[1,1,1,0,0],[0,0,0,1,1],[1,1,1,1,1],[1,0,0,1,1],[0,1,1,0,0]]) == 1 - - w1 = Word((1, 5, 2, 6)); w2 = Word([4, 2, 6, 7, 3, 10, 10]); w3 = Word(1, 3, 6, 2, 4) - @test all(isword.([w1, w2, w3])) == true - @test all(isabstractword.([w1, w2, w3])) == true - @test isabstractword(1, 2, 3) == true - @test all(isabstractword.([(1, 2, 3), [1, 2, 3]])) == true - @test any(isword.([(1, 2, 3), [4, 5, 6]])) == false - @test isword(1, 2, 3) == false - @test aresamearrays(Alphabet("123"), deepsym([1, 2, 3])) - @test aresamearrays(Alphabet([1, 2, 3]), deepsym([1, 2, 3])) - @test aresamearrays(Alphabet(["1", "2", "3"]), deepsym([1, 2, 3])) - # TODO: write test for CodeUniverse struct - @test aresamearrays([i for i in CodeUniverseIterator(["a", "b", "c"], 3)], get_all_words(["a", "b", "c"], 3)) # implicitly tests CodeUniverseIterator - @test aresamearrays(collect(CodeUniverseIterator(["a", "b", "c"], 4)), Tuple[(:a, :a, :a, :a), (:b, :a, :a, :a), (:c, :a, :a, :a), (:a, :b, :a, :a), (:b, :b, :a, :a), (:c, :b, :a, :a), (:a, :c, :a, :a), (:b, :c, :a, :a), (:c, :c, :a, :a), (:a, :a, :b, :a), (:b, :a, :b, :a), (:c, :a, :b, :a), (:a, :b, :b, :a), (:b, :b, :b, :a), (:c, :b, :b, :a), (:a, :c, :b, :a), (:b, :c, :b, :a), (:c, :c, :b, :a), (:a, :a, :c, :a), (:b, :a, :c, :a), (:c, :a, :c, :a), (:a, :b, :c, :a), (:b, :b, :c, :a), (:c, :b, :c, :a), (:a, :c, :c, :a), (:b, :c, :c, :a), (:c, :c, :c, :a), (:a, :a, :a, :b), (:b, :a, :a, :b), (:c, :a, :a, :b), (:a, :b, :a, :b), (:b, :b, :a, :b), (:c, :b, :a, :b), (:a, :c, :a, :b), (:b, :c, :a, :b), (:c, :c, :a, :b), (:a, :a, :b, :b), (:b, :a, :b, :b), (:c, :a, :b, :b), (:a, :b, :b, :b), (:b, :b, :b, :b), (:c, :b, :b, :b), (:a, :c, :b, :b), (:b, :c, :b, :b), (:c, :c, :b, :b), (:a, :a, :c, :b), (:b, :a, :c, :b), (:c, :a, :c, :b), (:a, :b, :c, :b), (:b, :b, :c, :b), (:c, :b, :c, :b), (:a, :c, :c, :b), (:b, :c, :c, :b), (:c, :c, :c, :b), (:a, :a, :a, :c), (:b, :a, :a, :c), (:c, :a, :a, :c), (:a, :b, :a, :c), (:b, :b, :a, :c), (:c, :b, :a, :c), (:a, :c, :a, :c), (:b, :c, :a, :c), (:c, :c, :a, :c), (:a, :a, :b, :c), (:b, :a, :b, :c), (:c, :a, :b, :c), (:a, :b, :b, :c), (:b, :b, :b, :c), (:c, :b, :b, :c), (:a, :c, :b, :c), (:b, :c, :b, :c), (:c, :c, :b, :c), (:a, :a, :c, :c), (:b, :a, :c, :c), (:c, :a, :c, :c), (:a, :b, :c, :c), (:b, :b, :c, :c), (:c, :b, :c, :c), (:a, :c, :c, :c), (:b, :c, :c, :c), (:c, :c, :c, :c)]) - - @test sphere_covering_bound(5,7,3) == 215 - @test sphere_packing_bound(5,7,3) == 2693 - - @test aresamearrays(construct_ham_matrix(3,2), [0 0 0 1 1 1 1; 0 1 1 0 0 1 1; 1 0 1 0 1 0 1]) - @test aresamearrays(construct_ham_matrix(3,3), [0 0 0 0 0 0 0 0 1 1 1 1 1; 0 0 1 1 1 2 2 2 0 0 0 1 1; 1 2 0 1 2 0 1 2 0 1 2 0 1]) - - @test isperfect(11, 6, 5, 3) == true - @test isperfect(23, 12, 7, 2) == true - @test isperfect(23, 12, 7, 3) == false - @test isperfect(11, 6, 5, 4) == false - @test isgolayperfect(11, 6, 5, 3) == true - @test isgolayperfect(23, 12, 7, 2) == true - @test isgolayperfect(23, 12, 7, 3) == false - @test isgolayperfect(11, 6, 5, 4) == false - - @test length(get_codewords(5, 5, 3)) ∈ [74:74...] # implicitly tests Word struct, and `get_codewords_random` - @test length(get_codewords(4, 7, 3; m = 1)) ∈ [256:308...] - @test length(get_codewords_greedy(5, 5, 3)) == 74 - randq, randn = rand(1:8, 2) - @test length(get_all_words(randq, randn)) == big(randq)^randn - @test aresamearrays(get_codewords([1 0 1 0; 0 1 1 1], 2), [[0, 0, 0, 0], [1, 0, 1, 0], [0, 1, 1, 1], [1, 1, 0, 1]]) - @test aresamearrays(get_codewords([1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1], 2), [[0, 0, 0, 0, 0, 0], [1, 0, 0, 1, 1, 0], [0, 1, 0, 1, 0, 1], [1, 1, 0, 0, 1, 1], [0, 0, 1, 0, 1, 1], [1, 0, 1, 1, 0, 1], [0, 1, 1, 1, 1, 0], [1, 1, 1, 0, 0, 0]]) - - # note that the following tests are not to be checked using `aresamearrays` because their order matters - @test syndrome([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == [0 0 0] - @test parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] - @test normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] - @test equivalent_code([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] - @test parity_check(normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3), 3) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] - @test isincode([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == true - @test isincode([1, 0, 2, 2, 1, 2, 1], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == false + @test isprimepower(36) == false + @test isperfectpower(9) == true + @test isprimepower(9) == true + @test isperfectpower(5) == false + @test isprimepower(2) == true + @test isprimepower(18) == false + + @test rref([1 0 1 1 1; 1 1 1 0 1; 0 1 1 1 1], 2) == [1 0 0 1 0; 0 1 0 1 0; 0 0 1 0 1] + @test rref([1 0 1 0 1 0; 0 1 0 0 1 0; 1 1 1 1 1 1], 2) == + [1 0 1 0 1 0; 0 1 0 0 1 0; 0 0 0 1 1 1] + @test rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5; colswap = false) == + [1 0 3 0 2 2; 0 1 2 0 1 1; 0 0 0 1 0 4] + @test rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5; colswap = true) == + [1 0 0 3 2 2; 0 1 0 2 1 1; 0 0 1 0 0 4] + @test rref([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] + @test rref([0 0 0 0 0; 1 0 1 0 1; 0 1 0 1 1; 1 1 1 1 0], 2) == + [1 0 1 0 1; 0 1 0 1 1; 0 0 0 0 0; 0 0 0 0 0] + @test rref([1 1 1 0; 1 1 0 1; 0 0 1 1], 2) == [1 1 0 1; 0 0 1 1; 0 0 0 0] + + @test multiplication_table(2, 3) == Polynomial[ + Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) + Polynomial([0]) Polynomial([1]) Polynomial([2]) Polynomial([0, 1]) Polynomial([1, 1]) Polynomial([2, 1]) Polynomial([0, 2]) Polynomial([1, 2]) Polynomial([2, 2]) + Polynomial([0]) Polynomial([2]) Polynomial([1]) Polynomial([0, 2]) Polynomial([2, 2]) Polynomial([1, 2]) Polynomial([0, 1]) Polynomial([2, 1]) Polynomial([1, 1]) + Polynomial([0]) Polynomial([0, 1]) Polynomial([0, 2]) Polynomial([0, 0, 1]) Polynomial([0, 1, 1]) Polynomial([0, 2, 1]) Polynomial([0, 0, 2]) Polynomial([0, 1, 2]) Polynomial([0, 2, 2]) + Polynomial([0]) Polynomial([1, 1]) Polynomial([2, 2]) Polynomial([0, 1, 1]) Polynomial([1, 2, 1]) Polynomial([2, 0, 1]) Polynomial([0, 2, 2]) Polynomial([1, 0, 2]) Polynomial([2, 1, 2]) + Polynomial([0]) Polynomial([2, 1]) Polynomial([1, 2]) Polynomial([0, 2, 1]) Polynomial([2, 0, 1]) Polynomial([1, 1, 1]) Polynomial([0, 1, 2]) Polynomial([2, 2, 2]) Polynomial([1, 0, 2]) + Polynomial([0]) Polynomial([0, 2]) Polynomial([0, 1]) Polynomial([0, 0, 2]) Polynomial([0, 2, 2]) Polynomial([0, 1, 2]) Polynomial([0, 0, 1]) Polynomial([0, 2, 1]) Polynomial([0, 1, 1]) + Polynomial([0]) Polynomial([1, 2]) Polynomial([2, 1]) Polynomial([0, 1, 2]) Polynomial([1, 0, 2]) Polynomial([2, 2, 2]) Polynomial([0, 2, 1]) Polynomial([1, 1, 1]) Polynomial([2, 0, 1]) + Polynomial([0]) Polynomial([2, 2]) Polynomial([1, 1]) Polynomial([0, 2, 2]) Polynomial([2, 1, 2]) Polynomial([1, 0, 2]) Polynomial([0, 1, 1]) Polynomial([2, 0, 1]) Polynomial([1, 2, 1]) + ] + + @test aresamearrays( + list_span([2, 1, 1], [1, 1, 1], 3), + [ + [0, 0, 0], + [1, 1, 1], + [2, 2, 2], + [2, 1, 1], + [0, 2, 2], + [1, 0, 0], + [1, 2, 2], + [2, 0, 0], + [0, 1, 1], + ], + ) + + @test islinear([[0, 0, 0], [1, 1, 1], [1, 0, 1], [1, 1, 0]], 2) == false + @test islinear([[0, 0, 0], [1, 1, 1], [1, 0, 1], [0, 1, 0]], 2) == true + + @test code_distance([ + [0, 0, 0, 0, 0], [1, 0, 1, 0, 1], [0, 1, 0, 1, 0], [1, 1, 1, 1, 1] + ]) == 2 + @test code_distance([ + [0, 0, 0, 0, 0], + [1, 1, 1, 0, 0], + [0, 0, 0, 1, 1], + [1, 1, 1, 1, 1], + [1, 0, 0, 1, 1], + [0, 1, 1, 0, 0], + ]) == 1 + + w1 = Word((1, 5, 2, 6)) + w2 = Word([4, 2, 6, 7, 3, 10, 10]) + w3 = Word(1, 3, 6, 2, 4) + @test all(isword.([w1, w2, w3])) == true + @test all(isabstractword.([w1, w2, w3])) == true + @test isabstractword(1, 2, 3) == true + @test all(isabstractword.([(1, 2, 3), [1, 2, 3]])) == true + @test any(isword.([(1, 2, 3), [4, 5, 6]])) == false + @test isword(1, 2, 3) == false + @test aresamearrays(Alphabet("123"), deepsym([1, 2, 3])) + @test aresamearrays(Alphabet([1, 2, 3]), deepsym([1, 2, 3])) + @test aresamearrays(Alphabet(["1", "2", "3"]), deepsym([1, 2, 3])) + # TODO: write test for CodeUniverse struct + @test aresamearrays( + [i for i in CodeUniverseIterator(["a", "b", "c"], 3)], + get_all_words(["a", "b", "c"], 3), + ) # implicitly tests CodeUniverseIterator + @test aresamearrays( + collect(CodeUniverseIterator(["a", "b", "c"], 4)), + Tuple[ + (:a, :a, :a, :a), + (:b, :a, :a, :a), + (:c, :a, :a, :a), + (:a, :b, :a, :a), + (:b, :b, :a, :a), + (:c, :b, :a, :a), + (:a, :c, :a, :a), + (:b, :c, :a, :a), + (:c, :c, :a, :a), + (:a, :a, :b, :a), + (:b, :a, :b, :a), + (:c, :a, :b, :a), + (:a, :b, :b, :a), + (:b, :b, :b, :a), + (:c, :b, :b, :a), + (:a, :c, :b, :a), + (:b, :c, :b, :a), + (:c, :c, :b, :a), + (:a, :a, :c, :a), + (:b, :a, :c, :a), + (:c, :a, :c, :a), + (:a, :b, :c, :a), + (:b, :b, :c, :a), + (:c, :b, :c, :a), + (:a, :c, :c, :a), + (:b, :c, :c, :a), + (:c, :c, :c, :a), + (:a, :a, :a, :b), + (:b, :a, :a, :b), + (:c, :a, :a, :b), + (:a, :b, :a, :b), + (:b, :b, :a, :b), + (:c, :b, :a, :b), + (:a, :c, :a, :b), + (:b, :c, :a, :b), + (:c, :c, :a, :b), + (:a, :a, :b, :b), + (:b, :a, :b, :b), + (:c, :a, :b, :b), + (:a, :b, :b, :b), + (:b, :b, :b, :b), + (:c, :b, :b, :b), + (:a, :c, :b, :b), + (:b, :c, :b, :b), + (:c, :c, :b, :b), + (:a, :a, :c, :b), + (:b, :a, :c, :b), + (:c, :a, :c, :b), + (:a, :b, :c, :b), + (:b, :b, :c, :b), + (:c, :b, :c, :b), + (:a, :c, :c, :b), + (:b, :c, :c, :b), + (:c, :c, :c, :b), + (:a, :a, :a, :c), + (:b, :a, :a, :c), + (:c, :a, :a, :c), + (:a, :b, :a, :c), + (:b, :b, :a, :c), + (:c, :b, :a, :c), + (:a, :c, :a, :c), + (:b, :c, :a, :c), + (:c, :c, :a, :c), + (:a, :a, :b, :c), + (:b, :a, :b, :c), + (:c, :a, :b, :c), + (:a, :b, :b, :c), + (:b, :b, :b, :c), + (:c, :b, :b, :c), + (:a, :c, :b, :c), + (:b, :c, :b, :c), + (:c, :c, :b, :c), + (:a, :a, :c, :c), + (:b, :a, :c, :c), + (:c, :a, :c, :c), + (:a, :b, :c, :c), + (:b, :b, :c, :c), + (:c, :b, :c, :c), + (:a, :c, :c, :c), + (:b, :c, :c, :c), + (:c, :c, :c, :c), + ], + ) + + @test sphere_covering_bound(5, 7, 3) == 215 + @test sphere_packing_bound(5, 7, 3) == 2693 + + @test aresamearrays( + construct_ham_matrix(3, 2), [0 0 0 1 1 1 1; 0 1 1 0 0 1 1; 1 0 1 0 1 0 1] + ) + @test aresamearrays( + construct_ham_matrix(3, 3), + [0 0 0 0 0 0 0 0 1 1 1 1 1; 0 0 1 1 1 2 2 2 0 0 0 1 1; 1 2 0 1 2 0 1 2 0 1 2 0 1], + ) + + @test isperfect(11, 6, 5, 3) == true + @test isperfect(23, 12, 7, 2) == true + @test isperfect(23, 12, 7, 3) == false + @test isperfect(11, 6, 5, 4) == false + @test isgolayperfect(11, 6, 5, 3) == true + @test isgolayperfect(23, 12, 7, 2) == true + @test isgolayperfect(23, 12, 7, 3) == false + @test isgolayperfect(11, 6, 5, 4) == false + + @test length(get_codewords(5, 5, 3)) ∈ [74:74...] # implicitly tests Word struct, and `get_codewords_random` + @test length(get_codewords(4, 7, 3; m = 1)) ∈ [256:308...] + @test length(get_codewords_greedy(5, 5, 3)) == 74 + randq, randn = rand(1:8, 2) + @test length(get_all_words(randq, randn)) == big(randq)^randn + @test aresamearrays( + get_codewords([1 0 1 0; 0 1 1 1], 2), + [[0, 0, 0, 0], [1, 0, 1, 0], [0, 1, 1, 1], [1, 1, 0, 1]], + ) + @test aresamearrays( + get_codewords([1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1], 2), + [ + [0, 0, 0, 0, 0, 0], + [1, 0, 0, 1, 1, 0], + [0, 1, 0, 1, 0, 1], + [1, 1, 0, 0, 1, 1], + [0, 0, 1, 0, 1, 1], + [1, 0, 1, 1, 0, 1], + [0, 1, 1, 1, 1, 0], + [1, 1, 1, 0, 0, 0], + ], + ) + + # note that the following tests are not to be checked using `aresamearrays` because their order matters + @test syndrome( + [0, 2, 1, 2, 0, 1, 0], + transpose(parity_check( + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3 + )), + 3, + ) == [0 0 0] + @test parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3) == + [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] + @test normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] + @test equivalent_code( + [1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3 + ) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1] + @test parity_check( + normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3), 3 + ) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1] + @test isincode( + [0, 2, 1, 2, 0, 1, 0], + transpose(parity_check( + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3 + )), + 3, + ) == true + @test isincode( + [1, 0, 2, 2, 1, 2, 1], + transpose(parity_check( + [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3 + )), + 3, + ) == false end # end runtests