Tridiag / libspline

A small C library for natural and clamped cubic spline interpolation. libspline fits a C2-continuous piecewise cubic polynomial to irregularly spaced data using a tridiagonal linear system and the Thomas algorithm.

Features

• Natural cubic spline with zero second derivatives at endpoints
• Clamped cubic spline with user-supplied endpoint slopes
• O(n) spline construction via Thomas algorithm
• O(m log n) evaluation for m query points using binary search
• Minimal, portable C API with no external dependencies beyond libm

Files

• spline.c - library implementation
• spline.h - public header and API declarations
• test_spline.c - unit test suite
• Makefile - build and test targets

Build

From the project root:

make clean
make

This builds:

• libspline.so or libspline.dll (shared library)
• test_spline executable

Run tests

make test

The test suite verifies interpolation accuracy, boundary conditions, continuity, extrapolation behavior, error handling, and minimal input cases.

API

spline_fit

SplineStatus spline_fit(
    const double       *x,
    const double       *y,
    size_t              n,
    const SplineParams *params,
    SplineFit         **out
);

Fit a cubic spline through n points (x[i], y[i]). The x values must be strictly increasing.

• params == NULL uses SPLINE_NATURAL boundary conditions.
• For clamped boundary conditions, pass a SplineParams struct with bc = SPLINE_CLAMPED, and set d0, dn.
• On success, *out is set to a newly allocated SplineFit *.

spline_eval

SplineStatus spline_eval(
    const SplineFit *s,
    const double    *x_new,
    size_t           m,
    int              extrapolate,
    double          *y_out
);

Evaluate the fitted spline at m query points.

• extrapolate != 0: extrapolate using the endpoint polynomial piece.
• extrapolate == 0: out-of-range points produce NAN.

spline_free

void spline_free(SplineFit *s);

Free the spline handle.

spline_get_n

size_t spline_get_n(const SplineFit *s);

Return the number of knots in the fitted spline.

spline_status_str

const char *spline_status_str(SplineStatus status);

Return a static string describing a status code.

Example

#include "spline.h"
#include <stdio.h>

int main(void)
{
    double x[] = {0.0, 1.0, 2.0, 3.0};
    double y[] = {0.0, 1.0, 4.0, 9.0};
    size_t n = 4;

    SplineFit *s = NULL;
    if (spline_fit(x, y, n, NULL, &s) != SPLINE_OK) {
        fprintf(stderr, "spline_fit failed\n");
        return 1;
    }

    double xq[] = {0.5, 1.5, 2.5};
    double yq[3];
    spline_eval(s, xq, 3, 1, yq);

    for (size_t i = 0; i < 3; i++) {
        printf("x=%.2f y=%.6f\n", xq[i], yq[i]);
    }

    spline_free(s);
    return 0;
}

Algorithm

The spline fit solves a tridiagonal system for the second derivative values M[i] at each knot. The system is:

lo[i] * M[i-1] + diag[i] * M[i] + up[i] * M[i+1] = rhs[i]

for interior knots, where:

• lo[i] = h[i-1]
• diag[i] = 2 * (h[i-1] + h[i])
• up[i] = h[i]
• rhs[i] = 6 * ((y[i+1] - y[i]) / h[i] - (y[i] - y[i-1]) / h[i-1])

Boundary conditions are handled by special first and last rows.

flowchart TD
    A[Input knots x[], y[]] --> B[Compute interval widths h[i]]
    B --> C[Build tridiagonal system]
    C --> D[Set boundary conditions]
    D --> E[Thomas forward sweep]
    E --> F[Back substitution for M[i]]
    F --> G[Compute cubic coefficients a, b, c, d]
    G --> H[Return fitted spline handle]

    subgraph evaluation [Spline evaluation]
      I[Query x_new] --> J[Find segment via binary search]
      J --> K[Evaluate cubic polynomial]
      K --> L[Return y_out]
    end

    H --> I

Loading

Notes

• The library is intentionally simple and focuses on correctness and portability.
• test_spline.c contains regression tests that can be extended for new cases.
• Use make clean to remove generated artifacts.
• GitHub publishing is supported via .gitignore, .gitattributes, and a GitHub Actions CI workflow in .github/workflows/ci.yml.

License

No license is included in this repository. Add one if you plan to redistribute or publish the code.