Session 4: the polynomial class

1. Session 4: the polynomial class Constructing and manipulating objects

2. Objectives • The class must have objects that behave as, and can be manipulated as, mathematical polynomials in a natural way • construction – structure is important • arithmetic operations (+, -, *, ^n OK, but "/")? • calculus – derivative and integral • display – printing and plotting • smooth morphing into functions when used as such • other convenience features • S3 class version only, at this stage. S Programming in R

3. Constructor function polynomial <- local({ ### basic constructor function - polynomials are functions .clip <- function(p) { ## ancillary p <- as.numeric(p) if((j <- length(p)) == 0) return(0) while(j > 1 && p[j] == 0) j <- j - 1 p[1:j] } S Programming in R

4. Constructor, cont'd function(coef = c(0,1), rat = all(coef %% 1 == 0)) { val <- local({ coef <- .clip(coef) .Coef <- rev(coef) rat <- rat function(x) { p <- rep(0, length(x)) for(a in .Coef) p <- a + x * p p } }) oldClass(val) <- "polynomial" val } }) S Programming in R

5. The "rat" flag • A display feature, only • "rat" flag is on, by default, if initially the coefficients are all whole numbers • Operations on 'rat' polynomials (including whole numbers) exclusively pass on the flag to any polynomials generated from them • A single non-rat operand turns off the rat flag (by default) • All arithmetic nevertheless done in floating point – rational form only applies to display. S Programming in R

6. Coercion to numeric > as.numeric function (x, ...) UseMethod("as.double") <environment: namespace:base> Hence as.double.polynomial <- function(x, ...) coef(x) converts to class "numeric". S Programming in R

7. Arithmetic operations • Uses methods for a "group generic" function, 'Ops' (which does not exist as a function as all) • Methods are functions of two variables, e1 and e2 • Dispatch to polynomial occurs if either has class "polynomial" • Several methods may be incorporated into one method function – another convenience • Take code in stages S Programming in R

8. Ops.polynomial <- function(e1, e2) { if(missing(e2)) return(switch(.Generic, "+" = e1, "-" = polynomial(-coef(e1), .getRat(e1)))) rat <- (if(class(e1) == "polynomial") .getRat(e1) else all(e1 %% 1 == 0)) && (if(class(e2) == "polynomial") .getRat(e2) else all(e2 %% 1 == 0)) e1 <- as.numeric(e1) e2 <- as.numeric(e2) S Programming in R

9. e1.op.e2 <- switch(.Generic, "+" = , "-" = { if((j <- length(e1) - length(e2)) < 0) e1 <- c(e1, rep(0, -j)) else if(j > 0) e2 <- c(e2, rep(0, j)) NextMethod(.Generic) }, "*" = .poly.mult(e1, e2), "/" = , "%/%" = .poly.quo.rem(e1, e2)$quotient, "%%" = .poly.quo.rem(e1, e2)$remainder, "^" = if(length(e2) > 1 || e2 < 0 || (e2 %% 1) != 0) { stop("positive integer powers only") } else { S Programming in R

10. switch(as.character(e2), "0" = 1, "1" = e1, { m <- e1 for(i in 2:e2) m <- .poly.mult(m, e1) m }) }, stop(.Generic, " is not supported for polynomials")) polynomial(e1.op.e2, rat = rat) } S Programming in R

11. A brief example > x <- polynomial() > p <- 1 - (1-x)^5 > p 5*x - 10*x^2 + 10*x^3 - 5*x^4 + x^5 > dp <- derivative(p) > dp 5 - 20*x + 30*x^2 - 20*x^3 + 5*x^4 > ip <- integral(p) > ip 5/2*x^2 - 10/3*x^3 + 5/2*x^4 - x^5 + 1/6*x^6 S Programming in R

12. Plotting > plot(p) > lines(derivative(p), col = "red") > tangents(p, x = 0:3/2, col = "blue") > tangents(derivative(p), x = 0:3/2, col = "blue") S Programming in R

13. S Programming in R

14. Polynomials and functions • Objects are functions and behave as such • To make an efficient version, use explicit conversion • May be evaluated at numerical arguments OR • May be evaluated at polynomial arguments • Solves two problems: numerical evaluation and change of origin/scale • Polynomial iteration is a by product, e.g. branching processes. S Programming in R

15. Example > p 5*x - 10*x^2 + 10*x^3 - 5*x^4 + x^5 > p(1:4) [1] 1 2 33 244 > p(-4:1) [1] -3124 -1023 -242 -31 0 1 > p(x-1) - 31 + 80*x - 80*x^2 + 40*x^3 - 10*x^4 + x^5 > p(x+1) 1 + x^5 S Programming in R

16. The function > as.function(p) function (x) { p <- -5 + x * 1 p <- 10 + x * p p <- -10 + x * p p <- 5 + x * p 0 + x * p } <environment: 019F634C> S Programming in R

17. How is this coercion done? as.function.polynomial <- function(x, ...) { b <- coef(x) p <- as.name("p") x <- as.name("x") b0 <- as.numeric(b) while ((an <- length(b0)) <= 1) b0 <- c(b, 0) statement <- call("{") statement[[i <- 2]] <- call("<-", p, call("+", b0[an-1], call("*", x, b[an]))) for (ai in rev(b0)[-(1:2)]) statement[[i <- i + 1]] <- call("<-", p, call("+", ai, call("*", x, p))) statement[[i]] <- statement[[i]][[3]] fun <- function(x) {} body(fun) <- statement fun } S Programming in R

18. Extras • Two classes of polynomial • polyInterp(x, y) – interpolation polynomial • polyZeros(z) – polynomial with specified zeros (roots) • Finding roots: • solve(p, b = 0) – roots of p(x) = b(x) • New generics • derivative(), integral(), "coef<-"() • New methods for old generics: • as.character, plot, lines, points, as.double, ... S Programming in R

19. Final notes • A personal exploratory example that has proved more useful (for real jobs) than expected (e.g. exploring polynomial regressions, Branching processes, &c) • Shows the power and flexibility of the object oriented features of the S language • With S4 classes there are some extra security features, (e.g. checking validity) but not very much • Extending to polynomials in more than one variable looks not to be entirely straightforward in an elegant way – recursively defined objects are not well supported. S Programming in R