stla/hypergeomPQ.jl

## hypergeomPQ.jl
# this one works only for alpha = 1 !
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico, lambda)
#function _N(dico::Dict{Integer,Integer}, lambda::Vector{Integer})
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if(size(lambda,1) == 0)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(alpha, a, b, kappa, i)
    if ((isempty(kappa)) || (kappa[1] == 0))
        return 1
    end
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha*kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

function betaratio(kappa::Vector{Integer}, mu::Vector{Integer}, k::Int64, alpha::T)  where T<:Real
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(map(x -> x/(x+alpha-1), u))
    prod2 = prod(map(x -> (x+alpha)/x, v))
    prod3 = prod(map(x -> (x+alpha)/x, w))
    return alpha * prod1 * prod2 * prod3
end

function Hypergeo(m, alpha, a, b, x)
    function summation(i, z, j)
        function jack(k, beta, c, t, mu, Nmu)
            for i in (max(k,1) : sum(mu .> 0))
                if(size(mu,1) == i || mu[i] > mu[i+1])
                    d = Nmu
                    gamma = beta * betaratio(kappa, mu, i, alpha)
                    mu[i] = mu[i] - 1
                    println("$mu")
                    Nmu = _N(dico, mu) - 1
                    println("$mu")
                    if mu[i] > 0
                        jack(i, gamma, c+1, t, mu, Nmu)
                    else
                        if Nkappa > 1
                            f = 0
                            if(any(mu .> 0))
                                f = J[Nmu, t-1]
                            else
                                f = 1
                            end
                            J[Nkappa,t] = J[Nkappa,t] + gamma*x[t]^(c+1) * f
                        end
                    end
                    mu[i] = mu[i] + 1
                    Nmu = d
                end
            end
            if k == 0
                if Nkappa > 1
                    J[Nkappa,t] = J[Nkappa,t] + J[Nkappa,t-1]
                end
            else
                J[Nkappa,t] = J[Nkappa,t] + beta * x[t]^c * J[Nmu,t-1]
            end
        end # end jack
        r = 0
        if i == 1
            r = j
        else
            r = min(kappa[i-1],j)
        end
        Nkappa = 0
        for kappai in 1:r
            if(size(kappa,1) >= i)
                kappa[i] = kappai
            else
                push!(kappa, kappai)
            end
            println("$kappa")
#            kappa = vcat(kappa, [kappai])
            Nkappa = _N(dico, kappa) - 1
            println("Nkappa: $Nkappa")
            z = z * _T(alpha, a, b, kappa, i)
            if Nkappa > 1 && (size(kappa,1)==1 || kappa[2] == 0)
                J[Nkappa,1] = x[1] * (1+alpha*(kappa[1]-1)) * J[Nkappa-1,1]
            end
            for t in 2:n
                jack(0, 1.0, 0, t, kappa, Nkappa)
            end
            s <- s + z * J[Nkappa,n]
            if j > kappai && i < n
                summation(i+1, z, j-kappai)
            end
        end
        kappa[i] = 0
    end # end summation
    #
    #
    #
    n = size(x, 1)
    (dico, Pmn) = DictParts(m, n)
    s = 1
    J = zeros(Pmn, n)
    J[1,:] = cumsum(x)
    kappa = Int64[]
    summation(1, 1.0, m)
    return s
end

## hypergeomPQ10.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                x = Array{Int64}(undef,3,l)
                for j in 1:l
                    x[1:3,j] = [fin+j, manque-j, j]
                end
                NewLast = hcat(NewLast, x)
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        kappa::Vector{<:Integer}
    ) where {R<:Real}
    isComplex = eltype(a) <: Complex || eltype(b) <: Complex
    if isempty(kappa) || kappa[1] == 0
        return isComplex ? Complex(R(1)) : R(1)
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha * kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> count(>=(j), kappa), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f) ./ (l+h)) # "/" was an error !
    out = prod1/prod2 * prod3 * prod4
    return isinf(out) || isnan(out) ? (isComplex ? Complex(R(0)) : R(0)) : out
end

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, count(>=(i), lambda))
        end
    end
    return out
end

# ------------------------------------------------------------------------------
function hypergeomI(
        m::I,
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        n::I,
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {I<:Integer, R<:Real}
    function summation(
            i::Int64,
            z::T,
            j::I,
            kappa::Vector{I}
        ) where {R<:Real, T<:Union{R,Complex{R}}}
        function go(
                kappai::I,
                zz::T,
                s::T
            ) where {R<:Real, T<:Union{R,Complex{R}}}
            if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
                return s
            else
                kappap = vcat(kappa, [kappai])
                t = _T(alpha, a, b, filter(x -> x>0, kappap))
                zp = zz * x * (n - i + alpha * (kappai-1)) * t
                sp = j > kappai && i <= n ?
                     s + summation(i+1, zp, j-kappai, kappap) : s
                spp = sp + zp
                go(kappai + 1, zp, spp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), z, T(0))
    end # end summation --------------------------------------------------------
    isComplex =
        eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex
    the_one = isComplex ? Complex(R(1)) : R(1)
    return the_one + summation(0, the_one, m, I[])
end

# ------------------------------------------------------------------------------
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Vector{<:Union{R,Complex{R}}},
        alpha::Union{Nothing,R}=nothing
    )::Union{R,Complex{R}} where R<:Real
    if isnothing(alpha)
        alpha = R(2)
    end
    n = length(x)
    if all(x .== x[1])
        return hypergeomI(m, alpha, a, b, n, x[1])
    end
    function jack(
            k::I,
            beta::Union{R,Complex{R}},
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{<:Union{R,Complex{R}},2},
            kappa::Vector{I},
            nkappa::I
        ) where {I<:Integer}
        for i in ((max(k,1) : count(>(0), mu)))
            u = mu[i]
            if length(mu) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                filter!(>(0), mup)
                if length(mup) >= i && u > 1 # "!isempty(mup) &&" useless
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if !isempty(mup) # && mup[1] > 0 #any(mup .> 0)
                            jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack -------------------------------------------------------------
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
                jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), T(z), T(0))
    end # end summation --------------------------------------------------------
    the_one = eltype(a) <: Complex || eltype(b) <: Complex ||
                eltype(x) <: Complex ? Complex(R(1)) : R(1)
    (dico, Pmn) = DictParts(m, n)
    J = zeros(typeof(the_one), Pmn, n)
    J[1,:] = cumsum(x)
    return the_one + summation(0, the_one, m, Int64[], J)
end

# univariate -------------------------------------------------------------------
# both work
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {R<:Real}
    return hypergeomI(m, R(2), a, b, 1, x)
end

function hhypergeom(
        m::Integer,
        a::Vector{<:Union{R,T}},
        b::Vector{<:Union{R,T}},
        x::Union{R,T}
    ) where {R<:Real,T<:Complex{R}}
    return hypergeomI(m, R(2), a, b, 1, x)
end

## hypergeomPQ11.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            @inbounds manque = record[2]
            @inbounds l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                x = Array{Int64}(undef,3,l)
                for j in 1:l
                    @inbounds x[1:3,j] = [fin+j, manque-j, j]
                end
                NewLast = hcat(NewLast, x)
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    @inbounds dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        kappa::Vector{<:Integer}
    ) where {R<:Real}
    isComplex = eltype(a) <: Complex || eltype(b) <: Complex
    if isempty(kappa) || kappa[1] == 0
        return isComplex ? Complex(R(1)) : R(1)
    end
    i = lastindex(kappa)
    @inbounds c = kappa[i] - 1 - (i-1)/alpha
    @inbounds d = alpha * kappa[i] - i
    @inbounds s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> count(>=(j), kappa), s)
    g = e .+ 1
    ss = 1:(i-1)
    @inbounds f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f) ./ (l+h)) # "/" was an error !
    out = prod1/prod2 * prod3 * prod4
    return isinf(out) || isnan(out) ? (isComplex ? Complex(R(0)) : R(0)) : out
end

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    @inbounds muk = mu[k]
    t = k - alpha * muk
    @inbounds u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    @inbounds v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    @inbounds w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, count(>=(i), lambda))
        end
    end
    return out
end

# ------------------------------------------------------------------------------
function hypergeomI(
        m::I,
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        n::I,
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {I<:Integer, R<:Real}
    function summation(
            i::Int64,
            z::T,
            j::I,
            kappa::Vector{I}
        ) where {R<:Real, T<:Union{R,Complex{R}}}
        function go(
                kappai::I,
                zz::T,
                s::T
            ) where {R<:Real, T<:Union{R,Complex{R}}}
            if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
                return s
            else
                kappap = vcat(kappa, [kappai])
                t = _T(alpha, a, b, filter(x -> x>0, kappap))
                zp = zz * x * (n - i + alpha * (kappai-1)) * t
                sp = j > kappai && i <= n ?
                     s + summation(i+1, zp, j-kappai, kappap) : s
                spp = sp + zp
                go(kappai + 1, zp, spp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), z, T(0))
    end # end summation --------------------------------------------------------
    isComplex =
        eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex
    the_one = isComplex ? Complex(R(1)) : R(1)
    return the_one + summation(0, the_one, m, I[])
end

# ------------------------------------------------------------------------------
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Vector{<:Union{R,Complex{R}}},
        alpha::Union{Nothing,R}=nothing
    )::Union{R,Complex{R}} where R<:Real
    if isnothing(alpha)
        alpha = R(2)
    end
    n = length(x)
    if all(x .== x[1])
        return hypergeomI(m, alpha, a, b, n, x[1])
    end
    function jack(
            k::I,
            beta::Union{R,Complex{R}},
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{<:Union{R,Complex{R}},2},
            kappa::Vector{I},
            nkappa::I
        ) where {I<:Integer}
        for i in ((max(k,1) : count(>(0), mu)))
            @inbounds u = mu[i]
            if length(mu) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                filter!(>(0), mup)
                if length(mup) >= i && u > 1 # "!isempty(mup) &&" useless
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if !isempty(mup) # && mup[1] > 0 #any(mup .> 0)
                            @inbounds jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            @inbounds jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                @inbounds jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            @inbounds jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack -------------------------------------------------------------
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
                @inbounds jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), T(z), T(0))
    end # end summation --------------------------------------------------------
    the_one = eltype(a) <: Complex || eltype(b) <: Complex ||
                eltype(x) <: Complex ? Complex(R(1)) : R(1)
    (dico, Pmn) = DictParts(m, n)
    J = zeros(typeof(the_one), Pmn, n)
    J[1,:] = cumsum(x)
    return the_one + summation(0, the_one, m, Int64[], J)
end

# univariate -------------------------------------------------------------------
# both work
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {R<:Real}
    return hypergeomI(m, R(2), a, b, 1, x)
end

function hhypergeom(
        m::Integer,
        a::Vector{<:Union{R,T}},
        b::Vector{<:Union{R,T}},
        x::Union{R,T}
    ) where {R<:Real,T<:Complex{R}}
    return hypergeomI(m, R(2), a, b, 1, x)
end

## hypergeomPQ2.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico, lambda)
#function _N(dico::Dict{Integer,Integer}, lambda::Vector{Integer})
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if(size(lambda,1) == 0)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(alpha, a, b, kappa, i)
    if ((isempty(kappa)) || (kappa[1] == 0))
        return 1
    end
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha*kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

#function betaratio(kappa::Vector{Integer}, mu::Vector{Integer}, k::Int64, alpha::T)  where T<:Real
function betaratio(kappa, mu, k, alpha)
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(map(x -> x/(x+alpha-1), u))
    prod2 = prod(map(x -> (x+alpha)/x, v))
    prod3 = prod(map(x -> (x+alpha)/x, w))
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Array{Int64,1})
    out = Int64[]
    for i in 1:lambda[1]
        push!(out, size(filter(x -> x>=i, lambda), 1))
    end
    return out
end

function summation(a, b, x, dico, n, alpha, i, z, j, kappa, jarray)
    function go(kappai, zp, s)
        if i == n || i == 0 && kappai > j || i > 0 && kappai > min(last(kappa), j)
            return s
        end
        kappap = vcat(kappa, [kappai])
        nkappa = _N(dico, kappap) - 1
        zpp = zp * _T(alpha, a, b, kappap, lastindex(kappap)) # !!! _T à modifier
        lkappap = size(kappap, 1)
        if nkappa > 1 && (lkappap == 1 || kappap[2] == 0)
            entry = jarray[nkappa-1,1]
            jarray[nkappa,1] = x[1] * (1 + alpha*(kappap[1]-1)) * entry
        end
        for t in 2:n
            jarray = jack(alpha, x, dico, 0, 1.0, 0, t, kappap, jarray, kappap, nkappa)
        end
        entryp = jarray[nkappa,n]
        sp = s + zpp * entryp
        if j > kappai && i <= n
            spp = summation(a, b, x, dico, n, alpha, i+1, zpp, j-kappai, kappap, jarray)
            go(kappai+1, zpp, sp + spp)
        else
            go(kappai+1, zpp, sp)
        end
    end # end go
    go(1, z, 0.0)
end

function jack(alpha, x, dico, k, beta, c, t, mu, jarray, kappa, nkappa)
    for i in ((max(k,1) : sum(mu .> 0)))
        u = mu[i]
        if(size(mu,1) == i || u > mu[i+1])
            gamma = beta * betaratio(kappa, mu, i, alpha)
            mup = copy(mu)
            mup[i] = u - 1
            mup = mup[mup .> 0]
            nmu= _N(dico, mup) - 1
            if !isempty(mup) && size(mup,1) >= i && u > 1
                jarray = jack(alpha, x, dico, i, gamma, c+1, t, mup, jarray, kappa, nkappa)
            else
                if nkappa > 1
                    entryp = jarray[nkappa, t]
                    if any(mup .> 0)
                        entry = jarray[nmu, t-1]
                        jarray[nkappa,t] = entryp + gamma * entry * x[t]^(c+1)
                    else
                        jarray[nkappa,t] = entryp + gamma * x[t]^(c+1)
                    end
                end
            end
        end
    end
    entry1 = jarray[nkappa,t]
    if k == 0
        if nkappa > 1
            entry2 = jarray[nkappa, t-1]
            jarray[nkappa,t] = entry1 + entry2
        end
    else
        entry2 = jarray[_N(dico, mu)-1, t-1]
        jarray[nkappa,t] = entry1 + beta * x[t]^c * entry2
    end
    return jarray
end

function hypergeom(m, alpha, a, b, x)
    n = size(x, 1)
    (dico, Pmn) = DictParts(m, n)
    J = zeros(Pmn, n)
    J[1,:] = cumsum(x)
    kappa = Int64[]
    return 1.0 + summation(a, b, x, dico, n, alpha, 0, 1.0, m, kappa, J)
end

## hypergeomPQ3.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if(size(lambda,1) == 0)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(alpha::T, a::Vector{T}, b::Vector{T}, kappa::Vector{I}) where {I<:Integer, T<:Real}
    if ((isempty(kappa)) || (kappa[1] == 0))
        return 1
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha*kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

function betaratio(kappa::Vector{I}, mu::Vector{I}, k::I, alpha::T) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(map(x -> x/(x+alpha-1), u))
    prod2 = prod(map(x -> (x+alpha)/x, v))
    prod3 = prod(map(x -> (x+alpha)/x, w))
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{Integer})
    out = Int64[]
    for i in 1:lambda[1]
        push!(out, sum(lambda .>= i))# size(filter(x -> x>=i, lambda), 1))
    end
    return out
end


function hypergeom(m::Integer, alpha::T, a::Vector{T}, b::Vector{T}, x::Vector{T}) where T<:Real
    function jack(k::I, beta::T, c::I, t::Integer, mu::Vector{I}, jarray::Array{T,2}, kappa::Vector{I}, nkappa::I) where I<:Integer
        for i in ((max(k,1) : sum(mu .> 0)))
            u = mu[i]
            if(size(mu,1) == i || u > mu[i+1])
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                mup = mup[mup .> 0]
                #nmu= _N(dico, mup) - 1
                if !isempty(mup) && size(mup,1) >= i && u > 1
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        entryp = jarray[nkappa, t]
                        if any(mup .> 0)
                            entry = jarray[_N(dico, mup) - 1, t-1]
                            jarray[nkappa,t] = entryp + gamma * entry * x[t]^(c+1)
                        else
                            jarray[nkappa,t] = entryp + gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] = jarray[nkappa,t] + jarray[nkappa, t-1]
            end
        else
            entry = jarray[_N(dico, mu)-1, t-1]
            jarray[nkappa,t] = jarray[nkappa,t] + beta * x[t]^c * entry
        end
        return jarray
    end # end jack
    function summation(i::I, z::T, j::I, kappa::Vector{I}, jarray::Array{T,2}) where I<:Integer
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j || i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            #lkappap = size(kappap, 1)
            if nkappa > 1 && (size(kappap, 1) == 1 || kappap[2] == 0)
                entry = jarray[nkappa-1,1]
                jarray[nkappa,1] = x[1] * (1 + alpha*(kappap[1]-1)) * entry
            end
            for t in 2:n
                jarray = jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go
        go(I(1), T(z), T(0))
    end # end summation
    n = size(x, 1)
    (dico, Pmn) = DictParts(m, n)
    J = zeros(T, Pmn, n)
    J[1,:] = cumsum(x)
    return T(1) + summation(0, T(1), m, Int64[], J)
end

## hypergeomPQ4.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::T,
        a::Vector{T},
        b::Vector{T},
        kappa::Vector{I}
    ) where {I<:Integer, T<:Real}
    if isempty(kappa) || kappa[1] == 0
        return T(1)
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha*kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    # prod1 = prod(map(x -> x/(x+alpha-1), u))
    # prod2 = prod(map(x -> (x+alpha)/x, v))
    # prod3 = prod(map(x -> (x+alpha)/x, w))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, sum(lambda .>= i)) # size(filter(x -> x>=i, lambda), 1))
        end
    end
    return out
end


function hypergeom(
        m::Integer,
        alpha::T,
        a::Vector{T},
        b::Vector{T},
        x::Vector{T}
    ) where T<:Real
    function jack(
            k::I,
            beta::T,
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{T,2},
            kappa::Vector{I},
            nkappa::I
        ) where I<:Integer
        for i in ((max(k,1) : sum(mu .> 0)))
            u = mu[i]
            if size(mu,1) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                mup = mup[mup .> 0]
                #nmu= _N(dico, mup) - 1
                if !isempty(mup) && size(mup,1) >= i && u > 1
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if any(mup .> 0)
                            jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            jarray[nkappa,t] +=  beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where I<:Integer
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (size(kappap, 1) == 1 || kappap[2] == 0)
                jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go
        go(I(1), T(z), T(0))
    end # end summation
    n = size(x, 1)
    (dico, Pmn) = DictParts(m, n)
    J = zeros(T, Pmn, n)
    J[1,:] = cumsum(x)
    return T(1) + summation(0, T(1), m, Int64[], J)
end

## hypergeomPQ5.jl
# this one allows complex numbers

function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::R,
        a::Vector{T},
        b::Vector{T},
        kappa::Vector{I}
    ) where {R<:Real, I<:Integer, T<:Number}
    if isempty(kappa) || kappa[1] == 0
        return T(1)
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha * kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, sum(lambda .>= i))
        end
    end
    return out
end


function hypergeom(
        m::Integer,
        alpha::R,
        a::Vector{T},
        b::Vector{T},
        x::Vector{T}
    ) where {R<:Real, T<:Number}
    function jack(
            k::I,
            beta::T,
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{T,2},
            kappa::Vector{I},
            nkappa::I
        ) where I<:Integer
        for i in ((max(k,1) : sum(mu .> 0)))
            u = mu[i]
            if length(mu) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                mup = mup[mup .> 0]
                if !isempty(mup) && length(mup) >= i && u > 1
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if any(mup .> 0)
                            jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack -------------------------------------------------------------
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where I<:Integer
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
                jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), T(z), T(0))
    end # end summation --------------------------------------------------------
    n = length(x)
    (dico, Pmn) = DictParts(m, n)
    J = zeros(T, Pmn, n)
    J[1,:] = cumsum(x)
    return T(1) + summation(0, T(1), m, Int64[], J)
end

## hypergeomPQ6.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::R,
        a::Vector{T},
        b::Vector{T},
        kappa::Vector{I}
    ) where {R<:Real, I<:Integer, T<:Number}
    if isempty(kappa) || kappa[1] == 0
        return T(1)
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha * kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, sum(lambda .>= i))
        end
    end
    return out
end

function hypergeomI(
        m::I,
        alpha::R,
        a::Vector{T},
        b::Vector{T},
        n::I,
        x::T
    ) where {I<:Integer, R<:Real, T<:Number}
    function summation(i::Int64, z::T, j::I, kappa::Vector{I})
        function go(kappai::I, zz::T, s::T)
            if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
                return s
            else
                kappap = vcat(kappa, [kappai])
                t = _T(alpha, a, b, filter(x -> x>0, kappap))
                zp = zz * x * (n - i + alpha * (kappai-1)) * t
                sp = j > kappai && i <= n ?
                     s + summation(i+1, zp, j-kappai, kappap) : s
                spp = sp + zp
                go(kappai + 1, zp, spp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), z, T(0))
    end # end summation --------------------------------------------------------
    return 1 + summation(0, T(1), m, I[])
end


function hypergeom(
        m::Integer,
        alpha::R,
        a::Vector{T},
        b::Vector{T},
        x::Vector{T}
    ) where {R<:Real, T<:Number}
    n = length(x)
    if all(x .== x[1])
        return hypergeomI(m, alpha, a, b, n, x[1])
    end
    function jack(
            k::I,
            beta::T,
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{T,2},
            kappa::Vector{I},
            nkappa::I
        ) where I<:Integer
        for i in ((max(k,1) : sum(mu .> 0)))
            u = mu[i]
            if length(mu) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                mup = mup[mup .> 0]
                if !isempty(mup) && length(mup) >= i && u > 1
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if any(mup .> 0)
                            jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack -------------------------------------------------------------
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where I<:Integer
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
                jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), T(z), T(0))
    end # end summation --------------------------------------------------------
    (dico, Pmn) = DictParts(m, n)
    J = zeros(T, Pmn, n)
    J[1,:] = cumsum(x)
    return T(1) + summation(0, T(1), m, Int64[], J)
end

## hypergeomPQ7.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::R,
        a::Vector{T},
        b::Vector{T},
        kappa::Vector{I}
    ) where {R<:Real, I<:Integer, T<:Number}
    if isempty(kappa) || kappa[1] == 0
        return T(1)
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha * kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

# _T(1//1, [1//1], [2//1],[1,1,1]) = -1//0

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, sum(lambda .>= i))
        end
    end
    return out
end

# ------------------------------------------------------------------------------
function hypergeomI(
        m::I,
        alpha::R,
        a::Vector{T},
        b::Vector{T},
        n::I,
        x::T
    ) where {I<:Integer, R<:Real, T<:Number}
    function summation(i::Int64, z::T, j::I, kappa::Vector{I})
        function go(kappai::I, zz::T, s::T)
            if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
                return s
            else
                kappap = vcat(kappa, [kappai])
                #println("kappap: $kappap")
                t = _T(alpha, a, b, filter(x -> x>0, kappap))
                #println("t: $t") # -1//0 when alpha=1//1 a = [1///1] b = [2//1] kappap = [1,1,1]
                www = zz * x * (n - i + alpha * (kappai-1))
                #println("www: $www") # 0//1
                zp = zz * x * (n - i + alpha * (kappai-1)) * t
                #println("zp: $zp")
                sp = j > kappai && i <= n ?
                     s + summation(i+1, zp, j-kappai, kappap) : s
                spp = sp + zp
                go(kappai + 1, zp, spp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), z, T(0))
    end # end summation --------------------------------------------------------
    return 1 + summation(0, T(1), m, I[])
end
# hypergeomI(3, 1//1, [1//1], [2//1], 1, 1//2) error division per zero

# ------------------------------------------------------------------------------
function hypergeom(
        m::Integer,
        a::Vector{T},
        b::Vector{T},
        x::Vector{T},
        alpha::Union{Nothing,R}=nothing
    ) where {R<:Real, T<:Union{R,Complex{R}}}
    if isnothing(alpha)
        alpha = R(2)
    end
    n = length(x)
    if all(x .== x[1])
        return hypergeomI(m, alpha, a, b, n, x[1])
    end
    function jack(
            k::I,
            beta::T,
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{T,2},
            kappa::Vector{I},
            nkappa::I
        ) where I<:Integer
        for i in ((max(k,1) : sum(mu .> 0)))
            u = mu[i]
            if length(mu) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                mup = mup[mup .> 0]
                if !isempty(mup) && length(mup) >= i && u > 1
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if any(mup .> 0)
                            jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack -------------------------------------------------------------
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where I<:Integer
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
                jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), T(z), T(0))
    end # end summation --------------------------------------------------------
    (dico, Pmn) = DictParts(m, n)
    J = zeros(T, Pmn, n)
    J[1,:] = cumsum(x)
    return T(1) + summation(0, T(1), m, Int64[], J)
end

# univariate -------------------------------------------------------------------
# function hypergeom(
#         m::Integer,
#         a::Vector{T},
#         b::Vector{T},
#         x::T,
#         alpha::Union{Nothing,R}=nothing # useless - that does not depend on alpha !
#     ) where {R<:Real, T<:Union{R,Complex{R}}}
#     if isnothing(alpha)
#         alpha = R(1)
#     end
#     return hypergeomI(m, alpha, a, b, 1, x)
# end
function hypergeom(
        m::Integer,
        a::Vector{T},
        b::Vector{T},
        x::T
    ) where {R<:Real, T<:Union{R,Complex{R}}}
    return hypergeomI(m, R(2), a, b, 1, x)
end

## hypergeomPQ8.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        kappa::Vector{I}
    ) where {R<:Real, I<:Integer, T<:Union{Real,Complex{Real}}}
    if isempty(kappa) || kappa[1] == 0
        return T(1)
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha * kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f)/(l+h))
    return prod1/prod2 * prod3 * prod4
end

# _T(1//1, [1//1], [2//1],[1,1,1]) = -1//0

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, sum(lambda .>= i))
        end
    end
    return out
end

# ------------------------------------------------------------------------------
function hypergeomI(
        m::I,
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        n::I,
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {I<:Integer, R<:Real}
    function summation(i::Int64, z::T, j::I, kappa::Vector{I}) where {R<:Real,T<:Union{R,Complex{R}}}
        function go(kappai::I, zz::T, s::T) where {R<:Real,T<:Union{R,Complex{R}}}
            if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
                return s
            else
                kappap = vcat(kappa, [kappai])
                #println("kappap: $kappap")
                t = _T(alpha, a, b, filter(x -> x>0, kappap))
                #println("t: $t") # -1//0 when alpha=1//1 a = [1///1] b = [2//1] kappap = [1,1,1]
                www = zz * x * (n - i + alpha * (kappai-1))
                #println("www: $www") # 0//1
                zp = zz * x * (n - i + alpha * (kappai-1)) * t
                #println("zp: $zp")
                sp = j > kappai && i <= n ?
                     s + summation(i+1, zp, j-kappai, kappap) : s
                spp = sp + zp
                go(kappai + 1, zp, spp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), z, T(0))
    end # end summation --------------------------------------------------------
    return 1.0 + summation(0, 1.0, m, I[])
end
# hypergeomI(3, 1//1, [1//1], [2//1], 1, 1//2) error division per zero

# ------------------------------------------------------------------------------
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Vector{T},
        alpha::Union{Nothing,R}=nothing
    )::Union{R,Complex{R}} where {R<:Real, T<:Union{R,Complex{R}}}
    if isnothing(alpha)
        alpha = R(2)
    end
    n = length(x)
    if all(x .== x[1])
        return hypergeomI(m, alpha, a, b, n, x[1])
    end
    function jack(
            k::I,
            beta::Union{R,Complex{R}},
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{<:Union{R,Complex{R}},2},
            kappa::Vector{I},
            nkappa::I
        ) where {I<:Integer}
        for i in ((max(k,1) : sum(mu .> 0)))
            u = mu[i]
            if length(mu) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                mup = mup[mup .> 0]
                if !isempty(mup) && length(mup) >= i && u > 1
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if any(mup .> 0)
                            jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack -------------------------------------------------------------
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
        #z = convert(typeof(x),z)
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
                jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), T(z), T(0))
    end # end summation --------------------------------------------------------
    the_one = eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex ? Complex(R(1)) : R(1)
    the_type = typeof(the_one)
    (dico, Pmn) = DictParts(m, n)
    #println("theone: $the_one")
    J = zeros(the_type, Pmn, n)
    J[1,:] = cumsum(x)
    return convert(the_type,1) + summation(0, the_one, m, Int64[], J)
end

# univariate -------------------------------------------------------------------
# both work
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {R<:Real}
    return hypergeomI(m, R(2), a, b, 1, x)
end

function hhypergeom(
        m::Integer,
        a::Vector{<:Union{R,T}},
        b::Vector{<:Union{R,T}},
        x::Union{R,T}
    ) where {R<:Real,T<:Complex{R}}
    return hypergeomI(m, R(2), a, b, 1, x)
end

## hypergeomPQ9.jl
function DictParts(m::Integer, n::Integer)
    D = Dict{Int64,Int64}()
    Last = hcat([0,m,m])
    fin = 0
    for i in 1:n
        NewLast = Array{Int64}(undef,3,0)
        for record in eachcol(Last)
            manque = record[2]
            l = min(manque, record[3])
            if l > 0
                D[record[1]+1] = fin + 1
                for j in 1:l
                    NewLast = hcat(NewLast, [fin+j, manque-j, j])
                end
                fin += l
            end
        end
        Last = NewLast
    end
    return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
    lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
    if isempty(lambda)
        return 1
    end
    dico[_N(dico, lambda[1:end-1])] + last(lambda)
end # il faut retirer 1 au résultat

function _T(
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        kappa::Vector{<:Integer}
    ) where {R<:Real}
    isComplex = eltype(a) <: Complex || eltype(b) <: Complex
    if isempty(kappa) || kappa[1] == 0
        return isComplex ? Complex(R(1)) : R(1)
    end
    i = lastindex(kappa)
    c = kappa[i] - 1 - (i-1)/alpha
    d = alpha * kappa[i] - i
    s = collect(1:(kappa[i]-1))
    e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
    g = e .+ 1
    ss = 1:(i-1)
    f = alpha .* kappa[ss] - collect(ss) .- d
    h = f .+ alpha
    l = h .* f
    prod1 = prod(a .+ c)
    prod2 = prod(b .+ c)
    prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
    prod4 = prod((l-f) ./ (l+h)) # "/" was an error !
    out = prod1/prod2 * prod3 * prod4
    return isinf(out) || isnan(out) ? (isComplex ? Complex(R(0)) : R(0)) : out
end

function betaratio(
        kappa::Vector{I},
        mu::Vector{I},
        k::I,
        alpha::T
    ) where {I<:Integer, T<:Real}
    muk = mu[k]
    t = k - alpha * muk
    u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
    v = map(i -> t - i + alpha * mu[i], 1:(k-1))
    muPrime = dualPartition(mu)
    w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
    prod1 = prod(u ./ (u .+ alpha .- 1))
    prod2 = prod((v .+ alpha) ./ v)
    prod3 = prod((w .+ alpha) ./ w)
    return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where T<:Integer
    out = T[]
    if !isempty(lambda)
        for i in 1:lambda[1]
            push!(out, sum(lambda .>= i))
        end
    end
    return out
end

# ------------------------------------------------------------------------------
function hypergeomI(
        m::I,
        alpha::R,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        n::I,
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {I<:Integer, R<:Real}
    function summation(
            i::Int64,
            z::T,
            j::I,
            kappa::Vector{I}
        ) where {R<:Real, T<:Union{R,Complex{R}}}
        function go(
                kappai::I,
                zz::T,
                s::T
            ) where {R<:Real, T<:Union{R,Complex{R}}}
            if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i],j)
                return s
            else
                kappap = vcat(kappa, [kappai])
                t = _T(alpha, a, b, filter(x -> x>0, kappap))
                zp = zz * x * (n - i + alpha * (kappai-1)) * t
                sp = j > kappai && i <= n ?
                     s + summation(i+1, zp, j-kappai, kappap) : s
                spp = sp + zp
                go(kappai + 1, zp, spp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), z, T(0))
    end # end summation --------------------------------------------------------
    isComplex =
        eltype(a) <: Complex || eltype(b) <: Complex || eltype(x) <: Complex
    the_one = isComplex ? Complex(R(1)) : R(1)
    return the_one + summation(0, the_one, m, I[])
end

# ------------------------------------------------------------------------------
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Vector{<:Union{R,Complex{R}}},
        alpha::Union{Nothing,R}=nothing
    )::Union{R,Complex{R}} where R<:Real
    if isnothing(alpha)
        alpha = R(2)
    end
    n = length(x)
    if all(x .== x[1])
        return hypergeomI(m, alpha, a, b, n, x[1])
    end
    function jack(
            k::I,
            beta::Union{R,Complex{R}},
            c::I,
            t::Integer,
            mu::Vector{I},
            jarray::Array{<:Union{R,Complex{R}},2},
            kappa::Vector{I},
            nkappa::I
        ) where {I<:Integer}
        for i in ((max(k,1) : sum(mu .> 0)))
            u = mu[i]
            if length(mu) == i || u > mu[i+1]
                gamma = beta * betaratio(kappa, mu, i, alpha)
                mup = copy(mu)
                mup[i] = u - 1
                mup = mup[mup .> 0]
                if length(mup) >= i && u > 1 # "!isempty(mup) &&" useless
                    jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
                else
                    if nkappa > 1
                        if !isempty(mup) && mup[1] > 0 #any(mup .> 0)
                            jarray[nkappa,t] +=
                              gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
                        else
                            jarray[nkappa,t] += gamma * x[t]^(c+1)
                        end
                    end
                end
            end
        end
        if k == 0
            if nkappa > 1
                jarray[nkappa,t] += jarray[nkappa, t-1]
            end
        else
            jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
        end
        return jarray
    end # end jack -------------------------------------------------------------
    function summation(
            i::I,
            z::T,
            j::I,
            kappa::Vector{I},
            jarray::Array{T,2}
        ) where {I<:Integer, R<:Real, T<:Union{R,Complex{R}}}
        function go(kappai::I, zp::T, s::T)
            if i == n || i == 0 && kappai > j ||
                i > 0 && kappai > min(last(kappa), j)
                return s
            end
            kappap = vcat(kappa, [kappai])
            nkappa = _N(dico, kappap) - 1
            zpp = zp * _T(alpha, a, b, kappap)
            if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
                jarray[nkappa,1] =
                  x[1] * (1 + alpha*(kappap[1]-1)) * jarray[nkappa-1,1]
            end
            for t in 2:n
                jarray =
                    jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
            end
            sp = s + zpp * jarray[nkappa,n]
            if j > kappai && i <= n
                spp = summation(i+1, zpp, j-kappai, kappap, jarray)
                go(kappai+1, zpp, sp + spp)
            else
                go(kappai+1, zpp, sp)
            end
        end # end go -----------------------------------------------------------
        go(I(1), T(z), T(0))
    end # end summation --------------------------------------------------------
    the_one = eltype(a) <: Complex || eltype(b) <: Complex ||
                eltype(x) <: Complex ? Complex(R(1)) : R(1)
    (dico, Pmn) = DictParts(m, n)
    J = zeros(typeof(the_one), Pmn, n)
    J[1,:] = cumsum(x)
    return the_one + summation(0, the_one, m, Int64[], J)
end

# univariate -------------------------------------------------------------------
# both work
function hypergeom(
        m::Integer,
        a::Vector{<:Union{R,Complex{R}}},
        b::Vector{<:Union{R,Complex{R}}},
        x::Union{R,Complex{R}}
    )::Union{R,Complex{R}} where {R<:Real}
    return hypergeomI(m, R(2), a, b, 1, x)
end

function hhypergeom(
        m::Integer,
        a::Vector{<:Union{R,T}},
        b::Vector{<:Union{R,T}},
        x::Union{R,T}
    ) where {R<:Real,T<:Complex{R}}
    return hypergeomI(m, R(2), a, b, 1, x)
end
	# this one works only for alpha = 1 !
	function DictParts(m::Integer, n::Integer)
	D = Dict{Int64,Int64}()
	Last = hcat([0,m,m])
	fin = 0
	for i in 1:n
	NewLast = Array{Int64}(undef,3,0)
	for record in eachcol(Last)
	manque = record[2]
	l = min(manque, record[3])
	if l > 0
	D[record[1]+1] = fin + 1
	for j in 1:l
	NewLast = hcat(NewLast, [fin+j, manque-j, j])
	end
	fin += l
	end
	end
	Last = NewLast
	end
	return (D, fin)
	end

	function _N(dico, lambda)
	#function _N(dico::Dict{Integer,Integer}, lambda::Vector{Integer})
	lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
	if(size(lambda,1) == 0)
	return 1
	end
	dico[_N(dico, lambda[1:end-1])] + last(lambda)
	end # il faut retirer 1 au résultat

	function _T(alpha, a, b, kappa, i)
	if ((isempty(kappa)) \|\| (kappa[1] == 0))
	return 1
	end
	c = kappa[i] - 1 - (i-1)/alpha
	d = alpha*kappa[i] - i
	s = collect(1:(kappa[i]-1))
	e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
	g = e .+ 1
	ss = 1:(i-1)
	f = alpha .* kappa[ss] - collect(ss) .- d
	h = f .+ alpha
	l = h .* f
	prod1 = prod(a .+ c)
	prod2 = prod(b .+ c)
	prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
	prod4 = prod((l-f)/(l+h))
	return prod1/prod2 * prod3 * prod4
	end

	function betaratio(kappa::Vector{Integer}, mu::Vector{Integer}, k::Int64, alpha::T) where T<:Real
	muk = mu[k]
	t = k - alpha * muk
	u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
	v = map(i -> t - i + alpha * mu[i], 1:(k-1))
	muPrime = dualPartition(mu)
	w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
	prod1 = prod(map(x -> x/(x+alpha-1), u))
	prod2 = prod(map(x -> (x+alpha)/x, v))
	prod3 = prod(map(x -> (x+alpha)/x, w))
	return alpha * prod1 * prod2 * prod3
	end

	function Hypergeo(m, alpha, a, b, x)
	function summation(i, z, j)
	function jack(k, beta, c, t, mu, Nmu)
	for i in (max(k,1) : sum(mu .> 0))
	if(size(mu,1) == i \|\| mu[i] > mu[i+1])
	d = Nmu
	gamma = beta * betaratio(kappa, mu, i, alpha)
	mu[i] = mu[i] - 1
	println("$mu")
	Nmu = _N(dico, mu) - 1
	println("$mu")
	if mu[i] > 0
	jack(i, gamma, c+1, t, mu, Nmu)
	else
	if Nkappa > 1
	f = 0
	if(any(mu .> 0))
	f = J[Nmu, t-1]
	else
	f = 1
	end
	J[Nkappa,t] = J[Nkappa,t] + gammax[t]^(c+1) f
	end
	end
	mu[i] = mu[i] + 1
	Nmu = d
	end
	end
	if k == 0
	if Nkappa > 1
	J[Nkappa,t] = J[Nkappa,t] + J[Nkappa,t-1]
	end
	else
	J[Nkappa,t] = J[Nkappa,t] + beta * x[t]^c * J[Nmu,t-1]
	end
	end # end jack
	r = 0
	if i == 1
	r = j
	else
	r = min(kappa[i-1],j)
	end
	Nkappa = 0
	for kappai in 1:r
	if(size(kappa,1) >= i)
	kappa[i] = kappai
	else
	push!(kappa, kappai)
	end
	println("$kappa")
	# kappa = vcat(kappa, [kappai])
	Nkappa = _N(dico, kappa) - 1
	println("Nkappa: $Nkappa")
	z = z * _T(alpha, a, b, kappa, i)
	if Nkappa > 1 && (size(kappa,1)==1 \|\| kappa[2] == 0)
	J[Nkappa,1] = x[1] * (1+alpha(kappa[1]-1)) J[Nkappa-1,1]
	end
	for t in 2:n
	jack(0, 1.0, 0, t, kappa, Nkappa)
	end
	s <- s + z * J[Nkappa,n]
	if j > kappai && i < n
	summation(i+1, z, j-kappai)
	end
	end
	kappa[i] = 0
	end # end summation
	#
	#
	#
	n = size(x, 1)
	(dico, Pmn) = DictParts(m, n)
	s = 1
	J = zeros(Pmn, n)
	J[1,:] = cumsum(x)
	kappa = Int64[]
	summation(1, 1.0, m)
	return s
	end
	# this one allows complex numbers

	function DictParts(m::Integer, n::Integer)
	D = Dict{Int64,Int64}()
	Last = hcat([0,m,m])
	fin = 0
	for i in 1:n
	NewLast = Array{Int64}(undef,3,0)
	for record in eachcol(Last)
	manque = record[2]
	l = min(manque, record[3])
	if l > 0
	D[record[1]+1] = fin + 1
	for j in 1:l
	NewLast = hcat(NewLast, [fin+j, manque-j, j])
	end
	fin += l
	end
	end
	Last = NewLast
	end
	return (D, fin)
	end

	function _N(dico::Dict{I,I}, lambda::Vector{I}) where I<:Integer
	lambda = filter(x -> x>0, lambda) # ne pas faire filter! car ça se propage
	if isempty(lambda)
	return 1
	end
	dico[_N(dico, lambda[1:end-1])] + last(lambda)
	end # il faut retirer 1 au résultat

	function _T(
	alpha::R,
	a::Vector{T},
	b::Vector{T},
	kappa::Vector{I}
	) where {R<:Real, I<:Integer, T<:Number}
	if isempty(kappa) \|\| kappa[1] == 0
	return T(1)
	end
	i = lastindex(kappa)
	c = kappa[i] - 1 - (i-1)/alpha
	d = alpha * kappa[i] - i
	s = collect(1:(kappa[i]-1))
	e = d .- alpha .* s + map(j -> sum(kappa .>= j), s)
	g = e .+ 1
	ss = 1:(i-1)
	f = alpha .* kappa[ss] - collect(ss) .- d
	h = f .+ alpha
	l = h .* f
	prod1 = prod(a .+ c)
	prod2 = prod(b .+ c)
	prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
	prod4 = prod((l-f)/(l+h))
	return prod1/prod2 * prod3 * prod4
	end

	function betaratio(
	kappa::Vector{I},
	mu::Vector{I},
	k::I,
	alpha::T
	) where {I<:Integer, T<:Real}
	muk = mu[k]
	t = k - alpha * muk
	u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
	v = map(i -> t - i + alpha * mu[i], 1:(k-1))
	muPrime = dualPartition(mu)
	w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
	prod1 = prod(u ./ (u .+ alpha .- 1))
	prod2 = prod((v .+ alpha) ./ v)
	prod3 = prod((w .+ alpha) ./ w)
	return alpha * prod1 * prod2 * prod3
	end

	function dualPartition(lambda::Vector{T}) where T<:Integer
	out = T[]
	if !isempty(lambda)
	for i in 1:lambda[1]
	push!(out, sum(lambda .>= i))
	end
	end
	return out
	end


	function hypergeom(
	m::Integer,
	alpha::R,
	a::Vector{T},
	b::Vector{T},
	x::Vector{T}
	) where {R<:Real, T<:Number}
	function jack(
	k::I,
	beta::T,
	c::I,
	t::Integer,
	mu::Vector{I},
	jarray::Array{T,2},
	kappa::Vector{I},
	nkappa::I
	) where I<:Integer
	for i in ((max(k,1) : sum(mu .> 0)))
	u = mu[i]
	if length(mu) == i \|\| u > mu[i+1]
	gamma = beta * betaratio(kappa, mu, i, alpha)
	mup = copy(mu)
	mup[i] = u - 1
	mup = mup[mup .> 0]
	if !isempty(mup) && length(mup) >= i && u > 1
	jarray = jack(i, gamma, c+1, t, mup, jarray, kappa, nkappa)
	else
	if nkappa > 1
	if any(mup .> 0)
	jarray[nkappa,t] +=
	gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c+1)
	else
	jarray[nkappa,t] += gamma * x[t]^(c+1)
	end
	end
	end
	end
	end
	if k == 0
	if nkappa > 1
	jarray[nkappa,t] += jarray[nkappa, t-1]
	end
	else
	jarray[nkappa,t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
	end
	return jarray
	end # end jack -------------------------------------------------------------
	function summation(
	i::I,
	z::T,
	j::I,
	kappa::Vector{I},
	jarray::Array{T,2}
	) where I<:Integer
	function go(kappai::I, zp::T, s::T)
	if i == n \|\| i == 0 && kappai > j \|\|
	i > 0 && kappai > min(last(kappa), j)
	return s
	end
	kappap = vcat(kappa, [kappai])
	nkappa = _N(dico, kappap) - 1
	zpp = zp * _T(alpha, a, b, kappap)
	if nkappa > 1 && (length(kappap) == 1 \|\| kappap[2] == 0)
	jarray[nkappa,1] =
	x[1] * (1 + alpha(kappap[1]-1)) jarray[nkappa-1,1]
	end
	for t in 2:n
	jarray =
	jack(I(0), T(1), I(0), t, kappap, jarray, kappap, nkappa)
	end
	sp = s + zpp * jarray[nkappa,n]
	if j > kappai && i <= n
	spp = summation(i+1, zpp, j-kappai, kappap, jarray)
	go(kappai+1, zpp, sp + spp)
	else
	go(kappai+1, zpp, sp)
	end
	end # end go -----------------------------------------------------------
	go(I(1), T(z), T(0))
	end # end summation --------------------------------------------------------
	n = length(x)
	(dico, Pmn) = DictParts(m, n)
	J = zeros(T, Pmn, n)
	J[1,:] = cumsum(x)
	return T(1) + summation(0, T(1), m, Int64[], J)
	end