mpoly

Specifying polynomials

mpoly is a simple collection of tools to help deal with multivariate polynomials symbolically and functionally in R. Polynomials are defined with the mp() function:

library(mpoly)
# Loading required package: stringr
mp("x + y")
# x  +  y
mp("(x + 4y)^2 (x - .25)")
# x^3  -  0.25 x^2  +  8 x^2 y  -  2 x y  +  16 x y^2  -  4 y^2

Term orders are available with the reorder function:

(p <- mp("(x + y)^2 (1 + x)"))
# x^3  +  x^2  +  2 x^2 y  +  2 x y  +  x y^2  +  y^2
reorder(p, varorder = c("y","x"), order = "lex")
# y^2 x  +  y^2  +  2 y x^2  +  2 y x  +  x^3  +  x^2
reorder(p, varorder = c("x","y"), order = "glex")
# x^3  +  2 x^2 y  +  x y^2  +  x^2  +  2 x y  +  y^2

Vectors of polynomials (mpolyList’s) can be specified in the same way:

mp(c("(x+y)^2", "z"))
# x^2  +  2 x y  +  y^2
# z

Polynomial parts

You can extract pieces of polynoimals using the standard [ operator, which works on its terms:

p[1]
# x^3
p[1:3]
# x^3  +  x^2  +  2 x^2 y
p[-1]
# x^2  +  2 x^2 y  +  2 x y  +  x y^2  +  y^2

There are also many other functions that can be used to piece apart polynomials, for example the leading term (default lex order):

LT(p)
# x^3
LC(p)
# [1] 1
LM(p)
# x^3

You can also extract information about exponents:

exponents(p)
# [[1]]
# x y 
# 3 0 
# 
# [[2]]
# x y 
# 2 0 
# 
# [[3]]
# x y 
# 2 1 
# 
# [[4]]
# x y 
# 1 1 
# 
# [[5]]
# x y 
# 1 2 
# 
# [[6]]
# x y 
# 0 2
multideg(p)
# x y 
# 3 0
totaldeg(p)
# [1] 3
monomials(p)
# x^3
# x^2
# 2 x^2 y
# 2 x y
# x y^2
# y^2

Polynomial arithmetic

Arithmetic is defined for both polynomials (+, -, * and ^)…

p1 <- mp("x + y")
p2 <- mp("x - y")

p1 + p2
# 2 x
p1 - p2
# 2 y
p1 * p2
# x^2  -  y^2
p1^2
# x^2  +  2 x y  +  y^2

… and vectors of polynomials:

(ps1 <- mp(c("x", "y")))
# x
# y
(ps2 <- mp(c("2x", "y+z")))
# 2 x
# y  +  z
ps1 + ps2
# 3 x
# 2 y  +  z
ps1 - ps2
# -1 x
# -1 z
ps1 * ps2
# 2 x^2
# y^2  +  y z

Some calculus

You can compute derivatives easily:

p <- mp("x + x y + x y^2")
deriv(p, "y")
# x  +  2 x y
gradient(p)
# y^2  +  y  +  1
# 2 y x  +  x

Function coercion

You can turn polynomials and vectors of polynomials into functions you can evaluate with as.function(). Here’s a basic example using a single multivariate polynomial:

f <- as.function(mp("x + 2 y")) # makes a function with a vector argument
# f(.) with . = (x, y)
f(c(1,1))
# [1] 3
f <- as.function(mp("x + 2 y"), vector = FALSE) # makes a function with all arguments
# f(x, y)
f(1, 1)
# [1] 3

Here’s a basic example with a vector of multivariate polynomials:

(p <- mp(c("x", "2 y")))
# x
# 2 y
f <- as.function(p) 
# f(.) with . = (x, y)
f(c(1,1))
# [1] 1 2
f <- as.function(p, vector = FALSE) 
# f(x, y)
f(1, 1)
# [1] 1 2

Whether you’re working with a single multivariate polynomial or a vector of them (mpolyList), if it/they are actually univariate polynomial(s), the resulting function is vectorized. Here’s an example with a single univariate polynomial.

f <- as.function(mp("x^2"))
# f(.) with . = x
f(1:3)
# [1] 1 4 9
(mat <- matrix(1:4, 2))
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4
f(mat) # it's vectorized properly over arrays
#      [,1] [,2]
# [1,]    1    9
# [2,]    4   16

Here’s an example with a vector of univariate polynomials:

(p <- mp(c("t", "t^2")))
# t
# t^2
f <- as.function(p)
# f(.) with . = (t)
f(1)
# [1] 1 1
f(1:3)
#      [,1] [,2]
# [1,]    1    1
# [2,]    2    4
# [3,]    3    9

You can use this to visualize a univariate polynomials like this:

f <- as.function(mp("(x-2) x (x+2)"))
# f(.) with . = x
x <- seq(-2.5, 2.5, .1)

library(ggplot2); theme_set(theme_classic())
qplot(x, f(x), geom = "line")

For multivariate polynomials, it’s a little more complicated:

f <- as.function(mp("x^2 - y^2")) 
# f(.) with . = (x, y)
s <- seq(-2.5, 2.5, .1)
df <- expand.grid(x = s, y = s)
df$f <- apply(df, 1, f)
qplot(x, y, data = df, geom = "tile", fill = f)

Algebraic geometry

Grobner bases are no longer implemented, see m2r

# polys <- mp(c("t^4 - x", "t^3 - y", "t^2 - z"))
# grobner(polys)

Homogenization and dehomogenization:

(p <- mp("x + 2 x y + y - z^3"))
# x  +  2 x y  +  y  -  z^3
(hp <- homogenize(p))
# x t^2  +  2 x y t  +  y t^2  -  z^3
dehomogenize(hp, "t")
# x  +  2 x y  +  y  -  z^3
homogeneous_components(p)
# x  +  y
# 2 x y
# -1 z^3

Special polynomials

mpoly can make use of many pieces of the polynom and orthopolynom packages with as.mpoly() methods. In particular, many special polynomials are available.

Chebyshev polynomials

You can construct Chebyshev polynomials as follows:

chebyshev(1)
# Loading required package: polynom
# 
# Attaching package: 'polynom'
# The following object is masked from 'package:mpoly':
# 
#     LCM
# x
chebyshev(2)
# -1  +  2 x^2
chebyshev(0:5)
# 1
# x
# 2 x^2  -  1
# 4 x^3  -  3 x
# 8 x^4  -  8 x^2  +  1
# 16 x^5  -  20 x^3  +  5 x

And you can visualize them:

library(tidyr); library(dplyr)
s <- seq(-1, 1, length.out = 201); N <- 5
(chebPolys <- chebyshev(0:N))
# 1
# x
# 2 x^2  -  1
# 4 x^3  -  3 x
# 8 x^4  -  8 x^2  +  1
# 16 x^5  -  20 x^3  +  5 x


df <- as.function(chebPolys)(s) %>% cbind(s, .) %>% as.data.frame
# f(.) with . = (x)
names(df) <- c("x", paste0("T_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)

Jacobi polynomials

s <- seq(-1, 1, length.out = 201); N <- 5
(jacPolys <- jacobi(0:N, 2, 2))
# 1
# 5 x
# 17.5 x^2  -  2.5
# 52.5 x^3  -  17.5 x
# 144.375 x^4  -  78.75 x^2  +  4.375
# 375.375 x^5  -  288.75 x^3  +  39.375 x
 
df <- as.function(jacPolys)(s) %>% cbind(s, .) %>% as.data.frame
# f(.) with . = (x)
names(df) <- c("x", paste0("P_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree) +
  coord_cartesian(ylim = c(-25, 25))

Legendre polynomials

s <- seq(-1, 1, length.out = 201); N <- 5
(legPolys <- legendre(0:N))
# 1
# x
# 1.5 x^2  -  0.5
# 2.5 x^3  -  1.5 x
# 4.375 x^4  -  3.75 x^2  +  0.375
# 7.875 x^5  -  8.75 x^3  +  1.875 x
 
df <- as.function(legPolys)(s) %>% cbind(s, .) %>% as.data.frame
# f(.) with . = (x)
names(df) <- c("x", paste0("P_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)

Hermite polynomials

s <- seq(-3, 3, length.out = 201); N <- 5
(hermPolys <- hermite(0:N))
# 1
# x
# x^2  -  1
# x^3  -  3 x
# x^4  -  6 x^2  +  3
# x^5  -  10 x^3  +  15 x

df <- as.function(hermPolys)(s) %>% cbind(s, .) %>% as.data.frame
# f(.) with . = (x)
names(df) <- c("x", paste0("He_", 0:N))
mdf <- df %>% gather(degree, value, -x)
qplot(x, value, data = mdf, geom = "path", color = degree)