Programa que resuelve sistemas de ecuaciones

Hace unos meses tuve que hacer para la uni un programa que resolviera sistemas de ecuaciones usando el metodo de Gauss-Jordan.
Esto fue lo que hice, esta bastante explicado, si se te ocurre alguna mejora coméntanosla.

El programa esta incompleto y aun tiene muchos fallos, pero es un buen comienzo para uno aun mejor, funciona bien para matrices "normales" pero en cuanto lo pones a resolver matrices con muchos ceros al principio falla. lo intentare solucionar.

#include <iostream>
#include <math.h>
using namespace std;

typedef long double BIG;

int f;          // nº de filas
int c;          // nº de columnas
int i,j,k;      // variables para los bucles

void imprimir ( BIG **M )  // Funcion que imprime la matriz
{
    for (int i=0; i<f; i++)
    {
        for (int j=0; j<c; j++)
            cout<&lt;
        M[i][j] <<"\t";
        cout << endl;
    }
    cout << endl;
};

int main()
{
    BIG Mt[3][3][4]= {{{1,1,-1,1},{8,3,-6,1},{-4,-1,3,1}}, // Matrices de prueba
        {{1,-9,5,33},{1,3,-1,-9},{1,-1, 1,5}},
        {{1,-9,5,33},{1,3,-1,-9},{1,3,-1,-8}}
    };
    int t,s,min,max; // variables usadas en el menu y en la obtencion de datos
top:
    cout<<"0->Matrices de prueba"<<endl<<"1->Matriz aleatoria"<<endl<<"2->Introducir manualmente"<<endl;
    cin >> t;
    if (t==0)
    {
        cout<<"Tipo de sistema:"<<endl<<"0->Compatible determinado"<<endl<<"1->Compatible indeterminado"<<endl<<"2->Incompatible"<<endl;
        cin >> s;
        f = 3;
        c = 4;
    }
    else if (t==1)
    {
        cout << "Numero de filas y columnas:";
        cin >> f >> c;
        cout << endl << "Rango de los numeros, min. y max.:";
        cin >> min >> max;
    }
    else if (t==2)
    {
        cout << "Numero de filas y columnas:";
        cin >> f >> c;
    }
    else goto top;

    // Creamos la Matriz
    BIG **M = (BIG**)malloc(sizeof(BIG*) * f);
    for (i=0; i<f; i++)
        M[i] = (BIG*)malloc(sizeof(BIG) * c);

    // Rellenamos la matriz
    for (i=0; i<f; i++)
        for (j=0; j<c; j++)
            if (t==0)        M[i][j] = Mt[s][i][j];   //Matriz predefinida
            else if (t==1)   M[i][j] = rand()%max+min; //Numeros aleatorios
            else            cin >> M[i][j];           //Forma manual

    cout<<"Matriz inicial:"<<endl;
    imprimir(M);

    cout<<"Simplificando..."<<endl;
    min = (f<c?f:c)-1;      // reusamos la variable min para guardar el minimo entre las filas y columnas
    for (j=0; j<min; j++)   // recorremos la diagonal principal (menos la ultima fila o columna)
    {

        // Pivote()
        int p=j;            // establecemos al primer numero como pivote
        for (i=j+1; i<f; i++) // comprovamos si existe un mejor candidato
            if (fabs(M[i][j]) < fabs(M[p][j]) && M[i][j])
                p=i;

        // intercambia()
        if (p!=j)           // solo si el pivote esta en una fila que no es la primera
        {
            BIG* tf = M[j];
            M[j]=M[p];
            M[p]=tf;
            //imprimir(M);
        }

        // Hacer ceros
        for (k=j+1; k<f; k++)           // todas las filas a partir del pivote)
        {
            // combina()
            BIG mul = -(M[k][j])/(M[j][j]);
            for (int i=j; i<c; i++)
                M[k][i] += mul*M[j][i]; // recorremos todas las columnas de esa fila y operamos
            //imprimir(M);
            //cout<<"----"<<endl;
        }
        //cout<<"__________"<<endl;
    }

    imprimir(M);

    // Eliminando filas sobrantes (0,0,0,0)
    j=i=0;
    while (!M[f-1][c-1])
    {
        while (!M[f-1][i] && i<c)
            i++;
        if (i==c)
        {
            cout << "Eliminadas filas innecesarias..." << endl;
            f--;
            imprimir(M);
        }
        else break;
    }

    // Resolver:
    cout << "Solucion: ";

    if (M[f-1][c-1] && !M[f-1][c-2])
        cout << "Sistema incompatible (Sin solucion)" << endl;

    else if (f<c-1)                 // Si el numero de ecuaciones es = al numero de incognitas.
    {
        cout << "Sistema compatible indeterminado (infinitas soluciones)" << endl;
        //int pa=c-1-f;             //cantidad de parametros

        for (j=f; j<c-1; j++)
            for (i=0; i<f; i++)
                M[i][j] *= -1;      //pasamos los parametros al otro lado del "=", por lo que cambia de signo.

        for (i=f-1; i>=0; i--)      //recorremos las filas
        {
            for (j=f-1; j>i; j--)   //recorremos las columnas de la submatriz que no contiene los parametros ni el termino independiente.
            {
                for (int k=1; k<=c-f; k++) //para cada parametro y el termino independiente
                    M[i][c-k] -= M[i][j]*M[i+1][c-k];//sustituimos variables que ya sabemos
                M[i][j]=0;
            };
            for (k=1; k<=c-f; k++)
                M[i][c-k] = M[i][c-k] / M[i][j];
            M[i][j]=1;              // dejamos la diagonal principal a 1
        }
    }
    else
    {
        cout<<"Sistema compatible determinado (solucion unica)"<<endl;
        for (i=f-1; i>=0; i--)      // i -> desde la ultima fila hasta 0
        {
            for (j=c-2; i<j; j--)   // j -> desde la penultima columna ()
            {
                M[i][c-1] -= M[i][j] * M[j][c-1];//sustituimos y restamos(pasamos el valor a la derecha del "=")
                M[i][j] = 0;
            };
            M[i][c-1] = M[i][c-1] / M[i][j];
            M[i][j] = 1;
            //cout<<"x"<<i<<" = "<<M[i][c-1]<<endl;// resultados de las variables (modo alternatico)
        }
    }
    imprimir(M);
    cout<<endl<<"Final del programa"<<endl;
    system("PAUSE");
}