Evaluate Polynomials

Evaluate one or several univariate polynomials at several locations, i.e. compute coef[1] + coef[2]*x + ... + coef[p+1]* x^p (in the simplest case where x is scalar and coef a vector).

Usage

polyn.eval(coef, x)

Arguments

coef

“numeric” vector or matrix. If a vector, x can be an array and the result matches x.
If coef is a matrix it specifies several polynomials of the same degree as rows, x must be a vector, coef[,k] is for \(x^{k-1}\) and the result is a matrix of dimension length(x) * nrow(coef).

Note that coef can also be complex or bigrational (as.bigq(.) from gmp, or arbitrary precision ("mpfr") from Rmpfr, or similar number-like objects for which basic arithmetic is defined.

x

“numeric” vector or array. Either x or coef must be a vector.

Details

The stable “Horner rule” is used for evaluation in any case.

When length(coef) == 1L, polyn.eval(coef, x) now returns a vector of length(x) whereas previously, it just gave the number coef independent of x.

Value

numeric vector or array, depending on input dimensionalities, see above.

Author

Martin Maechler, ages ago.

See also

For much more sophisticated handling of polynomials, use the polynom package, see, e.g., predict.polynomial. For multivariate polynomials (and also for nice interface to the orthopolynom package), consider the mpoly package.

Examples

polyn.eval(c(1,-2,1), x = 0:3)# (x - 1)^2
#> [1] 1 0 1 4
polyn.eval(c(0, 24, -50, 35, -10, 1), x = matrix(0:5, 2,3))# 5 zeros!
#>      [,1] [,2] [,3]
#> [1,]    0    0    0
#> [2,]    0    0  120
(cf <- rbind(diag(3), c(1,-2,1)))
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    1    0
#> [3,]    0    0    1
#> [4,]    1   -2    1
polyn.eval(cf, 0:5)
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    1
#> [2,]    1    1    1    0
#> [3,]    1    2    4    1
#> [4,]    1    3    9    4
#> [5,]    1    4   16    9
#> [6,]    1    5   25   16

x  <- seq(-3,7, by=1/4)
cf <- 4:1
(px <- polyn.eval(cf, x)) # is exact
#>  [1] -14.000  -9.922  -6.625  -4.016  -2.000  -0.484   0.625   1.422   2.000
#> [10]   2.453   2.875   3.359   4.000   4.891   6.125   7.797  10.000  12.828
#> [19]  16.375  20.734  26.000  32.266  39.625  48.172  58.000  69.203  81.875
#> [28]  96.109 112.000 129.641 149.125 170.547 194.000 219.578 247.375 277.484
#> [37] 310.000 345.016 382.625 422.922 466.000
if((gmpT <-"package:gmp" %in% search()) || require("gmp")) withAutoprint({
 pxq <- polyn.eval(coef = as.bigq(cf, 1), x=x)
 pxq
 stopifnot(pxq == px)
 if(!gmpT) detach("package:gmp")
})
#> Loading required package: gmp
#> 
#> Attaching package: ‘gmp’
#> The following objects are masked from ‘package:sfsmisc’:
#> 
#>     factorize, formatN, is.whole
#> The following objects are masked from ‘package:base’:
#> 
#>     %*%, apply, crossprod, matrix, tcrossprod
#> > pxq <- polyn.eval(coef = as.bigq(cf, 1), x = x)
#> > pxq
#> Big Rational ('bigq') object of length 41:
#>  [1] -14      -635/64  -53/8    -257/64  -2       -31/64   5/8      91/64   
#>  [9] 2        157/64   23/8     215/64   4        313/64   49/8     499/64  
#> [17] 10       821/64   131/8    1327/64  26       2065/64  317/8    3083/64 
#> [25] 58       4429/64  655/8    6151/64  112      8297/64  1193/8   10915/64
#> [33] 194      14053/64 1979/8   17759/64 310      22081/64 3061/8   27067/64
#> [41] 466     
#> > stopifnot(pxq == px)
#> > if (!gmpT) detach("package:gmp")

if((RmpfrT <-"package:Rmpfr" %in% search()) || require("Rmpfr")) withAutoprint({
 pxM <- polyn.eval(coef = mpfr(cf, 80), x=x) # 80 bits accuracy
 pxM
 stopifnot(pxM == px)
 if(!RmpfrT) detach("package:Rmpfr")
})
#> Loading required package: Rmpfr
#> Loading required package: gmp
#> 
#> Attaching package: ‘gmp’
#> The following objects are masked from ‘package:sfsmisc’:
#> 
#>     factorize, formatN, is.whole
#> The following objects are masked from ‘package:base’:
#> 
#>     %*%, apply, crossprod, matrix, tcrossprod
#> C code of R package 'Rmpfr': GMP using 64 bits per limb
#> 
#> Attaching package: ‘Rmpfr’
#> The following object is masked from ‘package:gmp’:
#> 
#>     outer
#> The following objects are masked from ‘package:stats’:
#> 
#>     dbinom, dchisq, dgamma, dnbinom, dnorm, dpois, dt, pgamma, pnorm
#> The following objects are masked from ‘package:base’:
#> 
#>     cbind, pmax, pmin, rbind
#> > pxM <- polyn.eval(coef = mpfr(cf, 80), x = x)
#> > pxM
#> 41 'mpfr' numbers of precision  80   bits 
#>  [1]        -14  -9.921875     -6.625  -4.015625         -2  -0.484375
#>  [7]      0.625   1.421875          2   2.453125      2.875   3.359375
#> [13]          4   4.890625      6.125   7.796875         10  12.828125
#> [19]     16.375  20.734375         26  32.265625     39.625  48.171875
#> [25]         58  69.203125     81.875  96.109375        112 129.640625
#> [31]    149.125 170.546875        194 219.578125    247.375 277.484375
#> [37]        310 345.015625    382.625 422.921875        466
#> > stopifnot(pxM == px)
#> > if (!RmpfrT) detach("package:Rmpfr")

stopifnot(identical(polyn.eval(12, x),      rep(12, length(x))),
          identical(polyn.eval(7, diag(3)), matrix(7, 3,3)))