Evaluating the Bell polynomials

Hello,

I just wrote a function, based on a recurrence relation that can be found on Wikipedia.

Do you have any suggestion whether it is useful for this package, or else if this code could be a useful addition elsewhere? Any (other) comment of course welcome too.

Cheers,

Evaluate the incomplete Bell polynomial B_{N,K}(x)
where x is an array containing the abscissa values [x_1, ... x{n-k+1}]

Definition and background information, see 
https://en.wikipedia.org/wiki/Bell_polynomials#Definitions
https://en.wikipedia.org/wiki/Bell_polynomials#Properties
https://en.wikipedia.org/wiki/Bell_polynomials#Table_of_values
"""
evaluate_incomplete_Bell_poly = function(N,K,all_x)
    # see formula C2
    # let us construct the triangular table
    if K > N
        return 0.
    end
    all_B_old = zeros(N+1,1)
    all_B_old[1] = 1.
    all_B_new = zeros(N+1,1)
    all_B_new[1] = 1.
    for k=1:K
        for n=1:N
            all_B_new[n+1] = sum([binomial(n-1,i-1)*all_x[i]*all_B_old[n-i+1] for i=1:n-k+1])
        end
        all_B_old .= all_B_new
    end
    return all_B_new[end]
end

## test the function 
x = randn(11,1)

@assert evaluate_incomplete_Bell_poly(0,0,x) ≈ 1.

@assert evaluate_incomplete_Bell_poly(1,0,x) ≈ 0.
@assert evaluate_incomplete_Bell_poly(3,0,x)≈ 0.

@assert evaluate_incomplete_Bell_poly(0,1,x) ≈ 0.
@assert evaluate_incomplete_Bell_poly(0,3,x) ≈ 0.

@assert evaluate_incomplete_Bell_poly(1,1,x) ≈ x[1]
@assert evaluate_incomplete_Bell_poly(2,2,x) ≈ x[1]^2
@assert evaluate_incomplete_Bell_poly(3,3,x) ≈ x[1]^3

@assert evaluate_incomplete_Bell_poly(6,4,x) ≈ 45*x[1]^2*x[2]^2 + 20*x[1]^3*x[3]
@assert evaluate_incomplete_Bell_poly(5,4,x) ≈ 10*x[1]^3*x[2]
@assert evaluate_incomplete_Bell_poly(6,3,x) ≈ 15*x[2]^3 + 60*x[1]*x[2]*x[3] + 15*x[1]^2*x[4]

# use some results in (Abbas, M.; Bouroubi, S. (2005). "On new identities for Bell's polynomial". Discrete Math. 293 (1–3)) as reference
# see https://vinar.vin.bg.ac.rs//bitstream/id/12791/4374.pdf
@assert evaluate_incomplete_Bell_poly(9,7,x) ≈ 378*x[1]^5*x[2]^2 + 84*x[1]^6*x[3]
@assert evaluate_incomplete_Bell_poly(10,7,x) ≈ 3150*x[1]^4*x[2]^3 + 2520*x[1]^5*x[2]*x[3] + 210*x[1]^6*x[4]
@assert evaluate_incomplete_Bell_poly(11,7,x) ≈ 17325*x[1]^3*x[2]^4 + 34650*x[1]^4*x[2]^2*x[3] + + 4620*x[1]^5*x[3]^2 + 6930*x[1]^5*x[2]*x[4] + 462*x[1]^6*x[5]

println("tests passed")

JuliaMath / Combinatorics.jl

Evaluating the Bell polynomials #139