Symbolics: get the coefficients of the polynomial binomial(x, c)

General Usage
1

How to get a polynomial expression from, e.g., binomial(x, 3)?

julia> using Symbolics

julia> @variables x
1-element Vector{Num}:
 x

julia> binomial(x, 3)
binomial(x, 3)

julia> simplify(ans)
binomial(x, 3)
2

binomial(x, 3) is not a polynomial. It corresponds to the quantity \frac{x!}{3! (x-3)!} (which could be extended to non-integer x via the gamma function, a transcendental function of x).

What polynomial do you want?

3

Trying to reverse engineer the question, maybe the meaning is for:

julia> @variables x
1-element Vector{Num}:
 x

julia> p = expand((1+x)^3)
1 + x^3 + 3x + 3(x^2)

That’s a polynomial, and is related to binomials (coefficients are the binomial numbers). To obtain the sequence of binomial coefficients with Symbolics, there must be a better way, but:

julia> pp = groebner_basis([p])[1]
(x^3) + 3(x^2) + 3x + 1

julia> using DynamicPolynomials

julia> map(coefficient, pp.p)
4-element Vector{Int64}:
 1
 3
 3
 1

An alternative is to use the Polynomials package:

julia> using Polynomials

julia> p = Polynomial([1,1])
Polynomial(1 + x)

julia> p^3
Polynomial(1 + 3*x + 3*x^2 + x^3)

julia> coeffs(p^3)
4-element Vector{Int64}:
 1
 3
 3
 1

All this… if the question was regarding this. For the Gamma function interpretation (from Wikipedia for Binomial Coefficient):

gamma_func
gamma_func1945×389 103 KB

1 Like
4

binomial(x, 3) should be \binom{x}{3} = \frac{x!}{3! (x-3)!} = {\frac{1}{3!} x (x - 1) (x - 2)} = {\frac{1}{3!} (x^3 - 3x^2 + 2x)}

That’s what I meant.

Hmmm, even the following fails to do as expected:

julia> factorial(x) / (factorial(3) * factorial(x - 3))
factorial(x) / (6factorial(x - 3))

julia> factorial(x) / factorial(x - 3)
factorial(x) / factorial(x - 3)

julia> simplify(ans)
factorial(x) / factorial(x - 3)
5

These simplifications are only valid if x is a natural number, by default Symbolics variables are Real. If x were negative, fractional, etc, the expression would look more like the one from Wikipedia above. If you create x as a Num{Integer} instead of a Num{Real}, I think you can define custom dispatches to get the expression you want, but I haven’t actually tried to do this.

1 Like
6

I’d say they’re perfectly valid as-is, but the other direction perhaps wouldn’t be valid, because a polynomial is defined over all reals, which include the integers.

Although I admit that these technical limitations may be reasonable in this case.

7

Oh, I mean the simplification x! == prod(1:x) is only valid for x∈𝒩

1 Like
8

How to do this? I can’t figure it out.

9

Me neither, I thought It was in the docs but I can’t find it. I think I have a fundamental misunderstanding of how Symbolics does this.

1 Like
10 
Bincoeffs(n)=coeffs(mapreduce( x->1//(x+1)*Polynomial([-x,1]), *, 0:n-1 ))
11

Oh right, this generalization. Yes, this is a indeed a polynomial for integer k. Right now, however, Julia only defines the binomial function for integers:

julia> methods(binomial)
# 3 methods for generic function "binomial":
[1] binomial(n::BigInt, k::UInt64) in Base.GMP at gmp.jl:684
[2] binomial(n::BigInt, k::Integer) in Base.GMP at gmp.jl:685
[3] binomial(n::T, k::T) where T<:Integer in Base at intfuncs.jl:1050

Arguably, it should have a method binomial(x::Number, k::Integer) (binomial(x,k) for non-integer x by stevengj · Pull Request #48124 · JuliaLang/julia · GitHub). If we add such a method:

function binomial(x::Number, k::Integer)
    k < 0 && throw(DomainError("k=$k must be non-negative"))
    k == 0 && return oneunit(x)/one(k)
    # we don't use prod(i -> (x-i+1), 1:k) / factorial(k),
    # and instead divide each term by i, to avoid spurious overflow.
    return prod(i -> (x-i+1)/i, OneTo(k))
end

along with a method Base.binomial(x::Num, k::Integer) = invoke(binomial, Tuple{Number, Integer}, x, k) to eliminate a method ambiguity in Symbolics.jl, then it works:

julia> binomial(x, 3)
x*((1//3)*x - (2//3))*((1//2)*x - (1//2))
2 Likes
12

The Julia factorial function is not also not defined for non-integer x. This definition definitely would require the Gamma function, which is not provided in the stdlib. At one point factorial(::Number) was provided by the SpecialFunctions.jl package (which was originally in Base but was split off for Julia 1.0), but was removed because it is type piracy.

So, if you want to use the factorial(x) for non-integer x, you need to explicitly call the gamma(x+1) function from SpecialFunctions.jl. However, looks like Symbolics.jl does not know how to simplify the ratio of two gamma functions into a polynomial.