#include <stdio.h>
#include <math.h>
#define Nmax 5
void MATGJ(int n,double a[Nmax+1][Nmax+2]);
/*-------------------------------------------------------------------*/

int main()
{
    /* 連立一次方程式[a]{x}={y}の解{x}を求める．a(i][j)=[a|y]を入力する */
    /* 入力出力結果                     */
    /* [a]={ 2][ 3][ 1} {y}={ 4} {x}={ 2} */
    /*     { 4][ 1][-3}     {-2}     {-1} */
    /*     {-1][ 2][ 2}     { 2]     { 3} */

    int i;
    int n=2;
    double a[Nmax+1][Nmax+2];

    a[0][0]= 2.0;a[0][1]= 3.0;a[0][2]= 1.0;a[0][3]= 4.0;
    a[1][0]= 4.0;a[1][1]= 1.0;a[1][2]=-3.0;a[1][3]=-2.0;
    a[2][0]=-1.0;a[2][1]= 2.0;a[2][2]= 2.0;a[2][3]= 2.0;

    printf("[a]=%e %e %e  {y}=%e\n",a[0][0],a[0][1],a[0][2],a[0][3]);
    printf("    %e %e %e      %e\n",a[1][0],a[1][1],a[1][2],a[1][3]);
    printf("    %e %e %e      %e\n",a[2][0],a[2][1],a[2][2],a[2][3]);

    MATGJ(n,a);

    printf("{x}=%e\n",a[0][3]);
    printf("    %e\n",a[1][3]);
    printf("    %e\n",a[2][3]);

    printf("リターンキーを押してください");
    getchar();
    return 0;
}
/*-------------------------------------------------------------------*/

void MATGJ(int n,double a[Nmax+1][Nmax+2])
{
    /* Gauss-Jordan法による連立一次方程式の解法 */
    int i,j,k,s;
    double p,d,max,dumy;

    for(k=0;k<=n;k++){
        max=0.0;
        s=k;
        for(j=k;j<=n;j++){
            if(fabs(a[j][k])>max){
                max=fabs(a[j][k]);
                s=j;
            }
        }
        for(j=0;j<=n+1;j++){
            dumy=a[k][j];
            a[k][j]=a[s][j];
            a[s][j]=dumy;
        }
        p=a[k][k];
        for(j=k;j<=n+1;j++){
            a[k][j]=a[k][j]/p;
        }
        for(i=0;i<=n;i++){
            if(i!=k){
                d=a[i][k];
                for(j=k;j<=n+1;j++){
                    a[i][j]=a[i][j]-d*a[k][j];
                }
            }
        }
    }
}
/*-------------------------------------------------------------------*/