polyn.eval.RdEvaluate 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).
polyn.eval(coef, x)“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.
“numeric” vector or array. Either x or coef must
be a vector.
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.
numeric vector or array, depending on input dimensionalities, see above.
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.
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, 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, 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)))