knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 3.5, message = FALSE, warning = FALSE) options(width = 80, tibble.width = Inf)

The `mvp`

package provides some functionality for fast
manipulation of multivariate polynomials, using the Standard Template
library of `C++`

, commonly known as the `STL`

. It is
comparable in speed to the `spray`

package for sparse arrays,
while retaining the symbolic capabilities of the `mpoly`

package
[@kahle2013]. The `mvp package`

uses the excellent print and
coercion methods of `mpoly`

. The `mvp`

package provides
improved speed over `mpoly`

, the ability to handle negative
powers, and a more sophisticated substitution mechanism.

`STL map`

classA `map`

is a sorted associative container that contains key-value
pairs with unique keys. It is interesting here because search and
insertion operations have logarithmic complexity. Multivariate
polynomials are considered to be the sum of a finite number of
*terms*, each multiplied by a coefficient. A *term* is something like
$x^2y^3z$. We may consider this term to be the map

{"x" -> 2, "y" -> 3, "z" -> 1}

where the map takes symbols to their (integer) power; it is understood
that powers are nonzero. A `mvp`

object is a map from terms to
their coeffients; thus $7xy^2 -3x^2yz^5$ would be

{{"x" -> 1, "y" -> 2} -> 7, {"x" -> 2, 'y" -> 1, "z" ->5} -> -3}

and we understand that coefficients are nonzero. In `C++`

the
declarations would be

typedef vector <signed int> mypowers; typedef vector <string> mynames; typedef map <string, signed int> term; typedef map <term, double> mvp;

Thus a `term`

maps a string to a (signed) integer, and a `mvp`

maps terms to doubles. One reason why the `map`

class is fast is
that the order in which the keys are stored is undefined: the compiler
may store them in the order which it regards as most propitious. This
is not an issue for the maps considered here as addition and
multiplication are commutative and associative.

Note also that constant terms are handled with no difficulty (constants are simply maps from the empty map to its value), as is the zero polynomial (which is simply an empty map).

Consider a simple multivariate polynomial $3xy+z^3+xy^6z$ and its representation in the following R session:

library("mvp",quietly=TRUE) (p <- as.mvp("3 x y + z^3 + x y^6 z"))

Coercion and printing are accomplished by the `mpoly`

package
(there is no way I could improve upon Kahle's work). Note carefully
that the printed representation of the mvp object is created by the
`mpoly`

package and the print method can rearrange both the terms
of the polynomial ($3xy+z^3+xy^6z = z^3+3xy+xy^6z$, for example) and
the symbols within a term ($3xy=3yx$, for example) to display the
polynomial in a human-friendly form.

However, note carefully that such rearranging does not affect the
mathematical properties of the polynomial itself. In the `mvp`

package, the order of the terms is not preserved (or even defined) in
the internal representation of the object; and neither is the order of
the symbols within a single term. Although this might sound odd, if
we consider a marginally more involved situation, such as

(M <- as.mvp("3 stoat goat^6 -4 + 7 stoatboat^3 bloat -9 float boat goat gloat^6")) dput(M)

it is not clear that any human-discernable ordering is preferable to
any other, and we would be better off letting the compiler decide a
propitious ordering. In any event, the `mpoly`

package can
specify a print order:

print(M,order="lex", varorder=c("stoat","goat","boat","bloat","gloat","float","stoatboat"))

The arithmetic operations `*`

, `+`

, `-`

and `^`

work
as expected:

(S1 <- rmvp(5,2,2,4)) (S2 <- rmvp(5,2,2,4)) S1 + S2 S1 * S2 S1^2

The package has two substitution functionalities. Firstly, we can substitute one or more variables for a numeric value. Define a mvp object:

(S3 <- as.mvp("x + 5 x^4 y + 8 y^2 x z^3"))

And then we may substitute $x=1$:

subs(S3, x = 1)

Note the natural R idiom, and that the return value is another mvp object. We may subsitute for the other variables:

subs(S3, x = 1, y = 2, z = 3)

(in this case, the default behaviour is to return a "dropped" version
of the resulting polynomial, that is, coerced to a scalar). The idiom
also allows one to substitute a variable for an `mvp`

object:

subs(as.mvp("a+b+c"), a="x^6")

Note carefully that `subs()`

depends on the order of substitution:

subs(as.mvp("a+b+c"), a="x^6",x="1+a") subs(as.mvp("a+b+c"), x="1+a",a="x^6")

Substitution works well with pipes:

as.mvp("a+b") %>% subs(a="a^2+b^2") %>% subs(b="x^6")

Differentiation is implemented. First we have the `deriv()`

method:

(S <- as.mvp("a + 5 a^5*b^2*c^8 -3 x^2 a^3 b c^3")) deriv(S, letters[1:3]) deriv(S, rev(letters[1:3])) # should be the same.

Also a slightly different form: `aderiv()`

, here used to evaluate
$\frac{\partial^6S}{\partial a^3\partial b\partial c^2}$:

aderiv(S, a = 3, b = 1, c = 2)

Again, pipes work quite nicely:

S %<>% aderiv(a=1,b=2) %>% subs(c="x^4") %>% `+`(as.mvp("o^99")) S

The package includes functionality to deal with Taylor and Laurent series:

(X <- as.mvp("1+x+x^2 y")^3) taylor(X,3) # terms with total power <= 3 taylor(X,3,"x") # terms with power of x <= 3 onevarpow(X,3,"x") # terms with power of x == 3

## second order taylor expansion of f(x)=sin(x+y) for x=1.1, about x=1: sinxpy <- horner("x+y",c(0,1,0,-1/6,0,+1/120,0,-1/5040)) # sin(x+y) dx <- as.mvp("dx") t2 <- sinxpy + aderiv(sinxpy,x=1)*dx + aderiv(sinxpy,x=2)*dx^2/2 (t2 %<>% subs(x=1,dx=0.1)) # (Taylor expansion of sin(y+1.1), left in symbolic form) (t2 %>% subs(y=0.3)) - sin(1.4) # numeric; should be small

The `mvp`

package handles negative powers, although the idiom is not perfect and I'm still working on it.
There is the `invert()`

function:

(p <- as.mvp("1+x+x^2 y")) invert(p)

In the above, `p`

is a regular multivariate polynomial which
includes negative powers. It obeys the same arithmetic rules as other
mvp objects:

p + as.mvp("z^6")

We can see the generating function for a chess knight:

```
knight(2)
```

How many ways are there for a 4D knight to return to its starting square after four moves? Answer:

constant(knight(4)^4)

I will show some timings using a particularly favourable example that
exploits the symbolic nature of the `mvp`

package.

library("spray") library("mpoly") n <- 100 k <- kahle(n, r = 3, p = 1:3, symbols = paste0("x", sprintf("%03d", 1:n)))

In the above, polynomial `k`

has 500 terms of the form $xy^2z^3$.
Coercing `k`

to `spray`

form would need a $500\times 500$
matrix for the indices, almost every element of which would be zero.
This makes the **spray** package slower:

library("microbenchmark") spray_k <- mvp_to_spray(k) microbenchmark( k^2, spray_k^2 )

In the above, the first line uses `mvp`

functionality, and the
second line uses `spray`

functionality. The speedup increases for
larger polynomials.

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.