Weekend 1 hour project - needed simplest possible "least squares fitting of a polynomial"

polyfit.c

//  polyfit
//  Created by leo on 6/25/22.
//  Copyright © 2022 leok7v.github.io. All rights reserved.

#include <stdbool.h>
#include <assert.h>
#include <math.h>
#include <float.h>

// by Legendre in 1805:
// https://en.wikipedia.org/wiki/Least_squares
// with decent math explanaition here:s
// https://web.archive.org/web/20220303000301/https://neutrium.net/mathematics/least-squares-fitting-of-a-polynomial/
// the code is from https://github.com/natedomin/polyfit

static bool polyfit(double x[], double y[], int n, double poly[], int order) {
    // poly[] array of polynomial coefficients [0..order] for
    // y = sum(poly[i] * pow(x, i))
    enum { max_order = 5 }; // max order
    double B[max_order + 1] = { 0 }; // the column vector
    double P[(max_order + 1) * 2 + 1] = { 0 }; // pow(x)
    double A[(max_order + 1) * 2 * (max_order + 1)] = { 0 }; // reduction matrix
    bool done = true; // assume it can be done
    assert(order <= max_order && order < n);
    if (order > max_order || n <= order) {
        done = false; // invalid parameters
    } else {
        const int k1 = order + 1;
        const int k2 = k1 * 2;
        for (int i = 0; i < n; i++) {
            double xi = x[i];
            double yi = y[i];
            double powx = 1;
            for (int j = 0; j < k1; j++) { B[j] += yi * powx; powx *= xi; }
        }
        P[0] = n;
        for (int i = 0; i < n; i++) {
            double xi = x[i];
            double powx = x[i];
            for (int j = 1; j < k2 + 1; j++) { P[j] += powx; powx *= xi; }
        }
        for (int i = 0; i < k1; i++) {
            double* ai = A + i * k2;
            for (int j = 0; j < k1; j++) { ai[j] = P[i + j]; }
            ai[i + k1] = 1;
        }
        // Move the identity matrix portion of the redux matrix to the left side
        // finding the inverse of the left side of the reduction matrix.
        for (int i = 0; i < k1; i++) {
            double* ai = A + i * k2;
            double xi = ai[i];
            assert(fabs(xi) >= DBL_EPSILON);
            if (fabs(xi) >= DBL_EPSILON) {
                for (int k = 0; k < k2; k++) { ai[k] /= xi; }
                for (int j = 0; j < k1; j++) {
                    if (j != i) {
                        double* aj = A + j * k2;
                        double yj = aj[i];
                        for (int k = 0; k < k2; k++) { aj[k] -= yj * ai[k]; }
                    }
                }
            } else { // singularity:
                done = false; // matrix determinant == 0 (close enough)
            }
        }
        if (done) { // calculate the polynomial coefficients
            for (int i = 0; i < k1; i++) {
                double* ai = A + i * k2;
                for (int j = 0; j < k1; j++) {
                    double xi = 0;
                    for (int k = 0; k < k1; k++) { xi += ai[k + k1] * B[k]; }
                    poly[i] = xi;
                }
            }
        }
    }
    return done;
}

// for i in [0..n] returns sum(poly[i] * pow(x, i))
static double extrapolate(double x, double poly[], int n) {
    double y = 0;
    double powx = 1;
    for (int i = 0; i < n + 1; i++) { y += powx * poly[i]; powx *= x; }
    return y;
}

static void polyfit_test() {
    enum { order = 3 };
    enum { count = 5 }; // number of elements
    double coefficients[order + 1] = { 0 };
    // These inputs should result in the following approximate coefficients:
    //         0.5           2.5           1.0        3.0
    //    y = (0.5 * x^3) + (2.5 * x^2) + (1.0 * x) + 3.0
    {
        double x[count] = { 12.0, 77.8, 44.1, 23.6, 108.2 };
        double y[count] = { 1239.00, 250668.38, 47792.19, 7991.13, 662740.98 };
        bool done = polyfit(x, y, count, coefficients, order);
        assert(done);
    }
    static const double precision = 0.01; // acceptable error
    assert(fabs(3.0 - coefficients[0]) <= precision);
    assert(fabs(1.0 - coefficients[1]) <= precision);
    assert(fabs(2.5 - coefficients[2]) <= precision);
    assert(fabs(0.5 - coefficients[3]) <= precision);
    for (double x = 1.0; x < 200.0; x += 1.0) {
        double predicted = extrapolate(x, coefficients, order);
        double expected = (0.5 * pow(x, 3)) + (2.5 * pow(x, 2)) + (1.0 * x) + 3.0;
        double error = sqrt((expected - predicted) * (expected - predicted));
        assert(error < 0.1);
    }
}

int main(int argc, const char * argv[]) {
    polyfit_test();
    return 0;
}