"Math.polynomial"<-
function(x, digits)
{
	x <- as.numeric(x)
	switch(.Generic,
		round = ,
		signif = as.polynomial(NextMethod(.Generic)),
		{
			warning("polynomial coerced to numeric vector")
			NextMethod(.Generic)
		}
		)
}
"Ops.polynomial"<-
function(e1, e2)
{
	if(missing(e2)) {
# Unary operators
		e2 <- e1
		e1 <- 0
	}
	e1 <- unclass(e1)
	e2 <- unclass(e2)
	switch(.Generic,
		"+" = poly.add(e1, e2),
		"-" = poly.subtract(e1, e2),
		"*" = poly.multiply(e1, e2),
		"/" = ,
		"%/%" = poly.divide(e1, e2),
		"%%" = poly.remainder(e1, e2),
		"^" = if(length(e2) > 1 || e2 < 0 || (e2 %% 1) != 0) {
			stop("non positive integer power requested")
		}
		else {
			switch(as.character(e2),
				"0" = as.polynomial(1),
				"1" = as.polynomial(e1),
				{
				  m <- e1
				  for(i in 2:e2)
				    m <- poly.multiply(m, e1)
				  m
				}
				)
		}
		,
		"==" = length(e1) == length(e2) && all(e1 == e2),
		"!=" = length(e1) != length(e2) || any(e1 != e2),
		"<" = length(e1) < length(e2) || (length(e1) == length(e2) && 
			e1[length(e1)] < e2[length(e2)]),
		">" = length(e1) > length(e2) || (length(e1) == length(e2) && 
			e1[length(e1)] > e2[length(e2)]),
		"<=" = length(e1) < length(e2) || (length(e1) == length(e2) && 
			e1[length(e1)] <= e2[length(e2)]),
		">=" = length(e1) > length(e2) || (length(e1) == length(e2) && 
			e1[length(e1)] >= e2[length(e2)]),
		stop(paste("\"", .Generic, "\"", 
			" is an undefined operation on polynomial arguments", 
			sep = "")))
}
"Summary.polynomial"<-
function(x, ...)
{
	warning("polynomial coerced to numeric vector")
	x <- as.numeric(x)
	NextMethod(.Generic)
}
"[.polynomial"<-
function(x, i, ...)
{
	x <- unclass(x)
	y <- NextMethod("[")
	if(length(y) == 0)
		return(numeric(0))
	if(any(duplicated(names(y))))
		y <- tapply(y, names(y), sum)
	my <- max(as.numeric(names(y)))
	r <- rep(0, my + 1)
	names(r) <- as.character(0:my)
	r[names(y)] <- y
	as.polynomial(r)
}
"[<-.polynomial"<-
function(x, ..., value)
{
	value <- as.numeric(value)
	x <- as.numeric(x)
	x <- NextMethod("[<-")
	as.polynomial(x)
}
"as.character.polynomial"<-
function(p)
make.character.polynomial(p)
"as.expression.polynomial"<-
function(x)
parse(text = as.character(x))
"as.function.polynomial"<-
function(p)
{
	f <- function(x, nam = as.character(x))
	NULL
	body <- paste("val <- ", make.character(p, style = "h"))
	body <- c(body, "if(!inherits(val, \"polynomial\")) names(val) <- nam", 
		"val")
	body <- parse(text = body)
	mode(body) <- "{"	# "}"
	f[[3]] <- body
	f
}
"as.polynomial"<-
function(p)
{
	if(is.polynomial(p))
		return(p)
	if(is.complex(p))
		p <- Re(p)
	p <- as.numeric(p)
	if(length(p) == 0)
		p <- 0
	if(any(p != 0)) {
		p <- rev(p)
		while(p[1] == 0) p <- p[-1]
		p <- rev(p)
	}
	else p <- 0
	names(p) <- 1:length(p) - 1
	class(p) <- "polynomial"
	p
}
"c.polynomial"<-
function(...)
{
	namlis <- unlist(lapply(list(...), function(x)
	names(as.polynomial(x))))
	y <- do.call("c", lapply(list(...), as.numeric))
	names(y) <- namlis
	if(any(duplicated(names(y))))
		y <- tapply(y, names(y), sum)
	my <- max(as.numeric(names(y)))
	r <- rep(0, my + 1)
	names(r) <- as.character(0:my)
	r[names(y)] <- y
	as.polynomial(r)
}
"change.origin"<-
function(p, o)
{
	if(!is.polynomial(p))
		stop(paste("\"", deparse(substitute(p)), "\"", 
			" is not a polynomial"))
	o <- unclass(o[1])
	r <- poly.value(p, o)
	m <- 1
	p <- deriv(p)
	while(p != 0) {
		r <- c(r, poly.value(p, o))
		m <- m + 1
		p <- as.polynomial(as.numeric(deriv(p))/m)
	}
	as.polynomial(r)
}
"deriv.polynomial"<-
function(expr, ...)
{
	expr <- unclass(expr)
	if(length(expr) == 1)
		return(as.polynomial(0))
	expr <- expr[-1]
	as.polynomial(expr * seq(along = expr))
}
"integral"<-
function(expr, ...)
UseMethod("integral")
"integral.default"<-
function(x, ...)
stop(paste("no integral method defined for objects of class", data.class(x)))
"integral.polynomial"<-
function(expr, ...)
{
	expr <- unclass(expr)
	expr <- as.numeric(expr)
	as.polynomial(c(0, expr/seq(along = expr)))
}
"is.polynomial"<-
function(p)
inherits(p, "polynomial")
"lines.polynomial"<-
function(p, length = 100, ...)
{
	pu <- par("usr")
	x <- seq(pu[1], pu[2], len = length)
	y <- poly.value(p, x, nam = NULL)
	out <- y <= pu[3] | y >= pu[4]
	y[out] <- NA
	lines(x, y, ...)
}
"make.character"<-
function(x, ...)
{
	UseMethod()
}
"make.character.default"<-
function(x, ...)
as.character(x)
"make.character.polynomial"<-
function(p, variable = "x", style = c("ordinary", "horner"))
{
	p <- unclass(p)
	lp <- length(p) - 1
	names(p) <- 0:lp
	p <- switch(match.arg(style),
		ordinary = {
			p <- p[p != 0]
			if(length(p) == 0)
				return(as.expression(0))
			sp <- c("- ", "+ ")[(p > 0) + 1]
			switch(sp[1],
				"- " = sp[1] <- "-",
				sp[1] <- "")
			np <- names(p)
			p <- as.character(abs(p))
			p[p == "1" & np != "0"] <- ""
			pow <- paste(variable, "^", np, sep = "")
			pow[np == "0"] <- ""
			pow[np == "1"] <- variable
			stars <- rep("*", length(p))
			stars[p == "" | pow == ""] <- ""
			p <- paste(sp, p, stars, pow, sep = "", collapse = " ")
		}
		,
		horner = {
			pn <- rev(p)
			p <- as.character(pn[1])
			for(i in pn[-1])
				p <- paste(i, "+", variable, "* (", p, ")")
		}
		)
	p
}
"monic"<-
function(p)
{
	p <- unclass(p)
	if(all(p == 0)) {
		warning("the zero polynomial has no monic form")
		return(as.polynomial(0))
	}
	p <- rev(p)
	while(p[1] == 0) p <- p[-1]
	p <- rev(p/p[1])
	as.polynomial(p)
}
"plot.polynomial"<-
function(p, xs = seq(0, 1, len = 100), type = "l", ...)
{
	deg <- length(p) - 1
	if(!missing(xs) || deg < 2 || deg > 45) {
		x <- xs
		Px <- poly.value(p, x)
		return(plot.default(x, Px, type = type, ...))
	}
# try to cover the "interesting" region
	oldopt <- options(warn = -1)
	on.exit(options(oldopt))
	xr <- range(0, Re(polyroot(p)), Re(polyroot(deriv(p))))
	options(oldopt)
	xr <- xr + c(-1, 1)/20 * max(diff(xr), 1)
	x <- seq(xr[1], xr[2], len = 100)
	Px <- poly.value(p, x)
	plot.default(x, Px, type = type, ...)
}
"points.polynomial"<-
function(p, length = 100, ...)
{
	pu <- par("usr")
	x <- seq(pu[1], pu[2], len = length)
	y <- poly.value(p, x, nam = NULL)
	out <- y <= pu[3] | y >= pu[4]
	y[out] <- NA
	points(x, y, ...)
}
"poly.add"<-
function(p, q)
{
	p <- unclass(p)
	q <- unclass(q)
	lp <- length(p)
	lq <- length(q)
	if(lp < lq)
		p <- c(p, rep(0, lq - lp))
	if(lp > lq)
		q <- c(q, rep(0, lp - lq))
	as.polynomial(p + q)
}
"poly.divide"<-
function(p, q)
{
	p <- rev(unclass(p))
	q <- rev(unclass(q))
	lq <- length(q)
	if(lq == 1)
		return(as.polynomial(rev(p/q)))
	r <- numeric(0)
	while(length(p) >= length(q)) {
		d <- p[1]/q[1]
		r <- c(d, r)
		p[1:lq] <- p[1:lq] - d * q
		p <- p[-1]
	}
	if(length(r) == 0)
		r <- 0
	as.polynomial(r)
}
"poly.from.roots"<-
function(...)
{
	poly.from.zeros(...)
}
"poly.from.values"<-
function(x, px, tol = sqrt(.Machine$double.eps))
{
	if(is.matrix(px)) {
		lis <- list()
		for(i in 1:dim(px)[2])
			lis[[paste("P", i, sep = "")]] <- poly.from.values(x, 
				px[, i], tol)
		return(lis)
	}
	if(any(toss <- duplicated(x))) {
		crit <- max(tapply(px, x, function(x)
		diff(range(x))))
		if(crit > tol)
			warning("some duplicated points are inconsistent")
		keep <- !toss
		px <- px[keep]
		x <- x[keep]
	}
	if(length(x) != length(px))
		stop("x and px(x) do not match in length!")
	r <- 0
	for(i in seq(along = x))
		r <- r + (px[i] * unclass(poly.from.zeros(x[ - i])))/prod(x[i] - 
			x[ - i])
	r[abs(r) < tol] <- 0
	r <- rev(r)
	while(length(r) >= 1 && !r[1]) r <- r[-1]
	as.polynomial(rev(r))
}
"poly.from.zeros"<-
function(...)
{
	r <- unlist(list(...))
	p <- 1
	for(x in r)
		p <- c(0, p) - c(x * p, 0)
	as.polynomial(p)
}
"poly.multiply"<-
function(p, q)
{
	p <- unclass(p)
	q <- unclass(q)
	if(length(q) == 1 || length(p) == 1)
		return(as.polynomial(p * q))
	m <- outer(p, q)
	m <- tapply(m, row(m) + col(m) - 2, sum)
	class(m) <- "polynomial"
	m
}
"poly.orth"<-
function(x, degree = length(unique(x)) - 1, nam = paste("P", 1:degree, sep = ""
	))
{
        if(!exists("poly", mode="function") ||
	   !exists("poly.raw", mode="function"))
	  stop("Need `poly()' and `poly.raw()' which are NOT yet available")
	x0 <- seq(min(x), max(x), length = degree + 1)
	xb <- poly(x, degree)
	xb <- poly.raw(x0, degree, coef = attr(xb, "coefs"))
	pl <- list()
	for(i in 1 + 1:degree)
		pl[[i - 1]] <- poly.from.values(x0[1:i], xb[1:i, i - 1])
	names(pl) <- nam
	pl
}
"poly.remainder"<-
function(p, q)
{
	p <- rev(unclass(p))
	q <- rev(unclass(q))
	lq <- length(q)
	if(lq == 1)
		return(polynomial(0))
	r <- numeric(0)
	while(length(p) >= length(q)) {
		d <- p[1]/q[1]
		r <- c(d, r)
		p[1:lq] <- p[1:lq] - d * q
		p <- p[-1]
	}
	as.polynomial(rev(p))
}
"poly.subtract"<-
function(p, q)
{
	p <- unclass(p)
	q <- unclass(q)
	lp <- length(p)
	lq <- length(q)
	if(lp < lq)
		p <- c(p, rep(0, lq - lp))
	if(lp > lq)
		q <- c(q, rep(0, lp - lq))
	as.polynomial(p - q)
}
"poly.value"<-
function(p, x, nam = as.character(x))
{
	if(!is.polynomial(p))
		stop(paste(deparse(substitute(p)), "is not a polynomial"))
	v <- 0
	p <- unclass(p)
	for(pj in rev(p))
		v <- x * v + pj
	if(!inherits(v, "polynomial"))
		names(v) <- nam
	v
}
"polynomial"<-
function(p)
as.polynomial(p)
"print.polynomial"<-
function(p0, variable = "x", quote, ...)
{
	ow <- max(35, unlist(options("width")) - 1)
	p <- make.character(signif(p0, digits = options("digits")[[1]]), 
		variable = variable)
	m2 <- 0
	np <- nchar(p)
	while(m2 < np) {
		m1 <- m2 + 1
		m2 <- m2 + ow
		if(m2 < np)
			while(m2 > m1 && substring(p, m2, m2) != " ") m2 <- m2 - 
				  1
		cat(substring(p, m1, m2), "\n")
	}
	invisible(p0)
}
