Need To Generalize code and solve determinant

[ahmedm512]

Dear all,
Thanks to "tpb" he provided me with code that read the equations and sort them in an array and I have added up operations that are done on this equations, now I am facing only 2 problems:

1-I need to make 'number one' optional in user input and there should be no multiples if user entered 2 of x1 they should be added up and if he entered 2 results they should be added up
ex: 2x1+3x2-x1-12=10 should be saved as 1x1+3x2=22

3x1+5x2-2x2-x1=12 should be saved as 2x1+3x2=12

x1-x2+3x3-1=10 should be saved as 1x1-1x2+3x3=11

2- solving determinant there is a problem with the code it gave weird numbers.

[tpb]

Here's a new way to read in the equations. It can handle missing coefficients (assumed to be 1) and terms on either side of the equals sign. I don't know anything about determinants so I can't help you with that.

#include <iostream>
#include <iomanip>
#include <sstream>
#include <string>
using namespace std;

const int MAX_EQS   = 20;
const int MAX_COEFS = 20;

struct Equation {
    int coef[MAX_COEFS];
    int result;
};

void read_eq(string &line, Equation &eq) {
    for (int i = 0; i < MAX_COEFS; i++)
        eq.coef[i] = 0;
    eq.result = 0;

    istringstream sin(line);
    int sign = 1, rhs = -1;
    char ch = 0;

    while (sin >> ch) {
        if (ch == '=')
            sign = rhs = 1;
        else if (ch == '+')
            sign = 1;
        else if (ch == '-')
            sign = -1;
        else {
            int coef = 1, subscript = -1;
            if (ch == 'x')
                sin >> subscript;
            else {
                sin.unget();
                sin >> coef >> ch;
                if (ch != 'x')
                    sin.unget();
                else
                    sin >> subscript;
            }
            if (subscript == -1)
                eq.result += coef * sign * rhs;
            else
                eq.coef[subscript] += coef * sign * -rhs;
        }
    }
}

void print_eq(Equation &eq) {
    bool first = true;
    for (int i = 0; i < MAX_COEFS; i++)
        if (eq.coef[i] != 0) {
            if (!first && eq.coef[i] > 0)
                cout << '+';
            if (eq.coef[i] == -1)
                cout << '-';
            else if (eq.coef[i] != 1)
                cout << eq.coef[i];
            cout << 'x' << i;
            first = false;
        }
    cout << '=' << eq.result << '\n';
}

int main() {
    Equation eq[MAX_EQS];
    int n = 0;
    string line;

    while (getline(cin, line))
        read_eq(line, eq[n++]);

    for (int i = 0; i < n; i++)
        print_eq(eq[i]);
}

[lastchance]

@ahmedm512,

You have a lot of problems with your determinant routine. Many of them stem from the fact that you are (apparently) counting your equations from 0 and your coefficients from 1. Thus in (one of) the base cases of your recursion, for a 1x1 matrix:
return eq[0].coef[1];
This is asking for trouble.

Here is a code for computing determinants by (recursively) expanding in cofactors, which is what you are trying to do. (Note that there are far better ways for large matrices - more about that below). Note first two things:
- the routine is templated, but I would expect it in reality to be using doubles, not ints (more below)
- the det() routine takes a ONE-DIMENSIONAL array as parameter (for considerable convenience); you can call it with &A[0][0] if you are sending a 2-d array. The way you have set up your problem you would have to write a wrapping routine to first set up a 1-d or 2-d array.

#include <iostream>
using namespace std;


//=====================================================================


template <typename T> T det( int N, T arr[] )              // Returns determinant of NxN array
{                                                          // Note: arr[] is mapped to 1-d: index = N * i + j
   if ( N == 1 ) return arr[0];                            // Base case: 1x1 matrix

   T result = 0;
   for ( int c = 0, sign = 1; c < N; c++, sign = -sign)    // Expand as sum of cofactors across top row
   {
      T *minor = new T[(N-1)*(N-1)];                       // Set up subarray to calculate minor for element [0][c]
      int a = N, m = 0;                                    // Indices in main array and subarray; incremented when used
      for ( int i = 1; i < N; i++ )                        // Omit row 0, on which you are expanding
      {
         for ( int j = 0; j < N; j++ )
         {
            if ( j != c ) minor[m++] = arr[a];             // Omit column for this cofactor
            a++;
         }
      }
      result += sign * arr[c] * det( N - 1, minor );       // Recursive call to evaluate minor (thence cofactor)
      delete [] minor;
   }

   return result;
}


//=====================================================================


int main ()
{
   const int N = 4;

   double A[N][N] = { {  1,  2,  3,  4 },
                      { -2,  5,  5,  7 },
                      {  1,  9, 10,  3 },
                      {  2,  2,  4,  3 } };

   cout << "Array:\n";
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) cout << '\t' << A[i][j];
      cout << '\n';
   }


   // The following assumes contiguous storage of successive rows;
   //    det() actually takes a ONE-DIMENSIONAL array of size N x N
   cout << "\nDeterminant is " << det( N, &A[0][0] ) << '\n';
}

Array:
	1	2	3	4
	-2	5	5	7
	1	9	10	3
	2	2	4	3

Determinant is 74

But now consider some other things. You should seriously rethink your code

(1) Your coefficients are all ints. OK, this may be true about the input. However, it does NOT mean that the solution is in integers.
Consider
x1-x2=1
x1+x2=2
The coefficients are all ints ... but the solution is x1=1.5, x2=0.5. Clearly NOT ints.

I strongly advise you to use doubles.


(2) What is your aim at the end of this? If you are evaluating the determinant to determine whether the matrix is singular then well and good. However, if your aim is to SOLVE those equations then finding the determinant first is a waste of time. There are far better solution methods.

The fact that you CAN write down the solution in terms of determinants does NOT mean that you should.


(3) No serious scientist would actually input a set of linear equations the way you are doing. He would simply input the matrix coefficients themselves as they perfectly express the set of equations. Fundamentally, you are setting up (and presumably aiming to solve) a set of equations of the form
Ax=b
where A is a matrix and x and b are (mathematical) vectors.

If you have been told to input equations this way then I guess you will have to stick to it. If you haven't then I suggest that you seriously rethink your input.

Even if you want your string-based input, you should still store coefficients within your code as a 1-d or 2-d array: nothing more.

[ahmedm512]

@tpb thanks for your help



@lastchance

thanks for your response.

(Note that there are far better ways for large matrices - more about that below)

I would consider a better a way of solving det than the way I used, I don't know much about it I was trying to solve it by any means. I need more efficient way.

(1) Thanks for the advice , I will use doubles.

(2) My aim at the end is to solve the equations but first I am asked to apply on them some operations as shown in my code.

There are far better solution methods.

I would be thrilled if you showed me one, as I was considering solving the the equations by solving the determinant of x1 and dividing it on the main determinant (Cramer rule).

(3)I was told to input equations this way.

counting your equations from 0 and your coefficients from 1.

I make it like this as it helps alot when applying operations.

I thought I figured this problem when I shifted the index but apparently not.

you should still store coefficients within your code as a 1-d or 2-d array

you mean each equation in a single array ?? (Show me what is right as I am confused )


Finally, to be specific and not waste your time, My aim is a function that solves determinants of the array that contain the equations which is eq[].coef[].(or the one that will be extracted to, If i understand right.)

Thanks in advance

[lastchance]

ahmedm512 wrote:

I was considering solving the the equations by solving the determinant of x1 and dividing it on the main determinant (Cramer rule).

Unless you are intending to give your answers as rational fractions (like 273/987) then you would be better off using Gaussian elimination. Maybe with pivoting to avoid zeroes on the diagonals.

Store your equation set as a single array.

counting your equations from 0 and your coefficients from 1.
I make it like this as it helps alot when applying operations.

I'm not convinced. I don't mind whether arrays start from 0 or 1 ... as long as they all start from the same thing!

(3)I was told to input equations this way.

I'd love to know by whom!

(1) Thanks for the advice , I will use doubles.

If you intend to actually solve the equations your answers will generally come out as non-integers, so unless you have a fraction class this is sensible.

[ahmedm512]

Unless you are intending to give your answers as rational fractions (like 273/987) then you would be better off using Gaussian elimination. Maybe with pivoting to avoid zeroes on the diagonals.

Could you show me an example

I'm not convinced. I don't mind whether arrays start from 0 or 1 ... as long as they all start from the same thing!

I got your point I am convinced :D

I'd love to know by whom!

By the university Dr.

[lastchance]

ahmedm512 wrote:

Could you show me an example

https://en.wikipedia.org/wiki/Gaussian_elimination

//======================================
// [m x n] matrices are represented as 1-d arrays
// Row i goes from 0 to m-1.
// Col j goes from 0 to n-1.
// The (i,j) component is stored at 1-d index
//    n * i + j
//
// Storage order as for C++ arrays (column index varies fastest);
//======================================

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <algorithm>
using namespace std;

const double EPS = 1.0e-30;

//======================================


template <class T> bool GaussElimination( int n, T *A, T *B, T *X )  // Solve AX=B
{                                                                   
   // Create copies to avoid overwriting
   T *AA = new T[n*n];
   T *BB = new T[n];
   copy( A, A + n * n, AA );
   copy( B, B + n    , BB );

   // Row operations for i = 0, ,,,, n - 2 (n-1 not needed)
   for ( int i = 0; i < n - 1; i++ )
   {
      // Pivoting; swap rows such that no elements below [i][i] have larger magnitude
      int r = i;                                 // Find row r below with largest element in column i
      int ii = n * i + i;
      T maxA = abs( AA[ii] );
      for ( int k = i + 1; k < n; k++ )       
      {
         int ki = n * k + i;
         if ( abs( AA[ki] ) > maxA )
         {
            r = k;
            maxA = abs( AA[ki] );
         }
      }
      if ( r != i )                              // Swap rows i and r (remains of these rows for AA)
      {
         for ( int j = i; j < n; j++ ) swap( AA[n*r+j], AA[n*i+j] );
         swap( BB[i], BB[r] );
      }

      // Row operations to make upper-triangular
      if ( abs( AA[ii] ) < EPS ) return false;   // Singular matrix

      for ( int r = i + 1; r < n; r++ )          // On lower rows
      {
         T multiple = AA[r*n+i] / AA[ii];        // Multiple of ith row necessary to clear element in ith column
         for ( int j = i; j < n; j++ ) AA[r*n+j] -= multiple * AA[i*n+j];
         BB[r] -= multiple * BB[i];
      }

   }

   // Back-substitute
   for ( int i = n - 1; i >= 0; i-- )
   {
      X[i] = BB[i];
      for ( int j = i + 1; j < n; j++ )  X[i] -= AA[n*i+j] * X[j];
      X[i] /= AA[n*i+i];
   }

   delete [] AA;
   delete [] BB;

   return true;
}

//======================================

template <class T> void matMul( int mA, int nA, int nB, T *A, T *B, T *C )     // Matrix multiply: A * B = C
{
   for ( int i = 0; i < mA; i++ )
   {
      for ( int j = 0; j < nB; j++ )
      {
         int ij = nB * i + j;
         C[ij] = 0.0;
         for ( int k = 0; k < nA; k++ ) C[ij] += A[nA*i+k] * B[nB*k+j];
      }
   }
}

//======================================

template <class T> void showMatrix( int m, int n, T *A )
{
   const double ZERO = 1.0e-10;
   const string SPACE = "  ";
   const int w = 12;
   const int p = 4;
   cout << fixed << setprecision( p );

   for ( int i = 0; i < m; i++ )
   {
      for ( int j = 0; j < n; j++ )
      {
         T val = A[n*i+j];   if ( abs( val ) < ZERO ) val = 0.0;
         cout << setw( w ) << val << SPACE;
      }
      cout << '\n';
   }
}

//======================================

int main()
{
   const int N = 4;
   double A[N*N] = {  1,  2,  3,  4,
                     -2,  5,  5,  7,
                      1,  9, 10,  3,
                      2,  2,  4,  3  };
   double B[N] = { 1, 2, 3, 4 };
   double X[N] = { 0 };

   // SOLVE
   bool ok = GaussElimination( N, A, B, X );
   if ( !ok ) { cout << "Can't solve\n";   return 1; }

   // SUMMARY
   cout << "\n*** Solve AX = B by Gaussian elimination ***\n";
   cout << "\nA:\n";   showMatrix( N, N, A );
   cout << "\nB:\n";   showMatrix( N, 1, B );
   cout << "\nX:\n";   showMatrix( N, 1, X );

   // CHECKS
   double *AX = new double[N*1];
   matMul( N, N, 1, A, X, AX  );
   cout << "\n\n*** CHECK SOLUTION ***\n";
   cout << "\nAX:\n";   showMatrix( N, 1, AX );
   cout << "\nB :\n";   showMatrix( N, 1, B  );
   delete [] AX;
}

*** Solve AX = B by Gaussian elimination ***

A:
      1.0000        2.0000        3.0000        4.0000  
     -2.0000        5.0000        5.0000        7.0000  
      1.0000        9.0000       10.0000        3.0000  
      2.0000        2.0000        4.0000        3.0000  

B:
      1.0000  
      2.0000  
      3.0000  
      4.0000  

X:
     -2.1622  
     -4.4595  
      4.6757  
     -0.4865  


*** CHECK SOLUTION ***

AX:
      1.0000  
      2.0000  
      3.0000  
      4.0000  

B :
      1.0000  
      2.0000  
      3.0000  
      4.0000  

[ahmedm512]

I will try adding this to my code.
Many thanks for your assistance I appreciate it.

[jonnin]

Generally speaking you can pick off the eigenvalues & vectors with an efficient approach and then if absolutely necessary you can build the det from those, but IMHO (and I did nearly 20 years of linear algebra controls code) the det is about useless (by the time you compute it, you already knew everything you needed about the matrix and could do anything you needed to do with it). Dets are kind of cute for 2x2 and 3x3, and after that, I strongly advise using eigens.

you may also want to set up your gauss elim to be able to return the inverse if you need it.

[ahmedm512]

One of the taskes asked from me is to print the value of det. And then solve the equations.
Could you show me an ex. to approch det. Using eigenvalues as i am asked to solve upto 100 equations.
Actually I couldnt get the value of det. Till now 😅, i couldnt apply @lastchance ex. on my code.
As this is my first project , I am missing some concepts and I try to understand first so I can replace values.

[lastchance]

You just have to invoke the function as det( N, A ) if A is a 1-d array (with N*N elements) or
det( N, &A[0][0] ) if A is a 2-d array with N elements per dimension.

It is up to you to provide that array A, and it depends on how the rest of your code is structured.

I fail to see what eigenvalues have got to do with this problem.

[jonnin]

A couple of quick links:
(the first one does not work all the time, but its neat)
https://math.stackexchange.com/questions/507641/show-that-the-determinant-of-a-is-equal-to-the-product-of-its-eigenvalues

http://www.math.northwestern.edu/~len/LinAlg/chap2.pdf

http://www.math.harvard.edu/archive/20_spring_05/handouts/ch05_notes.pdf

the code to get the eigens is nontrivial ... probably take a good 2 days to cook it all up and debug it if you don't use an established library... takes a lu decomp and SOR eigen solver and maybe a third piece. I used the old intel math kernel library back when I did this, so the bits and pieces you need are a bit fuzzy as I didnt write it so much as glued it together.

Whether this is useful for you or not, I do not know. Its probably a sledgehammer for a fly, if you are doing schoolwork :)

These may not be the best links either, I just glanced at it. But google around on alternative ways to get the det, see what pops out. I stumbled across this from the other direction back when.. I was looking for a better way to get the eigens and read a bit on getting them from the det, but then looked into that and getting the det (for large matrices) and it wasn't the right way (ended up using SOR).

[ahmedm512]

Thanks all for your concern, actually this is too advanced for my project but it is helpful for my knowledge, and I have learned through this forum more than I learned from the whole semester :D, I appreciate it.

@tpb @lastchance @jonnin