#include <stdio.h>     //printf
#include <malloc.h>    //alloc+free
#include <stdlib.h>    //system
#include <math.h>      //fabs
#define CLR "cls"
//********************* globale Daten ***************************************
int ZA,                //****** Zeilen-Anfang
    ZE,                //****** Zeilen-Ende
    SA,                //****** Spalten-Anfang
    SE,                //****** Spalten-Ende
    *I;                //****** Index-mitLaeufer fuer X-Werte (SpaltenTausch)
double *A,             //****** Koeffizienten-Matrix
      *X,              //****** LoesungsVektor
      *Y;              //****** LoesungenVektor
//********************* Proceduren Anmeldung ********************************
void max(void);        //Suchen von Max / Tauschen der Zeilen+Spalten
void piv(void);        //rechne Fuehre PivotZeile aus
int  nul(void);        //rechne Fuehre PivotZeile aus
void fhl(char);        //FehlerAusgabe
//************************ M M  A A  I I  N N *******************************
void main(void){
int z,s,fehler,rg=0,ZM;                    //laufVar
//***************************************************************************
//********************** E I N G A B E T E I L ******************************
//***************************************************************************
 system(CLR);   fhl(0);                   //****** EINGABE DER NxM
 do{fhl(2);scanf("%d",&SE);if(SE<1) printf("\a");}while(SE<1);
 do{fhl(1);scanf("%d",&ZE);if(ZE<1) printf("\a");}while(ZE<1);
 //******************** FEHLER 0 ********************************************
 if(ZE<SE){fhl(6);exit(0);}                //UnterBestimmung Fehler[1]
 if(ZE*SE>6400) {fhl(3);exit(0);}          //Dimension wird zu gross
 //**************************************************************************
 A=(double*)malloc(ZE*SE*sizeof(double));  //Hole Matrix
 X=(double*)malloc(SE*  sizeof(double));   //Hole X-Vektor
 I=(int  *)malloc(SE*  sizeof(int  ));     //Hole IndexVektor
 Y=(double*)malloc(ZE*  sizeof(double));   //Hole LoesungenVektor
 if((A==NULL)||(X==NULL)||(I==NULL)||(Y==NULL))  //beim Holen Fehler aufgetreten
    {fhl(3);exit(0);}
 for(s=0;s<SE;s++) *(I+s)=s;        //Index von X-Vektor festlegen
 ZA=SA=0;                           //StartWerte des Algos
 ZM=ZE;                             //merke die ZeilenAnzahl
 system(CLR);   fhl(4);             //***********************************
 for(z=0;z<ZE;z++){                 //****** EINGABE DER MATRIX *********
  for(s=0;s<SE;s++){                //***********************************
    printf("a(%d,%d)= ",z+1,s+1);   //***********************************
    scanf("%lf",(A+z*SE+s));}}      //***********************************
 system(CLR);   fhl(5);             //***********************************
 for(z=0;z<ZE;z++){                 //*** EINGABE DES LÖSUNGENVEKTORS ***
  printf("y(%d)= ",z+1);            //***********************************
  scanf("%lf",(Y+z));}              //***********************************
//********************** Rechen-Schleife ************************************
 do{
   fehler=nul();                    //Pruefe auf Null-Zeilen
   if(fehler!=0) continue;          //...
   max();                           //suche Max-PivoElement
   //printf("nach Max"); for(z=0;z<ZE;z++){printf("\n");for(s=0;s<SE;s++){
   //printf("%le  ",*(A+(z*SE+s)));} printf("%le  ",*(Y+z));} getchar();
   piv();                           //eliminiere Spalte
   rg++;                            //Rang der Matrix jetzt
//   printf("nach elim");for(z=0;z<ZE;z++){printf("\n");for(s=0;s<SE;s++){
//   printf("%le  ",*(A+(z*SE+s)));} printf("%le  ",*(Y+z));} getchar();
   ZA++;SA++;                       //Matrix wird kleiner
 }while((ZA<ZE)&&(SA<SE));          //Matrix Pruefen ob Ende erreicht
 //*************************** Fehler 1 *************************************
 if(fehler==2) {fhl(7);exit(0);}
 if(fehler==0) {fehler=nul();  }//erkenne die letzte Zeile
 if(fehler==2) {fhl(7);exit(0);}
 if(rg<SE)     {fhl(8);exit(0);}
//******************************* Errechnen des X-Vektors *******************
 for(s=SE-1;s>=0;s--){
  *(X+s)=*(Y+s)/(*(A+SE*s+s));       //bestimme das X-Element [X = Y div Koeff]
  for(z=SE-2;z>=0;z--){              //von allen restlichen Y's wird Koeff*X
    *(Y+z)-=(*(A+SE*z+s))*(*(X+s));}}//gleich abgezogen
//*************************** Ausgabe ***************************************
 system(CLR);

 if(ZE<ZM)  printf("\n\t\tGleichungsSystem war ueberbestimmt.\n");
 if(ZE==ZM) printf("\n\t\tGleichungsSystem war optimal bestimmt.\n");
 for(s=0;s<SE;s++){                //ueber alle X's
    for(z=0;*(I+z)!=s;z++);        //suche Index <S> == alte Ordnung
    if(*(X+z)<0.0)  printf("x(%d) = %.5lf \t== %.15le\n",s+1,*(X+z),*(X+z)); //Ausgabe
    if(*(X+z)>=0.0) printf("x(%d) =  %.5lf \t==  %.15le\n",s+1,*(X+z),*(X+z));} //Ausgabe
 free(A); free(X); free(I); free(Y);}  //Freigabe der 'Mallox'
//********************** E N D E   M A I N **********************************

//***************************************************************************
//** UP1 ********* Suchen MAX: und  Tauschen der Spalten und Zeilen *********
//***************************************************************************
void max(void){          //Suchen von Max / Tauschen der Zeilen+Spalten
int s,z,sm=SA,zm=ZA;     //default
double tmp=*(A+ZA*SE+SA);//default
 for(z=ZA;z<ZE;z++){
   for(s=SA;s<SE;s++){
     if(fabs(*(A+z*SE+s))>fabs(tmp)) {sm=s;zm=z;tmp=*(A+z*SE+s);}}}
 for(s=SA;s<SE;s++){     //**** Zeilen-Tausch Matrix A
    tmp=*(A+ZA*SE+s);    //merke das ZeilenElement aus 1.Zeile
    *(A+ZA*SE+s)=*(A+zm*SE+s); //Überschreibe das ZeilenElement 1.Zl.mit MAX
    *(A+zm*SE+s)=tmp;}   //MAX-Zeile-Elem.= mit 1.Zeile
 tmp=*(Y+ZA);          //**** Zeilen-Tausch Vektor Y
 *(Y+ZA)=*(Y+zm);      //Überschreibe das ZeilenElement 1.Zl.mit MAX
 *(Y+zm)=tmp;          //MAX-Zeile-Elem.= mit 1.Zeile
 for(z=0;z<ZE;z++){    //**** Spalten-Tausch alle SpaltenElemente
    tmp=*(A+z*SE+SA);  //merke das SpaltenElement aus 1.Spalte
    *(A+z*SE+SA)=*(A+z*SE+sm); //Überschreibe das SpaltenElement 1.Sp.mit MAX
    *(A+z*SE+sm)=tmp;} //MAX-Spalte-Elem.= mit 1.Spalte
 zm=*(I+SA);           //**** Spalten-Tausch Vektor- I  (X-Index)
 *(I+SA)=*(I+sm);      //nehme zm als tmp, da *I == int; brauch zm nicht mehr
 *(I+sm)=zm;}           //**********************

//***************************************************************************
//** UP2 ********* Fuehre die PivotZeile aus; mit Probe auf 0.0 **************
//***************************************************************************
void piv(void){
int s,z,l1,l2,l3;                   //laufVar
double fak;                         //Pivot-Faktor
 for(z=ZA+1;z<ZE;z++){              //ueber alle Zeilen
   fak=*(A+z*SE+SA)/(*(A+ZA*SE+SA));//1.Sp-Elem div Pivot
   l1=l2=l3=0;                      //Exponent setzen; auf NUL Merken
   if(*(Y+z)!=0.0) {                //*** Exponent Minuend
     l1=(int)fabs(log10(fabs((*(Y+z)))));}
   if(fak*(*(Y+ZA))!=0.0) {         //*** Exponent Subtrahend
     l2=(int)fabs(log10(fabs(fak*(*(Y+ZA)))));}
   *(Y+z)  -=fak*(*(Y+ZA));          //Y-Vektor Anpassen
   if(*(Y+z)!=0.0) {                //*** Exponent Ergebnis
     l3=(int)fabs(log10(fabs((*(Y+z)))));}
   if( (((l1-l2)*(l1-l2))<2) &&     //*** bei starker Exponent-
         (((l1-l3)*(l1-l3))>25)) *(Y+z)=0.0;    //*** Abweichung NULL setzen
   for(s=SA;s<SE;s++){               //ueber alle Spalten
     l1=l2=l3=0;                     //Exponent setzen auf NUL Merken
     if(*(A+z*SE+s)!=0.0) {          //*** Exponent Minuend
       l1=(int)fabs(log10(fabs((*(A+z*SE+s)))));}
     if(fak*(*(A+(ZA*SE+s)))!=0.0) { //*** Exponent Subtrahend
       l2=(int)fabs(log10(fabs(fak*(*(A+(ZA*SE+s))))));}
     //****** printf("%le %le  %d %d",*(A+z*SE+s),fak*(*(A+(ZA*SE+s))),l1,l2);
     *(A+z*SE+s)-=fak*(*(A+(ZA*SE+s))); //Ausfuehren des Pivotzeile
     if(*(A+z*SE+s)!=0.0) {             //*** Exponent Ergebnis
       l3=(int)fabs(log10(fabs((*(A+z*SE+s)))));}
     //****** printf("  %le %d\n",*(A+z*SE+s),l3);
     if( (((l1-l2)*(l1-l2))<2) &&                  //*** bei starker Exponent
         (((l1-l3)*(l1-l3))>25)) *(A+z*SE+s)=0.0;}}//*** Abweichung NULL setzen
} //Spalten abgleichen
//***************************************************************************
//** UP3 ********* Loeschen der Nul-Zeilen ***********************************
//***************************************************************************
int nul(void){
int s,z,ma,my,ret=0;                 //LaufVar
double tmp;                          //AblegeVar
 for(z=ZA;z<ZE;z++){                 //Über alle Zeilen
   ma=0;my=0;                        //Zaehl-Var loeschen
   for(s=SA;s<SE;s++){               //zaehle die Nullen in der Zeile
     if(*(A+z*SE+s)==0.0) ma++;}     //wenn gefunden dann Zaehle hoch
   if(*(Y+z)==0.0) my++;             //ist Y==0 oder nich?
   if((ma==(SE-SA))&&(my==0)){       //Unterbestimmung; grober Fehler
     SA=SE;return(2);}               //bereite das Ende vor und Tschuess
   if((ma==(SE-SA))&&(my==1)){       //sind alle Elemente der Zeile == 0
     ZE--;                           //ZeilenEnde neu anpassen
     for(s=SA;s<SE;s++){        //**** Zeilen-Tausch der Matrix A
      tmp=*(A+z*SE+s);          //merke das ZeilenElement aus aktuelle Zeile
      *(A+z*SE+s)=*(A+ZE*SE+s); //Überschreibe das ZeilenElement
      *(A+ZE*SE+s)=tmp;}        //MAX-Zeile-Elem.= mit 1.Zeile
     tmp=*(Y+z);                //**** Zeilen-Tausch Vektor Y
     *(Y+z)=*(Y+ZE);            //Überschreibe das ZeilenElement
     *(Y+ZE)=tmp;               //....
     z--;ret=1;}}               //aktuelle Zeile ist getauscht nochMal Zeile
return(ret);}                   //alles ok return(0);or(1);
//***************************************************************************
//** UP4 ********* Ausgabe von Saetzen und Fehlern ***************************
//***************************************************************************
void fhl(char nummer){
char *S[9]={           //Text-Ausgaben
"                            Gauss-Algo\n"                            //0
"                           ===========\n"                           //0
"                            _  _\n"                                 //0
"                        A * X= Y\n"                                 //0
"geg. sind Koeffizienten Matrix A (mit z-Zeilen & s-Spalten) und der\n"
"Loesungs-Vektor Y (mit z-Zeilen). Gesucht ist der Vektor X (mit s-Zeilen),\n"
"der dieses Gleichungssystem loest.\n\n"                              //0
"    Ú x[1]*a[1,1]   x[2]*a[1,2]  ...  x[s]*a[1,s] ¿ = Y[1]\n"       //0
"    ³ x[1]*a[2,1]   x[2]*a[2,2]  ...  x[s]*a[2,s] ³ = Y[2]\n"       //0
"    :  :      :      :     :     ...   :    :     : =  :  \n"       //0
"    À x[1]*a[z,1]   x[2]*a[z,2]  ...  x[s]*a[z,s] Ù = Y[z]\n\n",    //0
"\nGib die  ZeilenZahl von A  z:= ",                                 //1
"\nGib die SpaltenZahl von A  s:= ",                                 //2
"\a Matrix ist zu gross, zu wenig Speicher <EXIT>\n",                 //3
"\n gib die Werte der MATRIX A ein\n"                                //3
" Elemente-Format a(z,s), d.h. a(zeile,spalte)\n\n",                 //4
"\n gib die Werte des LoesungsVektors Y ein\n"                        //5
" Elemente-Format y(z), d.h. y(zeile)\n\n",                          //5
"\n\n\a\tDas GleichungsSystem ist  unterbestimmt, zur korrekten Aussage\n"
      "\tueber den  Vektor- X muessen  mindestens genausoviele Zeilen wie\n"
      "\tSpalten vorhanden sein. Das GleichungsSystem besitzt unendlich\n"
      "\tviele Loesungen.  <zu wenig Zeilen angegeben>\n",            //6
"\n\n\a\tDas Gleichungs- System kann  nicht geloest  werden. Die um\n"//7
      "\tY- Vektor  erweiterte  Koeffizienten- Matrix A  hat einen\n"//7
      "\thoeheren Rang, als nur  die Matrix A, dh. rg(A,Y) > rg(A),\n"//7
      "\tueber ein und den selben Vektor werden zwei oder mehr ver-\n"//7
      "\tschiedene Aussagen  gemacht. Fuer  z=2  und  s=1  ein Bsp.\n"//7
      "\t2*x=2 und 3*x=1\n",                                         //7
"\n\n\a\tDas GleichungsSystem kann nicht eindeutig geloest werden.\n" //8
      "\tEs gibt unendlich viele Loesungen. Die Gleichungen sind\n"   //8
      "\tlinear abhaengig.\n"};                                       //8
 printf("%s",S[nummer]);}
//***************************************************************************
//********************** ENDE PROGRAM GAUSS - ALGO **************************
//***************************************************************************
//***************************************************************************
//****** p r o g r a m m i e r t   v o n   d a v e   s u n   i n   14h ******
//***************************************************************************
//***************************************************************************
//************************** 30.11.1994 by Dave Sun *************************
//***************************************************************************