Polynomial Library in BsplineQuantReg

Alexandre Abbes

2026-08-20

Introduction

The BsplineQuantReg package provides a comprehensive set of tools for polynomial manipulation, implemented in pure R with a consistent interface. All polynomial functions follow the decreasing power order convention, where a polynomial \(p(x) = a_0 + a_1 x + a_2 x^2\) is represented as c(a_2, a_1, a_0).

This vignette covers: - Basic polynomial operations (addition, multiplication, evaluation, differentiation) - Piecewise polynomial (PP) form - Callable and non-callable PP objects - Extracting parameters from PP objects

Basic Polynomial Operations

Polynomial Multiplication

# (1 + x) * (1 + x) = 1 + 2x + x^2
p1 <- c(1, 1)  # 1 + x
p2 <- c(1, 1)  # 1 + x
product <- polymul(p1, p2)
print(product)  # c(1, 2, 1) → 1 + 2x + x^2
## [1] 1 2 1
# (x^2 + 2x + 1) * (x - 1) = x^3 + x^2 - x - 1
p3 <- c(1, 2, 1)  # x^2 + 2x + 1
p4 <- c(1, -1)    # x - 1
polymul(p3, p4)
## [1]  1  1 -1 -1

Polynomial Addition

# (1 + x) + (1 - x) = 2
polyadd(c(1, 1), c(1, -1))
## [1] 2 0
# (x^2 + 1) + (x + 1) = x^2 + x + 2
polyadd(c(1, 0, 1), c(1, 1))
## [1] 1 1 2

Polynomial Evaluation

# P(x) = 1 + x + x^2
p <- c(1, 1, 1)  # x^2 + x + 1
poly_eval(p, c(0, 1, 2))  # returns c(1, 3, 7)
## [1] 1 3 7
# Evaluate at many points
x <- seq(-2, 2, length.out = 10)
y <- poly_eval(c(1, 0, -1), x)  # 1 - x^2
plot(x, y, type = "l", main = "1 - x^2")

Polynomial Differentiation

# P(x) = x^2 → P'(x) = 2x
p <- c(1, 0, 0)  # x^2
polyderiv(p, der = 1)  # c(2, 0) → 2x
## [1] 2 0
# Second derivative: P''(x) = 2
polyderiv(p, der = 2)  # c(2)
## [1] 2
# Higher order derivatives
p <- c(1, 2, 3, 4)  # 4x^3 + 3x^2 + 2x + 1
polyderiv(p, der = 2)  # 24x + 6
## [1] 6 4

Polynomial Reduction

# Remove leading zeros
reduce_pol(c(0, 0, 1, 2, 1))  # c(1, 2, 1) → x^2 + 2x + 1
## [1] 1 2 1
reduce_pol(c(0, 0, 0, 5))     # c(5) → constant polynomial
## [1] 5

Taylor Basis Change

# Convert (x-1)^2 to expansion around 0
# (x-1)^2 = x^2 - 2x + 1
coeffs_a <- c(1, 0, 0)  # (x-1)^2 in basis centered at a=1
change_polynomial_base_taylor(coeffs_a, a = 1, b = 0)
## [1]  1 -2  1
# Returns c(1, -2, 1) → x^2 - 2x + 1

Piecewise Polynomial (PP) Form

The PP form represents a function defined piecewise by polynomials on different intervals.

Creating a PP Object

# Create a piecewise polynomial with two intervals
# Interval 1: x^2 on [0, 1]
# Interval 2: 2x - 1 on [1, 2]
coeff <- matrix(c(
  1, 0, 0,   # x^2
  0, 2, -1   # 2x - 1
), nrow = 2, byrow = TRUE)

knots <- c(0, 1, 2)

pp <- makpp(coeff, knots)
print(pp)
## Piecewise Polynomial (PP) (non-callable)
## ================================
##   $degree: 2 
##   $knot: 0 1 2 
##  coefficients dimension: 2 x 3 
## 
##   $coeff:
##     Intervals 1 : 1, 0, 0 
##     Intervals 2 : 0, 2, -1 
##  Usage: pp_eval(pp, x_values)
# Evaluate the PP form
x <- seq(0, 2, length.out = 100)
y <- evalpp(pp, x)

plot(x, y, type = "l", main = "Piecewise Polynomial")
abline(v = knots, col = "red", lty = 2)

Non-Callable PP Objects

# A non-callable PP object is a list with components:
# - coeff: matrix of polynomial coefficients
# - knot: knot positions
# - degree: polynomial degree

pp <- makpp(coeff, knots, callable = FALSE)
class(pp)  # "non_callable_pp"
## [1] "non_callable_pp"
# Access components
pp$coeff
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2   -1
pp$knot
## [1] 0 1 2
pp$degree
## [1] 2

Callable PP Objects

A callable PP object can be evaluated directly like a function.

# Create a callable PP
pp_call <- makpp(coeff, knots, callable = TRUE)
class(pp_call)  # "callable_pp" "function"
## [1] "callable_pp" "function"
# Evaluate directly
x <- seq(0, 2, length.out = 10)
y <- pp_call(x)

# Print shows information
print(pp_call)
## Callable Piecewise Polynomial (PP) Object
## ==========================================
##   Degree: 2 
##   Intervals: 2 
##   Knots: 3 
##  coefficients dimension: 2 x 3 
## 
##   $coeff:
##     Intervals 1 : 1, 0, 0 
##     Intervals 2 : 0, 2, -1 
##  Usage: pp(x_values) or evalpp(pp, x_values)
# Can be used in plots
plot(pp_call, xlim = c(0, 2))

displaying polynomial in a human readable form

Introduction

The show_poly() function displays a polynomial as a human-readable mathematical equation. It supports both canonical and local bases, and can convert between different expansion points using Taylor’s formula.

A polynomial \(P(x) = \sum_{k=0}^d c_k (x-a)^k\) represented by coefficients c(c_d, c_{d-1}, ..., c_0) (decreasing power order) can be displayed in any basis \((x-b)^k\).

Basic Usage

Canonical Basis (Default)

By default, show_poly() displays the polynomial in the canonical basis \(1, x, x^2, ...\):

# P(x) = 3x^3 - 2x^2 + x - 5
p <- c(3, -2, 1, -5)
show_poly(p)
## [1] "3*x^3-2*x^2+1*x^1-5"
# "3*x^3 - 2*x^2 + 1*x - 5"

Working with Local Bases

Display in Local Basis

To display a polynomial in the local basis \((x-a)^k\), use the b parameter:

# P(x) = 3(x-2)^3 - 2(x-2)^2 + (x-2) - 5
p <- c(3, -2, 1, -5)
show_poly(p, b = 2)
## [1] "3*(x-2)^3+16*(x-2)^2+29*(x-2)^1+13"
# "3*(x-2)^3 - 2*(x-2)^2 + 1*(x-2) - 5"

Converting Between Bases

The a parameter specifies the basis in which the coefficients are given. The function automatically converts using change_polynomial_base_taylor():

# Coefficients are in basis (x-2)^k, display in canonical basis
p <- c(3, -2, 1, -5)  # In basis (x-2)^k
show_poly(p, a = 2, b = 0)
## [1] "3*x^3-20*x^2+45*x^1-39"
# This expands to: 3x^3 - 18x^2 + 37x - 27

Examples with Different Bases

Positive and Negative Centers

# Center at positive value
p <- c(1, 0, 2)  # (x-3)^2 + 2
show_poly(p, b = 3)
## [1] "1*(x-3)^2+6*(x-3)^1+11"
# "(x-3)^2 + 2"

# Center at negative value
p <- c(1, 0, 2)  # (x+2)^2 + 2
show_poly(p, b = -2)
## [1] "1*(x+2)^2-4*(x+2)^1+6"
# "(x+2)^2 + 2"

Controlling Output Format

Number of Digits

p <- c(1/3, -2/7, 1/5)
show_poly(p, digits = 2)
## [1] "0.33*x^2-0.29*x^1+0.2"
# "0.33*x^2 - 0.29*x + 0.2"

show_poly(p, digits = 4)
## [1] "0.3333*x^2-0.2857*x^1+0.2"
# "0.3333*x^2 - 0.2857*x + 0.2"

Negative Coefficients

The function handles negative coefficients gracefully:

# All coefficients negative
p <- c(-3, -2, -1)  # -3x^2 - 2x - 1
show_poly(p)
## [1] "-3*x^2-2*x^1-1"
# "-3*x^2 - 2*x - 1"

# Mixed coefficients
p <- c(3, -2, 1, -5)  # 3x^3 - 2x^2 + x - 5
show_poly(p)
## [1] "3*x^3-2*x^2+1*x^1-5"
# "3*x^3 - 2*x^2 + 1*x - 5"

Working with B-splines

When working with B-splines, show_poly() helps visualize the polynomial pieces:

# Create a B-spline and display its polynomial form
sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75)

# Convert to PP
pp <- Bsplinetopp(basis, callable = FALSE)

# Display each piece
cat("Piece 1: ", show_poly(pp$coeff[1, ], b = 0), "\n")
## Piece 1:  -133.3333*x^3+91.6667*x^2-15*x^1+1
cat("Piece 2: ", show_poly(pp$coeff[2, ], b = 0), "\n")
## Piece 2:  28.4203*x^3-28.3333*x^2+4*x^1+1.15
cat("Piece 3: ", show_poly(pp$coeff[3, ], b = 0), "\n")
## Piece 3:  43.0325*x^3-2.7551*x^2-5.3265*x^1+0.5673

Integration with show_pp()

The show_poly() function is used internally by show_pp() and show_pp() to display piecewise polynomials:

# Create a PP with two intervals
knot <- c(0, 0.5, 1)
coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE)
pp <- makpp(coeff, knot)

# Display in local basis
show_pp(pp, local = TRUE, digits = 3)
##      [,1]                  [,2]                       
## [1,] "  [0.0000, 0.5000] " "1*x^2+2*x^1+0.5"          
## [2,] "  [0.5000, 1.0000] " "0*(x-0.5)^2+1*(x-0.5)^1-1"

Summary

Parameter Description Default
obj Polynomial coefficients (decreasing power order) Required
a Base of input coefficients \((x-a)^k\) 0
b Base for output display \((x-b)^k\) 0
digits Number of significant digits 4
verbose Print additional information FALSE

Displaying PP Equations

The package provides convenient functions to display PP objects as human-readable mathematical equations. The show_pp() and show_pp() functions format the polynomial expressions on each interval, with support for both canonical and local bases.

# Create a PP with two intervals
knot <- c(0, 0.5, 1)
coeff <- matrix(c(1, 2, 0.5, 0, 1, -1), nrow = 2, ncol = 3, byrow = TRUE)
pp <- makpp(coeff, knot)

# Display in canonical basis (1, x, x², ...)
show_pp(pp, local = FALSE)
##      [,1]                  [,2]             
## [1,] "  [0.0000, 0.5000] " "1*x^2+2*x^1+0.5"
## [2,] "  [0.5000, 1.0000] " "0*x^2+1*x^1-1.5"
# [0.000, 0.500] 0.5x^2 + 2x + 1
# [0.500, 1.000] -1x^2 + 1x + 0

# Display in local basis ((x-a)^i)
show_pp(pp, local = TRUE)
##      [,1]                  [,2]                       
## [1,] "  [0.0000, 0.5000] " "1*x^2+2*x^1+0.5"          
## [2,] "  [0.5000, 1.0000] " "0*(x-0.5)^2+1*(x-0.5)^1-1"
# [0.000, 0.500] 0.5(x-0)^2 + 2(x-0) + 1
# [0.500, 1.000] -1(x-0.5)^2 + 1(x-0.5) + 0

For spline objects, the conversion to PP is automatic:

# Create a B-spline and display its polynomial form
sn <- c(0, 0, 0, 0, 0.3, 0.6, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- c(1, -0.5, 2, 0.5, -1, 0.75)

show_pp(basis, local = TRUE, verbose = TRUE)
## Use makpp or make_spline to format obj
## NULL
# PP Information:
# Degree: 3
# Knots: 0, 0.3, 0.6, 1
# Number of intervals: 3
# 
# [0.000, 0.300] 1(x-0)^3 + 2.5(x-0)^2 + 1.5(x-0) + 1
# [0.300, 0.600] -4.63(x-0.3)^3 + ... 
# [0.600, 1.000] ...

Customization Options The display functions offer several customization options:

local: Control the basis (TRUE for local, FALSE for canonical)

digits: Number of significant digits (default: 4)

verbose: Display additional PP information

Callable PP Objects

A callable PP object can be evaluated directly like a function.

Create a callable PP

pp_call <- makpp(coeff, knots, callable = TRUE) class(pp_call) # “callable_pp” “function”

Evaluate directly

x <- seq(0, 2, length.out = 10) y <- pp_call(x)

Can be used in plots

plot(pp_call, xlim = c(0, 2))

Extrapolating values

When evalpp() is called with x values outside the knot range (x < knots[1] or x > knots[length(knots)]), it automatically performs extrapolation by extension using the border polynomial piece (the first or last polynomial in the PP structure). For example, if the knot vector is c(0, 1, 2), points with x < 0 will be evaluated using the polynomial from the first interval [0, 1], while points with x > 2 will use the polynomial from the last interval [1, 2]. This is equivalent to “extending” the first and last piece of the spline. A warning message is printed to alert the user when extrapolation occurs.

PP with two intervals: x^2 on [0,1], 2x-1 on [1,2]

coeff <- matrix(c(1, 0, 0, 0, 2, -1), nrow = 2, byrow = TRUE) knots <- c(0, 1, 2) pp <- makpp(coeff, knots, callable = TRUE)

Evaluate inside and outside the knot range

x_test <- c(-0.5, 0.5, 1.5, 2.5) y_test <- pp(x_test) data.frame(x = x_test, y = y_test) # The values at -0.5 and 2.5 are extrapolated


### Extracting Parameters


``` r
# For callable PP objects, use get_parameters()
params <- get_parameters(pp_call)
print(params$degree)
## [1] 2
print(params$knot)
## [1] 0 1 2
print(params$coeff)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2   -1
# For non-callable PP, access directly
pp=makpp(pp,callable=FALSE,verbose=TRUE)
## Already a 'non_callable_pp'
## 
pp$degree
## [1] 2
pp$knot
## [1] 0.0 0.5 1.0
pp$coeff
##      [,1] [,2] [,3]
## [1,]    1    2  0.5
## [2,]    0    1 -1.0

Advanced: PP Form Manipulation

Extracting a Single Interval Polynomial

# Get polynomial for interval i
pp=makpp(pp,callable=FALSE)
interval_poly <- pp$coeff[1, ]  # First interval
print(interval_poly)
## [1] 1.0 2.0 0.5
# Evaluate on a subinterval
x_sub <- seq(0, 1, length.out = 20)
y_sub <- poly_eval(interval_poly, x_sub)

Combining PP Objects

# Create two PP objects and combine them
# This is useful for constructing complex piecewise functions
coeff1 <- matrix(c(1, 0), nrow = 1)
coeff2 <- matrix(c(0, 1), nrow = 1)

pp1 <- makpp(coeff1, c(0, 1))
pp2 <- makpp(coeff2, c(1, 2))

# Combine coefficients and knots
combined_coeff <- rbind(coeff1, coeff2)
combined_knots <- c(0, 1, 2)
pp_combined <- makpp(combined_coeff, combined_knots)

# Evaluate combined function
x <- seq(0, 2, length.out = 100)
y_combined <- evalpp(pp_combined, x)
plot(x, y_combined, type = "l")
abline(v = c(0, 1, 2), col = "red", lty = 2)

Visualizing PP with Basis Functions

# Create a PP with multiple pieces
knots <- seq(0, 1, length.out = 5)
degree <- 2
n_pieces <- length(knots) - 1

# Generate random coefficients
coeff <- matrix(rnorm(n_pieces * (degree + 1)), nrow = n_pieces)

pp <- makpp(coeff, knots, callable = TRUE)

# Evaluate and plot
x <- seq(0, 1, length.out = 200)
y <- pp(x)

plot(x, y, type = "l", lwd = 2,
     main = "Random Piecewise Quadratic")
abline(v = knots, col = "red", lty = 2)
grid()

Comparison: PP vs B-spline

# Create a B-spline and convert to PP form
sn <- c(0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1)
basis <- Bspline_base(sn, degree = 3)
basis$coeff <- runif(basis$n_splines)

# Convert to PP
pp_from_bspline <- Bsplinetopp(basis)

# Both represent the same function
x <- seq(0, 1, length.out = 100)
y_bspline <- spline_eval(basis, x)
y_pp <- evalpp(pp_from_bspline, x)

# They should match
max(abs(y_bspline - y_pp))
## [1] 2.220446e-16

Summary