Linear equations C ver

linear.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

void vec_add_times(float *vec, float *other, float times, int n)
{
    while(n--) {
       vec[n] += other[n] * times;
    }
}

void vec_mul(float *vec, float multiplier, int n)
{
    while(n--) {
        vec[n] *= multiplier;
    }
}

float *vec_new(int n) {
    return (float *)malloc(sizeof(float) * n);
}

void vec_copy(float *dst, float *src, int n)
{
    memcpy(dst, src, sizeof(float)*n);
}

void vec_print(float *vec, int n)
{
    int i;
    printf("(");
    for(i=0; i<n-1; i++) {
        printf("%f, ", vec[i]);
    }
    printf("%f)", vec[n-1]);
}

int main(int argc, char *argv[])
{
    int n, i, j;
    float **mat;

    printf("Input n of unknowns:");
    scanf("%d", &n);
    mat = (float**)malloc(sizeof(float*) * n);

    printf("Input coefficients:\n");
    for(i=0; i<n; i++) {
        mat[i] = vec_new(i+1);
        for(j=0; j<n+1; j++) {
            scanf("%f", &mat[i][j]);
        }
    }

    float *tmpv = vec_new(i+1);
    for(i=0; i<n; i++) {
        if(!mat[i][i]) { // cur coef is zero
            for(j=i+1; j<n; j++) {
                if(mat[j][i]) { // if some other is not
                    float *v;
                    v = mat[i]; mat[i] = mat[j]; mat[j] = v; // swap them
                    break;
                }
            }
            if(j == n) continue; // all the coefs are zero
        }

        vec_mul(mat[i], 1.0/mat[i][i], n+1);
        
        for(j=i+1; j<n; j++) {
            vec_add_times(mat[j], mat[i], -mat[j][i], n+1);
        }
    }

    int rank;

    for(i=n-1; i>=0; i--) {
        for(j=0; j<n; j++) {
            if(mat[i][j]) {
                rank = i + 1;
                goto done;
            }
        }
        if(mat[i][n]) { // coef all zero but not b
            printf("No solution!\n");
            return 1;
        }
    }
    done:

    for(i=1; i<rank; i++) {
        for(j=0; j<i; j++) {
            vec_add_times(mat[j], mat[i], -mat[j][i], n+1);
        }
    }

    printf("Solution = \n");
    for(i=0; i<n; i++) {
        tmpv[i] = mat[i][n];
    }
    vec_print(tmpv, n);
    for(i=0; i<n-rank; i++) {
        for(j=0; j<rank; j++) {
            tmpv[j] = -mat[j][rank+i];
        }
        for(; j<n; j++) {
            tmpv[j] = 0.0;
        }
        tmpv[rank+i] = 1.0;
        printf("+ c%d * ", i);
        vec_print(tmpv, n);
    }
    return 0;
}