17

In R how would one represent polynomial expressions and do polynomial math with the numeric vector objects? For example:

x1 <- c(2,1)  # 2 + x
x2 <- c(-1,3)  # -1 + 3*x

And want:

x1 * x2 # to return -2 + 5*x + 3*x^2

Note: I answered a question this morning and then the poster apparently deleted it (making me wonder if it were homework.) So I am re-posting the question from memory.

animuson
  • 53,861
  • 28
  • 137
  • 147
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • 3
    Try the mpoly package? – hadley Jun 03 '13 at 13:00
  • That package does look more full-featured than 'polynom' in that it handles multivariate operations. I'm not intending to checkmark my own answer, so I invite better answers than the one I provided. – IRTFM Jun 03 '13 at 17:03

3 Answers3

16

One could multiply the coefficients directly using outer and then aggregate the results

x1 <- c(2,1)  # 2 + x
x2 <- c(-1,3)  # -1 + 3*x
tmp <- outer(x1, x2)
tapply(tmp, row(tmp) + col(tmp) - 1, sum)
# 1  2  3 
#-2  5  3

x1 <- c(2, 1) # 2 + x
x2 <- c(-1, 3, 2) # -1 + 3*x + 2*x^2
tmp <- outer(x1, x2)
tapply(tmp, row(tmp) + col(tmp) - 1, sum) # should give -2 + 5*x + 7*x^2 + 2*x^3
# 1  2  3  4 
#-2  5  7  2

as discussed in the comments the '-1' in the code isn't necessary. When coming up with the solution that helped me because it allowed me to map each location in the output of outer to where it would end up in the final vector. If we did a '-2' instead then it would map to the exponent on x in the resulting polynomial. But we really don't need it so something like the following would work just as well:

tmp <- outer(x1, x2)
tapply(tmp, row(tmp) + col(tmp), sum)
Dason
  • 60,663
  • 9
  • 131
  • 148
  • That's kewl. Kind of a binomial expansion triangle tilted 45 degrees counterclockwise. – IRTFM May 31 '13 at 20:06
  • This is very similar to the code in `polynom:::Ops.polynomial`, although they use `tapply(m, row(m) + col(m), sum)`... Why do both give the same result? – Ferdinand.kraft May 31 '13 at 20:20
  • I actually don't need the `-1 in the code but it helped me mentally think about the problem at first. I should probably just get rid of it. Those numbers are just being used to group together results so we could add/subtract any amount and it would stay the same. – Dason May 31 '13 at 20:41
  • If you subtracted 2 instead of 1 you could use the sum as the exponent for `x`. – IRTFM May 31 '13 at 21:52
  • Yeah - I was using it as the index of the resulting vector (thus only subtracting 1 instead of 2) - but either way it's unnecessary. – Dason May 31 '13 at 22:25
13

Use the polynom package:

 require(polynom)
# Loading required package: polynom
# From the example for as.polynomial
 p <- as.polynomial(c(1,0,3,0))
 p
# 1 + 3*x^2 

 x1 <- c(2,1)
 x2 <- c(-1,3)
 px1 <- as.polynomial(x1)
 px2 <- as.polynomial(x2)

 px1*px2
# -2 + 5*x + 3*x^2 
 prod.p <- .Last.value
 str(prod.p)
# Class 'polynomial'  num [1:3] -2 5 3
 unclass(prod.p)
# [1] -2  5  3
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • 2
    May I add that this does not allow symbolic manipulation of polynomial expressions, since the coefficients are stored as regular `double` values and precision is lost during arithmetic operations. If one needs to perform symbolic computations on such expressions, I'd suggest Mathematica instead. – Ferdinand.kraft May 31 '13 at 19:53
  • 6
    There is also an Ryacas package. – IRTFM May 31 '13 at 19:57
  • If R cannot do what is wanted, I think Maxima does polynomials and is free. Although www.wolframalpha.com is free too and I think is based on Mathematica. – Mark Miller May 31 '13 at 20:05
  • R definitely can do what David asked: and the answer *is* using the polynom package, as long as you do not need exact arithmetic. – Martin Mächler Jun 15 '13 at 08:42
  • 1
    For exact (or arbitrary precise) arithmetic, I'd also like to mention the `gmp` R package (based on the GNU MP ("GMP") library), for exact integer and rational arithmetic, *and* (my) package `Rmpfr`, based on the (GNU) MPFR library for arbitrary precision arithmetic. What really would be neat is a marriage of 'polynom' with gmp and Rmpfr. I'd be happy to collaborate on that. – Martin Mächler Jun 15 '13 at 08:44
2

Polynomial multiplication is the convolution of the coefficients

convolve(c(2,1),rev(c(-1,3)),type="open")
#[1] -2  5  3
A. Webb
  • 26,227
  • 1
  • 63
  • 95