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

float abs(float x){
       if (x>0){
          return x;
          }
       else {
            return (-1)*x;
            }
       }

int main(){
    
    int i,j,k,l,m,M,n,N,r,y;
   
printf("bitte Zeilenanzahl eingeben: \n");
scanf ("%d",&n);
printf("bitte Spaltenanzahl eingeben: \n");
scanf ("%d",&m);
      
      y=0;
      M=n;
      if (m<n) {M=m;}
      if (m!=n) {y++;}

    float D[n][m],E[n][m],A[n][m],c[n][m],a[n][m],L[n][M],R[M][m],S[n-1][n][n],P[n][n][n],q;
    
    //Matrizen S und P mit Einheitsmatrizen vorbelegen:
    for (k=1;k<=n-1;k++){
        for (i=1;i<=n;i++){
            for (j=1;j<=n;j++){
                if (i==j) {S[k][i][j]=1;    P[k][i][j]=1;        }
                else      {S[k][i][j]=0;    P[k][i][j]=0;        }
            }
        }
    }    
    // A eingeben:        
         for (i=1;i<=n;i++){
             for (j=1;j<=m;j++){
                 printf("bitte Eintrag A[%d][%d] eingeben:", i, j);
                 scanf("%f",&A[i][j]);
                 a[i][j]=A[i][j]; //für Aufgabenteil 4d 
                 c[i][j]=A[i][j]; //für Aufgabenteil 5b
             }
         }    

    // A auf Symmetrie überprüfen::        
    if (y==0){
         for (i=1;i<=n;i++){
             for (j=1;j<=n;j++){
                 if (c[i][j]!=c[j][i]){y++;}
             }
         } 
    }
    
    // A auf positive Definitheit überprüfen:
    if (y==0){
              //Gauß:     
        for (k=1;k<=n-1;k++){    
           for (i=k+1;i<=n;i++){
               c[i][k]=c[i][k]/c[k][k];
               for (j=k+1;j<=m;j++){
                   c[i][j]=c[i][j]-c[k][j]*c[i][k];
               }
           }
       }
       for (i=1;i<=n;i++){
           if (c[i][i]<=0){y++;}
       }
    }
    
if (y==0){              
            // Choleski:
                        
     float LSchlange[n][n],LSchlangeTransponiert[n][n],f;
     // LSchlange mit Nullen füllen:         
    for (i=1;i<=n;i++){
        for (j=1;j<=n;j++){
            LSchlange[i][j]=0;
        }
    }     
    // LSchlange verändern:     
     for (k=1;k<=n;k++){
         f=A[k][k];
         for (i=1;i<k;i++){
             f=f-LSchlange[k][i]*LSchlange[k][i];
         }
         LSchlange[k][k]=sqrt(f);
         for (j=k+1;j<=n;j++){
             f=A[j][k];
             for (l=1;l<k;l++){
                 f=f-LSchlange[j][l]*LSchlange[k][l];
             }
             LSchlange[j][k]=f/LSchlange[k][k];
         }
    }
    // Ausgabe von LSchlange:  
    printf ("\nLSchlange ergibt sich zu : \n \n");
    
    for (i=1;i<=n;i++){
        for (j=1;j<=n-1;j++){
            printf ("  %f\t", LSchlange[i][j]);
        }
        printf ("  %f\n", LSchlange[i][n]);
    } 
    
    // LSchlangeTransponiert definieren:
    for (i=1;i<=n;i++){
        for (j=1;j<=n;j++){
            LSchlangeTransponiert[i][j]=LSchlange[j][i];
        }
    }
    // Ausgabe von LSchlangeTransponiert:  
    printf ("\nLSchlangeTransponiert ergibt sich zu : \n \n");
    
    for (i=1;i<=n;i++){
        for (j=1;j<=n-1;j++){
            printf ("  %f\t", LSchlangeTransponiert[i][j]);
        }
        printf ("  %f\n", LSchlangeTransponiert[i][n]);
    }
       //Multiplikation von LSchlange und LSchlangeTransponiert:
   float B[n][n],z;
   for (i=1;i<=n;i++){
       for (j=1;j<=n;j++){
           z=0;
               for (k=1;k<=n;k++){
                   z=z+LSchlange[i][k]*LSchlangeTransponiert[k][j];
               }
           B[i][j]=z;
       }
   }
   //Ausgabe von LSchlange*LSchlangeTransponiert:  
        
        printf ("\nA-Schlange ergibt sich zu : \n \n");
    
    for (i=1;i<=n;i++){
        for (j=1;j<=n-1;j++){
            printf ("  %f\t", B[i][j]);
        }
        printf ("  %f\n", B[i][n]);
    }    
    
    //Berechnung von A-LSchlange*LSchlangeTransponiert:
    
    float C[n][n];
    
    for (i=1;i<=n;i++){
        for (j=1;j<=n;j++){
            C[i][j]=a[i][j]-B[i][j];
        }
    }
    
        //Ausgabe von A-LSchlange*LSchlangeTransponiert:
    
    printf ("\nDie Differenz von A und A-Schlange ergibt sich zu : \n \n");
    
    for (i=1;i<=n;i++){
        for (j=1;j<=m-1;j++){
            printf ("  %f\t", C[i][j]);
        }
        printf ("  %f\n", C[i][m]);
    } 
}  
      
else {
           //Gauß:
     
    // L und R mit Pivotisierung in A berechnen: 
       N=n;
       if (m<n){N=m+1;}     
       for (k=1;k<=N-1;k++){
           //Pivotierung:
           r=k;
           q=abs(A[k][k]);
           for (l=k+1;l<=n;l++){
               if (abs(A[l][k])>q){
               r=l;
               q=abs(A[l][k]);
               }
           }
       S[k][k][k]=0;
       S[k][r][r]=0;
       S[k][k][r]=1;
       S[k][r][k]=1;               
       for (i=1;i<=n;i++){
           for (j=1;j<=m;j++){
               D[i][j]=0;
               for (l=1;l<=n;l++){
                   D[i][j]=D[i][j]+S[k][i][l]*A[l][j];
               }
           }
       } 
       for (i=1;i<=n;i++){
           for (j=1;j<=m;j++){
               A[i][j]=D[i][j];
           }
       }    
       //Gauß:         
           for (i=k+1;i<=n;i++){
               A[i][k]=A[i][k]/A[k][k];
               for (j=k+1;j<=m;j++){
                   A[i][j]=A[i][j]-A[k][j]*A[i][k];
               }
           }
       } 
       
          // Ausgabe von A:  
          printf ("\nA ergibt sich zu : \n \n");    
    
    for (i=1;i<=n;i++){
        for (j=1;j<=m-1;j++){
            printf ("  %f\t", A[i][j]);
        }
        printf ("  %f\n", A[i][m]);
    }
     
        // Ausgabe von a:  
        printf ("\na ergibt sich zu : \n \n");    
    
    for (i=1;i<=n;i++){
        for (j=1;j<=m-1;j++){
            printf ("  %f\t", a[i][j]);
        }
        printf ("  %f\n", a[i][m]);
    }
            
    // P berechnen:
    for (k=1;k<=n-1;k++){
         for (i=1;i<=n;i++){
           for (j=1;j<=n;j++){
               P[k+1][i][j]=0;
               for (l=1;l<=n;l++){
                   P[k+1][i][j]=P[k+1][i][j]+P[k][i][l]*S[n-k][l][j];
               }
           }
         }
    } 
    // L mit Nullen füllen:         
    for (i=1;i<=n;i++){
        for (j=1;j<=M;j++){
            L[i][j]=0;
        }
    } 
    
        //R mit Nullen füllen:         
    for (i=1;i<=M;i++){
        for (j=1;j<=m;j++){
            R[i][j]=0;
        }
    }
               
    // L und R mit den Einträgen aus A füllen:    
    for (i=1;i<=n;i++){
        for (j=1;j<i;j++){
              L[i][j]=A[i][j];
        }
    }
    
    for (i=1;i<=M;i++){
        L[i][i]=1;
    }
        
    for (i=1;i<=M;i++){
        for (j=m;j>=i;j--){
            R[i][j]=A[i][j];
        }      
    }    
    // Ausgabe von L:  
    printf ("\nL ergibt sich zu : \n \n");
    
    for (i=1;i<=n;i++){
        for (j=1;j<=M-1;j++){
            printf ("  %f\t", L[i][j]);
        }
        printf ("  %f\n", L[i][M]);
    }    
   // Ausgabe von R:  
   printf ("\nR ergibt sich zu : \n \n");    
    
    for (i=1;i<=M;i++){
        for (j=1;j<=m-1;j++){
            printf ("  %f\t", R[i][j]);
        }
        printf ("  %f\n", R[i][m]);
    } 
    // Ausgabe von P:
    printf ("\nP ergibt sich zu : \n \n");    
    
    for (i=1;i<=n;i++){
        for (j=1;j<=n-1;j++){
            printf ("  %f\t", P[n][i][j]);
        }
        printf ("  %f\n", P[n][i][n]);
    }   
    //Aufgabe 4c)

          /*    float b[n],xSchlange[n],x[n];

    for (i=1;i<=n;i++){
            printf("\nb[%d] eingeben:",i);
            scanf("%f",&b[i]);
    }
    
           //Algorithmus:

    double summe; 
    for (i=1;i<=n;i++){
        xSchlange[i]=1/(L[i][i])*(b[i]-summe);   
        summe=0;     
        for (j=1;j<=i;j++){   
        summe=summe+xSchlange[j]*L[i+1][j];
        }
    }

    summe=0;
        for (i=n;i>=1;i--){
        x[i]=(1/R[i][i])*(xSchlange[i]-summe);
        summe=0;
                for (j=i;j<=n;j++){   
                summe=summe+x[j]*R[i-1][j];
                }
        }
 
        //Ausgabe von xSchlange
 
    for (i=1;i<n;i++){
        printf("%f  ",xSchlange[i]); 
        }  
        printf("%f  \n",xSchlange[n]);
        

        // Ausgabe der Lösung x:
        
    printf("Die Lösung von Ax=b ist:");    
    for (i=1;i<n;i++){
        printf("%f\t",x[i]);
        }
        printf("%f  \n",x[n]);
        */        
        
        // Aufgabe 4d

   //Multiplikation von L und R:

   float B[n][m],z;
   for (i=1;i<=n;i++){
       for (j=1;j<=m;j++){
           z=0;
               for (k=1;k<=M;k++){
                   z=z+L[i][k]*R[k][j];
               }
           B[i][j]=z;
       }
   } 
   
   //Multiplikation von P und A:

   for (i=1;i<=n;i++){
       for (j=1;j<=m;j++){
           z=0;
               for (k=1;k<=n;k++){
                   z=z+P[n][i][k]*a[k][j];
               }
           E[i][j]=z;
       }
   }
     
   //Ausgabe von L*R:  
        
        printf ("\nA-Schlange ergibt sich zu : \n \n");
    
    for (i=1;i<=n;i++){
        for (j=1;j<=m-1;j++){
            printf ("  %f\t", B[i][j]);
        }
        printf ("  %f\n", B[i][m]);
    }    
    
    //Berechnung von P*A-L*R:
    
    float C[n][n];
    
    for (i=1;i<=n;i++){
        for (j=1;j<=m;j++){
            C[i][j]=E[i][j]-B[i][j];
        }
    }
    
        //Ausgabe von A-L*R:
    
    printf ("\nDie Differenz von P*A und A-Schlange ergibt sich zu : \n \n");
    
    for (i=1;i<=n;i++){
        for (j=1;j<=m-1;j++){
            printf ("  %f\t", C[i][j]);
        }
        printf ("  %f\n", C[i][m]);
    }   
}      
return 0;
}