| Title: | A Collection of Functions to Implement a Class for Univariate Polynomial Manipulations |
|---|---|
| Description: | A collection of functions to implement a class for univariate polynomial manipulations. |
| Authors: | Bill Venables [aut, cre] (S original), Kurt Hornik [aut] (R port), Martin Maechler [aut] (R port) |
| Maintainer: | Bill Venables <[email protected]> |
| License: | GPL-2 |
| Version: | 1.4-1 |
| Built: | 2024-10-09 06:38:41 UTC |
| Source: | CRAN |
Takes a polynomial argument and constructs an R function to evaluate it at arbitrary points.
## S3 method for class 'polynomial' as.function(x, ...)## S3 method for class 'polynomial' as.function(x, ...)
x |
An object of class |
... |
further arguments to be passed to or from methods. |
This is a method for the generic function as.function.
The polynomial is evaluated within the function using the Horner scheme.
Note that you can use the model-oriented predict method for
polynomials for purpose of evaluation (without explicit coercion to
a function), see the example below.
A function to evaluate the polynomial p.
as.function,
predict.polynomial
pr <- (poly.calc(-1:1) - 2 * polynomial(c(1, 2, 1)))^2 pr ## 4 + 20*x + 33*x^2 + 16*x^3 - 6*x^4 - 4*x^5 + x^6 prf <- as.function(pr) prf ## function (x) ## 4 + x * (20 + x * (33 + x * (16 + x * (-6 + x * (-4 + x * (1)))))) ## <environment: 0x402440f0> prf(-3:3) ## 1024 64 0 4 64 144 64 predict(pr, -3:3) ## 1024 64 0 4 64 144 64pr <- (poly.calc(-1:1) - 2 * polynomial(c(1, 2, 1)))^2 pr ## 4 + 20*x + 33*x^2 + 16*x^3 - 6*x^4 - 4*x^5 + x^6 prf <- as.function(pr) prf ## function (x) ## 4 + x * (20 + x * (33 + x * (16 + x * (-6 + x * (-4 + x * (1)))))) ## <environment: 0x402440f0> prf(-3:3) ## 1024 64 0 4 64 144 64 predict(pr, -3:3) ## 1024 64 0 4 64 144 64
Calculate the coefficients of a polynomial relative to a new origin on the x axis.
change.origin(p, o)change.origin(p, o)
p |
an object of class |
o |
a numeric scalar representing the new origin on the original scale. |
Let be a given polynomial and consider
writing . This function calculates
the coefficients and returns the result as a polynomial.
A polynomial with coefficients relative to the re-located x axis.
pr <- poly.calc(1:5) pr ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 change.origin(pr, 3) ## 4*x - 5*x^3 + x^5pr <- poly.calc(1:5) pr ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 change.origin(pr, 3) ## 4*x - 5*x^3 + x^5
Calculates the derivative of a univariate polynomial.
## S3 method for class 'polynomial' deriv(expr, ...)## S3 method for class 'polynomial' deriv(expr, ...)
expr |
an object of class |
... |
further arguments to be passed to or from methods. |
This is a method for the generic function deriv.
Derivative of the polynomial.
pr <- poly.calc(1:5) pr ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 deriv(pr) ## 274 - 450*x + 255*x^2 - 60*x^3 + 5*x^4pr <- poly.calc(1:5) pr ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 deriv(pr) ## 274 - 450*x + 255*x^2 - 60*x^3 + 5*x^4
Compute the greatest common divisor (GCD) and least common multiple (LCM) of a collection of polynomials and polylists.
## S3 method for class 'polylist' GCD(...) ## S3 method for class 'polynomial' GCD(...) ## S3 method for class 'polylist' LCM(...) ## S3 method for class 'polynomial' LCM(...)## S3 method for class 'polylist' GCD(...) ## S3 method for class 'polynomial' GCD(...) ## S3 method for class 'polylist' LCM(...) ## S3 method for class 'polynomial' LCM(...)
... |
a list of objects of class |
pl <- polylist(poly.from.roots(-1), poly.from.roots(c(-1, -1)), poly.from.roots(1)) GCD(pl) GCD(pl[-3]) LCM(pl) LCM(pl, pl, pl[[2]])pl <- polylist(poly.from.roots(-1), poly.from.roots(c(-1, -1)), poly.from.roots(1)) GCD(pl) GCD(pl[-3]) LCM(pl) LCM(pl, pl, pl[[2]])
Find the integral of a univariate polynomial.
## S3 method for class 'polynomial' integral(expr, limits = NULL, ...)## S3 method for class 'polynomial' integral(expr, limits = NULL, ...)
expr |
an object of class |
limits |
numeric vector of length 2 giving the integration limits. |
... |
further arguments to be passed to or from methods. |
If limits is not given, the integral of p from 0 to
‘x’'. Otherwise, the integral with the given integration
limits.
p <- poly.calc(1:5) p ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 deriv(p) ## 274 - 450*x + 255*x^2 - 60*x^3 + 5*x^4 integral(deriv(p)) - 120 ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5p <- poly.calc(1:5) p ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5 deriv(p) ## 274 - 450*x + 255*x^2 - 60*x^3 + 5*x^4 integral(deriv(p)) - 120 ## -120 + 274*x - 225*x^2 + 85*x^3 - 15*x^4 + x^5
Add a polynomial to an existing plot usually as a line plot.
## S3 method for class 'polynomial' lines(x, len = 1000, xlim = NULL, ylim = NULL, ...)## S3 method for class 'polynomial' lines(x, len = 1000, xlim = NULL, ylim = NULL, ...)
x |
an object of class |
len |
size of vector at which evaluations are to be made. |
xlim, ylim
|
the range of x and y values with sensible defaults. |
... |
additional arguments as for the |
This is a method for the generic function lines.
Lines representing the given polynomial are added to an existing plot. Values outside the current plot region are not shown.
lines,
points,
points.polynomial,
plot,
plot.polynomial.
plot (poly.calc(-1:5)) lines (poly.calc( 2:4), lty = 2) points(poly.calc(-2:6), pch = 4)plot (poly.calc(-1:5)) lines (poly.calc( 2:4), lty = 2) points(poly.calc(-2:6), pch = 4)
Group method for functions in the Math group.
## S3 method for class 'polynomial' Math(x, ...)## S3 method for class 'polynomial' Math(x, ...)
x |
an object of class |
... |
further arguments to be passed to or from methods, such
as |
Most math group functions are disallowed with polynomial arguments.
The only exceptions are ceiling, floor,
round, trunc, and signif
which may be used to transform the coefficients accordingly.
A polynomial with transformed coefficients.
Ops.polynomial,
Summary.polynomial.
op <- options(digits=18) p <- poly.from.values(1:4, (2:5)^2) ## 1 + 2.00000000000001*x + x^2 p <- round(p) ## 1 + 2*x + x^2 options(op)op <- options(digits=18) p <- poly.from.values(1:4, (2:5)^2) ## 1 + 2.00000000000001*x + x^2 p <- round(p) ## 1 + 2*x + x^2 options(op)
Convert a polynomial to monic form by dividing by the leading coefficient.
monic(p)monic(p)
p |
A polynomial. A warning is issued if the polynomial is identically zero. |
Similar in effect to p/as.numeric(p[length(p)]) but with some
safeguards against leading zero coefficients.
A polynomial proportional to p with leading coefficient 1.
Allows arithmetic operators to be used for polynomial calculations, such as addition, multiplication, division, etc.
## S3 method for class 'polynomial' Ops(e1, e2)## S3 method for class 'polynomial' Ops(e1, e2)
e1 |
an object of class |
e2 |
an object of class |
A polynomial got by performing the operation on the two arguments.
Math.polynomial,
Summary.polynomial.
p <- polynomial(c(1, 2, 1)) ## 1 + 2*x + x^2 r <- poly.calc(-1 : 1) ## -1*x + x^3 (r - 2 * p)^2 ## 4 + 20*x + 33*x^2 + 16*x^3 - 6*x^4 - 4*x^5 + x^6p <- polynomial(c(1, 2, 1)) ## 1 + 2*x + x^2 r <- poly.calc(-1 : 1) ## -1*x + x^3 (r - 2 * p)^2 ## 4 + 20*x + 33*x^2 + 16*x^3 - 6*x^4 - 4*x^5 + x^6
Plots polynomials, optionally allowing the “interesting” region to be automatically determined.
## S3 method for class 'polynomial' plot(x, xlim = 0:1, ylim = range(Px), type = "l", len = 1000, ..., xlab = "x", ylab = "P(x)")## S3 method for class 'polynomial' plot(x, xlim = 0:1, ylim = range(Px), type = "l", len = 1000, ..., xlab = "x", ylab = "P(x)")
x |
an object of class |
xlim |
the range to be encompassed by the x axis. |
ylim |
the range to be encompassed by the y axis. |
type |
as for |
len |
number of x points drawn. |
... |
additional arguments as for |
xlab, ylab
|
graphical parameters. |
This is a method for the generic function plot.
A plot of the polynomial is produced on the currently active device. Unless otherwise specified, the domain is chosen to enclose the real parts of all zeros, stationary points and zero itself.
plot,
lines,
points,
lines.polynomial,
points.polynomial.
plot(p <- poly.calc(-1:5))plot(p <- poly.calc(-1:5))
Add a polynomial to an existing plot usually as a point plot.
## S3 method for class 'polynomial' points(x, length = 100, ...)## S3 method for class 'polynomial' points(x, length = 100, ...)
x |
an object of class |
length |
size of x vector at which evaluations are to be made. |
... |
additional arguments as for the points generic. |
This is a method for the generic function points.
Points representing the given polynomial are added to an existing plot. Values outside the current plot region are not shown.
plot,
lines,
points,
plot.polynomial,
lines.polynomial.
plot(poly.calc(-1:5)) lines(poly.calc(2:4), lty=2) points(poly.calc(-2:6), pch=4)plot(poly.calc(-1:5)) lines(poly.calc(2:4), lty=2) points(poly.calc(-2:6), pch=4)
Calculate either the monic polynomial with specified zeros, or the Lagrange interpolation polynomial through the (x,y) points.
poly.calc(x, y, tol=sqrt(.Machine$double.eps), lab=dimnames(y)[[2]])poly.calc(x, y, tol=sqrt(.Machine$double.eps), lab=dimnames(y)[[2]])
x |
numeric vector specifying either the zeros of the desired polynomial if this is the only non-missing argument, or the x-values for Lagrange interpolation. |
y |
numeric vector or matrix specifying the y-values for the
Lagrange interpolation polynomial. If |
tol |
An absolute value tolerance, below which coefficients are treated as zero. |
lab |
If |
If y is a matrix, the result is a list of polynomials using
each column separately.
If x only is given, repeated zeros are allowed. If x
and y are given, repeated values in the x vector must
have identical y values associated with them (up to
tol), otherwise the first y-value only is used and a warning
is issued.
Either a polynomial object, or a list of polynomials, as appropriate.
In the latter case the object is of class "polylist".
poly.calc(rep(1,3)) ## -1 + 3*x - 3*x^2 + x^3 poly.calc(0:4, (0:4)^2 + 1) ## 1 + x^2 poly.calc(0:4, cbind(0:4, (0:4)^2 + 1), lab = letters[1:2]) ## List of polynomials: ## $a: ## x ## ## $b: ## 1 + x^2poly.calc(rep(1,3)) ## -1 + 3*x - 3*x^2 + x^3 poly.calc(0:4, (0:4)^2 + 1) ## 1 + x^2 poly.calc(0:4, cbind(0:4, (0:4)^2 + 1), lab = letters[1:2]) ## List of polynomials: ## $a: ## x ## ## $b: ## 1 + x^2
Construct the orthogonal polynomials on a given vector, up to a specified degree.
poly.orth(x, degree = length(unique(x)) - 1, norm = TRUE)poly.orth(x, degree = length(unique(x)) - 1, norm = TRUE)
x |
a numeric vector of abscissae. When evaluated at |
degree |
maximum degree required. The default is one fewer than
the number of distinct values in |
norm |
a logical indicating whether the polynomials should be normalized. |
A list of class "polylist" of objects of class
"polynomial" of degree 1, 2, ..., degree.
x <- rep(1:4, 1:4) # x with repetitions for weighting x ## [1] 1 2 2 3 3 3 4 4 4 4 polx <- poly.orth(x, 3) # calculate orthogonal polynomials polx ## List of polynomials: ## [[1]] ## 0.3162278 ## ## [[2]] ## -0.9486833 + 0.3162278*x ## ## [[3]] ## 2.139203 - 1.863177*x + 0.3450328*x^2 ## ## [[4]] ## -5.831564 + 8.80369*x - 3.803194*x^2 + 0.4930066*x^3 v <- sapply(polx, predict, x) # orthonormal basis round(crossprod(v), 10) # check orthonormalityx <- rep(1:4, 1:4) # x with repetitions for weighting x ## [1] 1 2 2 3 3 3 4 4 4 4 polx <- poly.orth(x, 3) # calculate orthogonal polynomials polx ## List of polynomials: ## [[1]] ## 0.3162278 ## ## [[2]] ## -0.9486833 + 0.3162278*x ## ## [[3]] ## 2.139203 - 1.863177*x + 0.3450328*x^2 ## ## [[4]] ## -5.831564 + 8.80369*x - 3.803194*x^2 + 0.4930066*x^3 v <- sapply(polx, predict, x) # orthonormal basis round(crossprod(v), 10) # check orthonormality
Create and manipulate lists of polynomials.
polylist(...) as.polylist(x) is.polylist(x)polylist(...) as.polylist(x) is.polylist(x)
... |
R objects, polynomials or capable of coercion to polynomials. |
x |
an R object. |
polylist takes a list of arguments, tries to convert each into
a polynomial (see polynomial), and sets the class of the
list to "polylist".
as.polylist tries to coerce its arguments to a polylist, and
will do so for arguments which are polynomials or lists thereof.
is.polylist tests whether its argument is a polylist.
This class has several useful methods, such as taking derivatives
(deriv) and antiderivatives (integral),
printing and plotting, subscripting, computing sums and products of
the elements, and methods for c, rep, and
unique.
## Calculate orthogonal polynomials pl <- poly.orth(rep(1:4, 1:4), 3) pl plot(pl) deriv(pl) integral(pl) sum(pl) prod(pl) unique(rep(pl, 3)[c(8, 12)])## Calculate orthogonal polynomials pl <- poly.orth(rep(1:4, 1:4), 3) pl plot(pl) deriv(pl) integral(pl) sum(pl) prod(pl) unique(rep(pl, 3)[c(8, 12)])
Construct, coerce to, test for, and print polynomial objects.
polynomial(coef = c(0, 1)) as.polynomial(p) is.polynomial(p) ## S3 method for class 'polynomial' as.character(x, decreasing = FALSE, ...) ## S3 method for class 'polynomial' print(x, digits = getOption("digits"), decreasing = FALSE, ...)polynomial(coef = c(0, 1)) as.polynomial(p) is.polynomial(p) ## S3 method for class 'polynomial' as.character(x, decreasing = FALSE, ...) ## S3 method for class 'polynomial' print(x, digits = getOption("digits"), decreasing = FALSE, ...)
coef |
numeric vector, giving the polynomial coefficients in increasing order. |
p |
an arbitrary R object. |
x |
a |
decreasing |
a logical specifying the order of the terms; in increasing (default) or decreasing powers. |
digits |
the number of significant digits to use for printing. |
... |
potentially further arguments passed to and from other methods. |
polynomial constructs a polynomial from its coefficients,
i.e., p[1:k] specifies the polynomial
Internally, polynomials are simply numeric coefficient vectors of
class "polynomial". Several useful methods are available for
this class, such as coercion to character (as.character()) and
function (as.function.polynomial), extraction of
the coefficients (coef()), printing (using as.character),
plotting (plot.polynomial), and computing sums and
products of arbitrarily many polynomials.
as.polynomial tries to coerce its arguments to a polynomial.
is.polynomial tests whether its argument is a polynomial (in
the sense that it has class "polynomial".
polynomial(1:4) p <- as.polynomial(c(1,0,3,0)) p print(p, decreasing = TRUE) stopifnot(coef(p) == c(1,0,3)) polynomial(c(2,rep(0,10),1))polynomial(1:4) p <- as.polynomial(c(1,0,3,0)) p print(p, decreasing = TRUE) stopifnot(coef(p) == c(1,0,3)) polynomial(c(2,rep(0,10),1))
Evaluate a polynomial at a given numeric or polynomial argument.
## S3 method for class 'polynomial' predict(object, newdata, ...)## S3 method for class 'polynomial' predict(object, newdata, ...)
object |
A polynomial object to be evaluated. |
newdata |
Argument at which evaluation is requested. May be numeric or itself a polynomial |
... |
Not used by this method. |
This is a method for the generic function predict.
The polynomial is evaluated according to the Horner scheme for speed and numerical accuracy.
Evaluated object of the same class as newdata.
Find the zeros, if any, of a given polynomial.
## S3 method for class 'polynomial' solve(a, b, ...)## S3 method for class 'polynomial' solve(a, b, ...)
a |
A polynomial object for which the zeros are required. |
b |
a numeric value specifying an additional intercept. If
given, the zeros of |
... |
Not used by this method. |
This is a method for the generic function solve.
The zeros are found as the eigenvalues of the companion matrix, sorted according to their real parts.
A numeric vector, generally complex, of zeros.
polyroot,
poly.calc,
summary.polynomial
p <- polynomial(6:1) p ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5 pz <- solve(p) pz ## [1] -1.49180+0.0000i -0.80579-1.2229i -0.80579+1.2229i ## [4] 0.55169-1.2533i 0.55169+1.2533i ## To retrieve the original polynomial from the zeros: poly.calc(pz) ## Warning: imaginary parts discarded in coercion ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5p <- polynomial(6:1) p ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5 pz <- solve(p) pz ## [1] -1.49180+0.0000i -0.80579-1.2229i -0.80579+1.2229i ## [4] 0.55169-1.2533i 0.55169+1.2533i ## To retrieve the original polynomial from the zeros: poly.calc(pz) ## Warning: imaginary parts discarded in coercion ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5
Summarize a polynomial by describing its “key” points.
## S3 method for class 'polynomial' summary(object, ...)## S3 method for class 'polynomial' summary(object, ...)
object |
an object of class |
... |
Not used by this method. |
This is a method for the generic function summary.
A list of class "summary.polynomial" (which has its own
print method) containing information on zeros, stationary and
inflexion points.
p <- polynomial(6:1) p ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5 pz <- summary(p) pz ## [1] -1.49180+0.0000i -0.80579-1.2229i -0.80579+1.2229i ## [4] 0.55169-1.2533i 0.55169+1.2533i ## To retrieve the original polynomial from the zeros: poly.calc(pz) ## Warning: imaginary parts discarded in coercion ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5p <- polynomial(6:1) p ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5 pz <- summary(p) pz ## [1] -1.49180+0.0000i -0.80579-1.2229i -0.80579+1.2229i ## [4] 0.55169-1.2533i 0.55169+1.2533i ## To retrieve the original polynomial from the zeros: poly.calc(pz) ## Warning: imaginary parts discarded in coercion ## 6 + 5*x + 4*x^2 + 3*x^3 + 2*x^4 + x^5
Allows summary group generics to be used on polynomial arguments.
## S3 method for class 'polynomial' Summary(..., na.rm = FALSE)## S3 method for class 'polynomial' Summary(..., na.rm = FALSE)
... |
R objects, the first supplied of class
|
na.rm |
logical: should missing values be removed? |
For the sum and prod functions, the sum and product of
the given polynomials, respectively. For the other members of the
Summary group, an error is returned.
Note that one could order polynomials by divisibility, and
define min and max as the corresponding lattice meet and
join, i.e., the greatest common divisor and the least common multiple,
respectively. This is currently not provided: instead, functions
GCD and LCM should be called directly.
Math.polynomial,
Ops.polynomial