# 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)

# Print shows information

print(pp_call)

# 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