Created
October 27, 2022 11:49
-
-
Save martian0x80/48bb5c8799bdfea0c986414f812a8c59 to your computer and use it in GitHub Desktop.
Partition function, coin-change problem, generating integer partitions in (anti-)/lexicographic order.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using DataStructures | |
using DSP | |
#= | |
ZS1 and ZS2 based on ALGORITHMS explained in : | |
FAST ALGORITHMS FOR GENERATING | |
INTEGER PARTITIONS | |
ANTOINE ZOGHBIU and IVAN STOJMENOVIC* | |
Bell Northern Research, P.O . Box 35 I I , Station C , Mail Stop 091 , | |
Ottawa, Ontario KIY 4H7, Canada; | |
Computer Science Department, SITE, University of Ottawa, | |
Ottawa. Ontario KIN 6N5. Canada | |
=# | |
function ZS1(n::Int64) | |
partitions = Array{Array{Int64, 1}, 1}() | |
x_i = Dict{Int64, Int64}([i=>1 for i in 1:n]) | |
x_i[1] = n | |
m = 1 | |
h = 1 | |
push!(partitions, Vector{Int64}([x_i[1]])) | |
while x_i[1] ≠ 1 | |
if x_i[h] == 2 | |
m += 1 | |
x_i[h] = 1 | |
h = h - 1 | |
else | |
r = x_i[h] - 1 | |
t = m - h + 1 | |
x_i[h] = r | |
while t ≥ r | |
h = h + 1 | |
x_i[h] = r | |
t = t - r | |
end | |
if t == 0 | |
m = h | |
else | |
m = h + 1 | |
if t > 1 | |
h = h + 1 | |
x_i[h] = t | |
end | |
end | |
end | |
push!(partitions, Array{Int64, 1}([x_i[j] for j∈1:m])) | |
end | |
return partitions | |
end | |
function ZS2(n::Int64) | |
partitions = Array{Array{Int64, 1}, 1}() | |
x_i = Dict{Int64, Int64}([i=>1 for i in 1:n]) | |
push!(partitions, [1 for j in 1:n]) | |
x_i[0] = -1 | |
x_i[1] = 2 | |
h = 1 | |
m = n - 1 | |
push!(partitions, [x_i[j] for j in 1:m]) | |
while x_i[1] ≠ n | |
if m - h > 1 | |
h = h + 1 | |
x_i[h] = 2 | |
m = m - 1 | |
else | |
j = m - 2 | |
while x_i[j] == x_i[m-1] | |
x_i[j] = 1 | |
j = j - 1 | |
end | |
h = j + 1 | |
x_i[h] = x_i[m-1] + 1 | |
r = x_i[m] + x_i[m-1]*(m-h-1) | |
x_i[m] = 1 | |
if m - h > 1 | |
x_i[m-1] = 1 | |
end | |
m = h + r - 1 | |
end | |
push!(partitions, Array{Int64, 1}([x_i[j] for j∈1:m])) | |
end | |
return partitions | |
end | |
function partitions_recurrence(N::Int64, k::Int64, dict::Bool=false, ordered::Bool=false) | |
#= Generates p(n) (unrestricted partition function) ∀ n ∈ 1 to N, returns p(k) | |
Based on the recurrence relation: | |
p(n) = p(n-1) + p(n-2) + p(n-5) + p(n-7) + ... | |
Efficient for values under 10⁵. | |
=# | |
generalized_pentnums = k -> Tuple{Int64, Int64}([((3k^2-k)//2).num, ((3k^2+k)//2).num]) | |
if ordered | |
cached_p = OrderedDict{Int64, BigInt}([0 => 1, 1 => 1, 2 => 2]) | |
else | |
cached_p = Dict{Int64, BigInt}([0 => 1, 1 => 1, 2 => 2]) | |
end | |
for j in 3:N | |
p_j::BigInt = 0 | |
i = 1 | |
Λ = Tuple{Function, Function}((+, -)) | |
while true | |
exponents = generalized_pentnums(i) | |
# checks to see if p(N-x) is valid (i.e., N-x>0), where x is a generalized pentagonal number as described in the pentagonal number theorem. | |
if exponents[1] > j | |
# println("Breaking, $j-$(exponents[1]) < 0") | |
break | |
else | |
# println("p($j) $(Λ[1])= p($(j-exponents[1]))") | |
p_j = Λ[1](p_j, cached_p[j-exponents[1]]) | |
end | |
if exponents[2] > j | |
# println("Breaking, $j-$(exponents[2]) < 0") | |
break | |
else | |
# println("p($j) $(Λ[1])= p($(j-exponents[2]))") | |
p_j = Λ[1](p_j, cached_p[j-exponents[2]]) | |
end | |
Λ = reverse(Λ) | |
i += 1 | |
end | |
cached_p[j] = p_j | |
end | |
dict && return cached_p; | |
return cached_p[k]; | |
end | |
#= | |
The partitions of $n$ with denominations $S$ is encoded by the generating function $$\dfrac{1}{(q;q)_{s\in S}}=\prod_{s\in S}\dfrac{1}{\left(1-q^s\right)}=\sum_{s\geq 0}p(s)q^s$$ | |
I wrote a function that converts the generating function of the form $\dfrac{1}{1-q^i}$ to it's formal power series representation, i.e., $1+q^i+q^{2i}+q^{3i}+\cdots$. | |
Then computed the generalized euler function or [url=https://en.wikipedia.org/wiki/Q-Pochhammer_symbol]q-shifted factorial[/url] by performing the [url=https://en.wikipedia.org/wiki/Cauchy_product]Cauchy Product[/url] (Discrete Convolution) of all the series represented by arrays upto a certain degree ($N = 200$). | |
=# | |
function fn_series(n::Int64, length::Int64) :: Vector{Int8} | |
#= | |
Converts a generating function of the form 1/(1-x^n) to the formal power series it represents. | |
1/(1-x^i) = 1 + x^i + x^2i + x^3i + ... (length) | |
=# | |
series = zeros(Int8, length + 1) | |
series[1] = Int8(1) | |
i = 1 | |
while i < (length + 1) ÷ n | |
series[i*n + 1] = Int8(1) | |
i += 1 | |
end | |
return series | |
end | |
function partitions_generatingfn(N::Int64, denominations::Array{Int64, 1}) | |
#= | |
Solution to the coin-change problem, count the paritions with only specific summands'/denominations. | |
=# | |
if length(denominations) > 0 | |
return reduce((i, j) -> conv(i, j), Vector{Vector{Int8}}([fn_series(j, N) for j in denominations]))[N + 1] #repeated cauchy-product of the formal power series, works well upto 10000 | |
else | |
return reduce((i, j) -> conv(i, j), Vector{Vector{Int8}}([fn_series(j, N) for j in 1:N]))[N + 1] #some problems with this, doesn't work at all for values over 50 | |
end | |
end | |
#@time println(partitions_generatingfn(200, [1, 2, 5, 10, 20, 50, 100, 200])) | |
#@time println(partitions_recurrence(1000, 1000)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment