# 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)
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
# 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)
#
# 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)``````