Partition function, coin-change problem, generating integer partitions in (anti-)/lexicographic order.

partitions.jl

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))