R/polybase.R

In eestileib/ComplexPoly: A Collection of Functions to Implement a Class for Univariate Polynomial Manipulations

polynomial <-
function(coef = c(0, 1))
{
    a <- as.complex(coef)
    while((la <- length(a)) > 1 && a[la] == 0) a <- a[-la]
    structure(a, class = "polynomial")
}

as.polynomial <-
function(p)
    if(is.polynomial(p)) p else polynomial(p)
is.polynomial <-
function(p)
    inherits(p, "polynomial")

Ops.polynomial <-
function(e1, e2)
{
    if(missing(e2))
        return(switch(.Generic,
                      "+" = e1,
                      "-" = polynomial(NextMethod(.Generic)),
                      stop("unsupported unary operation")))
    e1 <- unclass(e1)
    e2 <- unclass(e2)
    l1 <- length(e1)
    l2 <- length(e2)
    e1.op.e2 <-
        switch(.Generic,
               "+" = ,
               "-" = {
                   e1 <- c(e1, rep.int(0, max(0, l2 - l1)))
                   e2 <- c(e2, rep.int(0, max(0, l1 - l2)))
                   NextMethod(.Generic)
               },
               "*" = if(l1 == 1 || l2 == 1) e1 * e2 else {
                   m <- outer(e1, e2)
                   as.vector(tapply(m, row(m) + col(m), sum))
               },
               "/" = {
                   if(l2 == 0)
                       stop("unsupported polynomial division")
                   if(l2 == 1)
                       e1 / e2
                   else {
                       p <- rev(e1)
                       q <- rev(e2)
                       r <- rep.int(0, length(p))
                       i <- 0
                       while(length(p) >= l2) {
                           i <- i + 1
                           d <- p[1]/q[1]
                           r[i] <- d
                           p[1:l2] <- p[1:l2] - d * q
                           p <- p[-1]
                       }
                       if(i == 0) 0 else r[i:1]
                   }
               },
               "^" = {
                   if(l2 != 1 || e2 < 0 || e2 %% 1 != 0)
                       stop("unsupported polynomial power")
                   switch(as.character(e2),
                          "0" = 1,
                          "1" = e1,
                      {
                          p <- q <- polynomial(e1)
                          for(i in 2:e2)
                              p <- p * q
                          as.complex(p)
                      })
               },
               "%%" = {
                   if(l2 == 1)
                       0
                   else {
                       p <- rev(e1)
                       q <- rev(e2)
                       while(length(p) >= l2) {
                           d <- p[1]/q[1]
                           p[1:l2] <- p[1:l2] - d * q
                           p <- p[-1]
                       }
                       if(length(p) == 0) 0 else rev(p)
                   }
               },
               "==" = return(l1 == l2 && all(e1 == e2)),
               "!=" = return(l1 != l2 || any(e1 != e2)),
               stop("unsupported operation on polynomials"))
    polynomial(e1.op.e2)
}

Summary.polynomial <-
function(..., na.rm = FALSE)
{
    ok <- switch(.Generic, sum = , prod = TRUE, FALSE)
    if(!ok)
        stop(gettextf("Generic '%s' not defined for \"%s\" objects.",
                      .Generic, .Class))
    switch(.Generic,
           "sum" = accumulate("+", as.polynomial(0), polylist(...)),
           "prod" = accumulate("*", as.polynomial(1), polylist(...)))
}

Math.polynomial <-
function(x, ...)
{
    switch(.Generic,
           round = ,
           signif = ,
           floor = ,
           ceiling = ,
           trunc = polynomial(NextMethod(.Generic)),
           stop(paste(.Generic, "unsupported for polynomials")))
}

as.character.polynomial <- function(x, decreasing = FALSE, ...)
{
    p <- unclass(x)
    lp <- length(p) - 1
    names(p) <- 0:lp
    p <- p[p != 0]

    if(length(p) == 0) return("0")

    if(decreasing) p <- rev(p)

    signs <- rep("+ ", length(p));
    signs[1] <- "";

    np <- names(p)
    p <- as.character(p)
    p[p == "1" & np != "0"] <- ""

    pow <- paste("x^", np, sep = "")
    pow[np == "0"] <- ""
    pow[np == "1"] <- "x"
    stars <- rep.int("*", length(p))
    stars[p == "" | pow == ""] <- ""
    paste(signs, "(", p, ")", stars, pow, sep = "", collapse = " ")
}

print.polynomial <-
function(x, digits = getOption("digits"), decreasing = FALSE, ...)
{
    p <- as.character.polynomial(signif(x, digits = digits),
                                 decreasing = decreasing)
    pc <- nchar(p)
    ow <- max(35, getOption("width"))
    m2 <- 0
    while(m2 < pc) {
        m1 <- m2 + 1
        m2 <- min(pc, m2 + ow)
        if(m2 < pc)
            while(substring(p, m2, m2) != " " && m2 > m1 + 1)
                m2 <- m2 - 1
        cat(substring(p, m1, m2), "\n")
    }
    invisible(x)
}

as.function.polynomial <-
function(x, ...)
{
    a <- rev(coef(x))
    w <- as.name("w")
    v <- as.name("x")
    ex <- call("{", call("<-", w, 0))
    for(i in seq_along(a)) {
        ex[[i + 2]] <- call("<-", w, call("+", a[1], call("*", v, w)))
        a <- a[-1]
    }
    ex[[length(ex) + 1]] <- w
    f <- function(x) NULL
    body(f) <- ex
    f
}

poly.orth <-
function(x, degree = length(unique(x)) - 1, norm = TRUE)
{
    at <- attr(poly(x, degree), "coefs")
    a <- at$alpha
    N <- at$norm2
    x <- polynomial()
    p <- list(polynomial(0), polynomial(1))
    for(j in 1:degree)
        p[[j + 2]] <-
            (x - a[j]) * p[[j + 1]] - N[j + 1]/N[j] * p[[j]]
    p <- p[-1]
    if(norm) {
        sqrtN <- sqrt(N[-1])
        for(j in 1 + 0:degree) p[[j]] <- p[[j]]/sqrtN[j]
    }
    class(p) <- "polylist"
    p
}