A solver for arbitrary systems of linear equations

slineq.c

/*
14. 6. 2021

This piece of C99 code is a set of rudimentary functions for manipulating 
linear functions of an aribtrary number of coefficients. This includes
a function for solving systems of an arbitrary amount of equations using
the substitution method.

Note that I aimed for the absolute best elegance of the code rather
than performance, which would require additional "if" checks, function 
arguments, restrictions et cetera.

The main function, which can be found at the bottom of the file, contains
a simple example program using the slineq_solve function.

Explanations of the inner workings of the functions are provided in the 
code itself.

You can do whatever you want with this code with or without crediting me,
github.com/jacobsebek
Improvements, insights and suggestions are encouraged!
*/

// The maximum number of coefficients in a linear equation
// (Thus the maximum number of unknowns + 1)
#define LINEQ_MAX_COEFFS 4

// The type used for the coefficient values
typedef float real_t;

// A structure denoting a linear equation, it represents
// coeffs[0] * x1 + coeffs[1] * x2 + ... = const.
typedef struct {
    real_t coeffs[LINEQ_MAX_COEFFS];
} lineq_s;

// A structure denoting a system of linear equations
typedef struct {
    lineq_s eqs[LINEQ_MAX_COEFFS];
} slineq_s;

// Solves a linear equation for a given term
lineq_s lineq_solve(lineq_s eq, const int ci) {

    const real_t coeff = eq.coeffs[ci];

    for (int i = 0; i < LINEQ_MAX_COEFFS; i++)
        eq.coeffs[i] /= coeff;

    return eq;
}

// Substitutes a term in a linear equation
lineq_s lineq_subst(lineq_s eq, lineq_s sub_eq, const int ci) {

    sub_eq = lineq_solve(sub_eq, ci);
    const real_t coeff = eq.coeffs[ci];

    for (int i = 0; i < LINEQ_MAX_COEFFS; i++)
        eq.coeffs[i] += -sub_eq.coeffs[i] * coeff;

    return eq;
}

// Solves an arbitrary system of linear equations
slineq_s slineq_solve(const slineq_s seqs) {

    slineq_s results = {0};

    // The first pass solves solves each equation for one coefficient
    for (int e = 0; e < LINEQ_MAX_COEFFS; e++) {

        lineq_s curr = seqs.eqs[e];

        // Substitute as many unknowns in the current equation
        for (int c = 0; c < LINEQ_MAX_COEFFS; c++)
            if (results.eqs[c].coeffs[c] != 0)
                curr = lineq_subst(curr, results.eqs[c], c);

        // Find an unknown that hasn't been solved yet
        // and is solvable from the current equation
        int r = 0;
        for (; r < LINEQ_MAX_COEFFS; r++)
            if (curr.coeffs[r] != 0 && results.eqs[r].coeffs[r] == 0)
                break;

        if (r == LINEQ_MAX_COEFFS) continue;

        // Add the solved unknown into results
        results.eqs[r] = lineq_solve(curr, r);
    }

    // The second pass solves the partially solved unknowns
    // After this step, we are left with a fully solved system
    for (int r = 0; r < LINEQ_MAX_COEFFS; r++) {

        // Substitute all the unknowns
        for (int c = 0; c < LINEQ_MAX_COEFFS; c++)
            if (c != r && results.eqs[c].coeffs[c] != 0)
                results.eqs[r] = lineq_subst(results.eqs[r], results.eqs[c], c);

    }

    return results;
}

#include <stdio.h>
int main() {

    slineq_s sys = {{
        {{1, 1, 1, -6}}, // a + b + c = 6 
        {{1, 2, 0, -5}}, // a + 2b = 5
        {{2, 0, 1, -5}}  // 2a + c = 5
    }};

    sys = slineq_solve(sys);

    // Expected output: "a: 1.00, b: 2.00, c: 3.00"
    printf("a: %.2f, b: %.2f, c: %.2f\n", 
        -lineq_solve(sys.eqs[0], 0).coeffs[3],
        -lineq_solve(sys.eqs[1], 1).coeffs[3],
        -lineq_solve(sys.eqs[2], 2).coeffs[3]
    );
    // Note that the arguments above could have been written in a simpler way:
    // "-sys.eqs[0].coeffs[3]"
    // because if the system is solvable, slineq_solve is guaranteed to return  
    // the result solved for 1x. When the system is not solvable, the system is
    // "solved" for 0x, making the result NaN, this is handled by lineq_solve.

    return 0;
}