I made a program to solve linear equations using Gauss Jordan elimination method. The program is working correctly but in some cases instead of giving the answer as 0, it returns a very small value.

#include<iostream>
#include<iomanip>
#include<cassert>
#define N 10
using namespace std;

//printing out the array
void print(float x[][N], int n){
    for(int i=0;i<n;i  ){
        for(int j=0;j<=n;j  )
            cout << setprecision(5) << setw(15) <<x[i][j];
        cout << endl;
    }
    cout << endl;
}

//to normalise the leading entries to 1
void normalize(float x[][N], int n, int i, int j){
    float fac = x[i][j];
    for(int k=0;k<=n;k  )
        x[i][k] /= fac;
}

//check if the leading entry is a zero
bool chk_zero(float x[][N], int n, int i){
    int j, k, c{0};
    if(x[i][i]==0){
    for(j=i;j<n-1;j  ){
        c=1;
        while(x[j c][i]==0 && i c<n){
            c  ;
            if(i c==n-1)
                assert(i c==n-1 && "Equation has no solution");
                return false;
            }

        for(k=0;k<=n;k  ){
            swap(x[i][k], x[i c][k]);
            }
        return true;

        }
    }
    return true;

}

//Gauss Jordan elimination method
void GaussJordan(float x[][N], int n){
    int i, j, k, c;
    float rat;

    for(i=0;i<n;i  ){
        //not taking the zero in pivot column case
        chk_zero(x, n, i);
        normalize(x, n, i, i);

        for(j=0;j<n;j  ){
            if (i != j){
               float fac{x[j][i]};
               for(k=0;k<=n;k  )
                   x[j][k] = x[j][k]-fac*x[i][k];
            }
        }
    }
}



int main(){
    float arr[][N] = { {0, 5, 1, 2},
                       {2, 11, 5, 3},
                       {1, 0, 0, 0} };

    int n = sizeof(arr)/sizeof(arr[0]);

    print(arr, n);

    GaussJordan(arr, n);

    cout << n << endl;

    print(arr, n);
    return 0;
}

the output I get is:

 1              0              0     2.3842e-08
 0              1              0            0.5
-0             -0              1           -0.5

The output I should get is:

 1              0              0              0
 0              1              0            0.5
-0             -0              1           -0.5

The value 2.3842e-08 should be a zero. Is it because of the precision of the floating point in C ?

If so, what should I do in order to round such low values to 0 without loss of data ?

Also, why is there a "-0" instead of a 0.

CodePudding user response：

The value 2.3842e-08 should be a zero. Is it because of the precision of the floating point in C ?

Yes, it is.

If so, what should I do in order to round such low values to 0 without loss of data ?

You can use an epsilon value to decide whether a float value is considered to be 0:

const float epsilon = 0.000001; // Depends on the use case of your application
if (std::abs(some_float) < epsilon) {
    // Treat some_float as 0
}

However, this approach has some cons. Read here for more.

Also, why is there a "-0" instead of a 0.

This is because of how float numbers are represented. Most C compilers use the IEEE 754 standard to implement floating point arithmetics.

Few things to add (unrelated to the question):

• You should avoid using using namespace std as it pollutes the global namespace. Especially for large codebases.
• Use std::array<> instead of raw arrays as they are more safe.

CodePudding user response：

The value 2.3842e-08 should be a zero. Is it because of the precision of the floating point in C ?

No. It is a feature of float point number representation in binary number system according IEEE-754 standard. It is implemented in hardware and is used by almost all the languages. All the floating point calculations are performed with errrors so you should set right rounding mode for FPU. See <fenv.h> for details.