R: Multinomial Formula in R

Question

I am working with the R programming language.

I am trying to write an R function that can fully expand a mathematical expression (in symbolic form) of the form: (a+b+c+....)^n

In the above expression, n is an integer. I want the final answer to be in terms of a,b,c... (i.e. I am not looking to perform arithmetic calculations, but rather expand an expression symbolically.

As I understand, this corresponds to the following mathematical equation: https://en.wikipedia.org/wiki/Multinomial_theorem

I was able to get the following function to work that expands a binomial expression:

expandBinom <- function(a, b, n) {
    sapply(0:n, function(k) {
        choose(n, k) * a^(n-k) * b^k
    })
}

The above function calculates the coefficients of the terms being expanded. For example, I can verify that (a+b)^2 = a^2 + 2*ab + b^2:

> expandBinom(1,1,2)
[1] 1 2 1

But I do not know how to generalize this above function to the multinomial case.

Is it possible to write a function like this for multinomial expressions - and actually have the expression appear in the form of (a^2 + 2ab^2 + c..) instead of just (1 2 1)

Thanks!

References:

• https://en.wikipedia.org/wiki/Pascal%27s_triangle

Answer 1

For symbolic computations

Probably you can try Ryacas package

library(Ryacas)

f <- function(n, ...) {
    as_r(yac_str(sprintf("Simplify(Expand((%s)^%s))", paste0(c(...), collapse = "+"), n)))
}

and you will obtain, for example

> f(5, "a", "b", "c")
expression(a^5 + 5 * a^4 * b + 5 * a^4 * c + 10 * a^3 * b^2 + 
    20 * a^3 * b * c + 10 * a^3 * c^2 + 10 * a^2 * b^3 + 30 * 
    a^2 * b^2 * c + 30 * a^2 * b * c^2 + 10 * a^2 * c^3 + 5 *
    a * b^4 + 20 * a * b^3 * c + 30 * a * b^2 * c^2 + 20 * a *
    b * c^3 + 5 * a * c^4 + b^5 + 5 * b^4 * c + 10 * b^3 * c^2 +
    10 * b^2 * c^3 + 5 * b * c^4 + c^5)

---

For numerical computations (coefficients only)

If you only care about the numerical coefficients of polynomial (e.g., a polynomial of x), you can try convolve + Reduce

h <- function(n, coeffs) {
    Reduce(\(x, y) convolve(x, rev(y), type = "open"), rep(list(coeffs), n))
}

or its recursion variant

h2 <- function(n, coeffs) {
    if (n == 1) {
        return(coeffs)
    }
    convolve(Recall(n - 1, coeffs), rev(coeffs), type = "open")
}

and you will see

# (1+x)^2 = 1 + 2x + x^2
> h(2, c(1, 1)) 
[1] 1 2 1

# (1+2x)^3 = 1 + 6x + 12x^2 + 8x^3
> h(3, c(1, 2)) 
[1]  1  6 12  8

# (1 + 4x + 5x^2) = 1 + 16x + 116x^2 + 496x^3 + 1366x^4 + 2480x^5 + 2900x^6 + 2000x^7 + 625x^8
> h(4, c(1, 4, 5)) 
[1]    1   16  116  496 1366 2480 2900 2000  625

Answer 2

You can use giacR (not on CRAN yet):

library(giacR)
giac <- Giac$new()
giac$execute(("expand((a + b + c)^5)"))
# "a^5+b^5+c^5+5*a*b^4+5*a*c^4+5*b*c^4+10*a^2*b^3+10*a^2*c^3+10*a^3*b^2+10*a^3*c^2+5*a^4*b+5*a^4*c+10*b^2*c^3+10*b^3*c^2+5*b^4*c+20*a*b*c^3+30*a*b^2*c^2+20*a*b^3*c+30*a^2*b*c^2+30*a^2*b^2*c+20*a^3*b*c"