R: Полиномиальная формула в R
Я работаю с языком программирования R.
Я пытаюсь написать функцию R, которая может полностью расширить математическое выражение (в символьной форме) вида: (a+b+c+....)^n
В приведенном выше выражении n является целым числом. Я хочу, чтобы окончательный ответ был в терминах a,b,c... (т.е. я не собираюсь выполнять арифметические вычисления, а скорее символически расширяю выражение.
Насколько я понимаю, это соответствует следующему математическому уравнению: https://en.wikipedia.org/wiki/Мультиномиальная_теорема
Мне удалось заставить работать следующую функцию, которая расширяет биномиальное выражение:
expandBinom <- function(a, b, n) {
sapply(0:n, function(k) {
choose(n, k) * a^(n-k) * b^k
})
}
Приведенная выше функция вычисляет коэффициенты расширяемых членов. Например, я могу проверить, что (a+b)^2 = a^2 + 2*ab + b^2:
> expandBinom(1,1,2)
[1] 1 2 1
Но я не знаю, как обобщить эту вышеприведенную функцию на полиномиальный случай.
Можно ли написать такую функцию для полиномиальных выражений, и на самом деле выражение будет отображаться в форме (a^2 + 2ab^2 + c..), а не просто (1 2 1)
Спасибо!
Использованная литература:
Ответы 4
Для символьных вычислений
Наверное можно попробовать Ryacas пакет
library(Ryacas)
f <- function(n, ...) {
as_r(yac_str(sprintf("Simplify(Expand((%s)^%s))", paste0(c(...), collapse = "+"), n)))
}
и вы получите, например
> 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)
Для численных расчетов (только коэффициенты)
Если вас интересуют только числовые коэффициенты многочлена (например, многочлен от x), вы можете попробовать convolve + Reduce
h <- function(n, coeffs) {
Reduce(\(x, y) convolve(x, rev(y), type = "open"), rep(list(coeffs), n))
}
или его рекурсивный вариант
h2 <- function(n, coeffs) {
if (n == 1) {
return(coeffs)
}
convolve(Recall(n - 1, coeffs), rev(coeffs), type = "open")
}
и ты увидишь
# (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
Вы можете использовать giacR (еще не в CRAN):
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"
Использование пакета mpoly:
mpoly::mp("(a + b + c)^5")
#> 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
Или, для показателей и коэффициентов:
library(partitions)
f <- function(vars, pow) {
m <- as.matrix(compositions(pow, length(vars)))
rownames(m) <- vars
list(
exponents = m,
coeffs = exp(lgamma(pow + 1) - colSums(array(lgamma(1:(pow + 1))[m + 1L], dim(m))))
)
}
f(letters[1:3], 5)
#> $exponents
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
#> a 5 4 3 2 1 0 4 3 2 1 0 3 2 1 0 2 1 0 1 0 0
#> b 0 1 2 3 4 5 0 1 2 3 4 0 1 2 3 0 1 2 0 1 0
#> c 0 0 0 0 0 0 1 1 1 1 1 2 2 2 2 3 3 3 4 4 5
#>
#> $coeffs
#> [1] 1 5 10 10 5 1 5 20 30 20 5 10 30 30 10 10 20 10 5 5 1
Получение только показателей и коэффициентов происходит намного быстрее.
library(mpoly)
system.time(mp("(a + b + c + d + e)^20"))
#> user system elapsed
#> 11.78 0.30 12.07
system.time(f(letters[1:5], 20))
#> user system elapsed
#> 0.02 0.00 0.00
Если вы не хотите использовать внешние пакеты, а только базовый R, вы можете определить пользовательскую функцию, как показано ниже (но она не будет хорошо масштабироваться, когда у вас есть высокий порядок n или много терминов в ..., поскольку применяются expand.grid + apply)
g <- function(n, ...) {
comb <- expand.grid(rep(list(c(...)), n))
d <- table(
apply(comb, 1, \(x) {
v <- sort(table(x))
gsub("\\^1((?=\\D)|$)",
"",
paste0(paste(names(v), v, sep = "^"), collapse = "*"),
perl = TRUE
)
})
)
sort(d)
}
и вы получите, например
> g(5, "x", "y")
x^5 y^5 x*y^4 y*x^4 x^2*y^3 y^2*x^3
1 1 5 5 10 10
> g(3, "a", "b", "c")
a^3 b^3 c^3 a*b^2 a*c^2 b*a^2 b*c^2 c*a^2 c*b^2 a*b*c
1 1 1 3 3 3 3 3 3 6
Приведенный выше код просто представляет термины в таблице, и вы определенно можете выполнять дополнительные действия для создания выражения «сумма».