//@—ĢŽUŒ^ƒf[ƒ^‚ɑ΂·‚éÅŽ‚Qæ‹ßŽ—‘―€ŽŪ 

#include<stdio.h>

#include<math.h>

#define  M        6      //  ÃÞ-Ā É š―ģ 

#define  N        3      //  Āšģž· É žÞ―ģ +1

#define  EPSILON  1.0E-5

#define  TRUE     1

#define  FALSE    0

double  data[M][2]={{ 0.0, 10.0 },

                   { 1.0, 9.05 },

                   { 2.0, 8.19 },

                   { 3.0, 7.41 },

                   { 4.0, 6.70 },

                   { 5.0, 6.07 }};

double  a[N][N+1];

double  c[N];

int sw=TRUE;

void main()

{

    void op1();

    void comp_a();

    void solve_eq();

    void op2();

    op1();

    comp_a();

    solve_eq();

    op2();

}

void op1()

{

    int  k;

    printf("\n\nŠÖ”f(x)‚Ė—ĢŽUŒ^ƒf[ƒ^‚ɑ΂·‚éÅŽ‚Qæ‹ßŽ—‘―€ŽŪg(x)‚ð‹‚ß‚é\n\n");

    printf(" iƒf[ƒ^”%dA‘―€ŽŪ‚ĖŽŸ”%dA‚Ėę‡j\n\n",M,N);

    printf("data     x        fx\n");

    for(k=0;k<M;k++)

        printf(" #%2d  %f %f\n",k+1,data[k][0],data[k][1]);

}

void comp_a()

{

    int   i,j,k;

    double powf(double,int);

    for(i=0;i<N;i++)

    {

        for(j=0;j<N;j++)

        {

            if(i==0&&j==0)

                a[i][j]=M;

            else

            {

                a[i][j]=0.0;

                for(k=0;k<M;k++)

                    a[i][j]=a[i][j]+powf(data[k][0],i+j);

            }

        }

        a[i][N]=0.0;

        for(k=0;k<M;k++)

            if(i==0)

                a[i][N]=a[i][N]+data[k][1];

            else

                a[i][N]=a[i][N]+powf(data[k][0],i)*data[k][1];

    }

    printf("\na[i][j]=\n");

    for(i=0;i<N;i++)

    {

        for(j=0;j<N;j++)

            printf(" %10.6f",a[i][j]);

        printf("    %10.6f\n",a[i][N]);

    }

}

double powf(double x,int n)

{

    int   i;

    double y=1.0;

    if(n==0)

        if(x!=0)

            return 1;

        else

            return 0;

    else

    {

        for(i=0;i<n;i++)

            y=y*x;

        return y;

    }

}

void solve_eq()

{

    int   i,j,k,ik;

    double ak,aik,w;

    for(k=0;k<N;k++)

    {

        if(fabs(a[k][k])<=EPSILON)

        {

            ik=k+1;

            while((ik<N) && (fabs(a[ik][k])<EPSILON))

                ik++;

            if(ik<N)

                for(j=k;j<N+1;j++)

                {

                    w=a[k][j];

                    a[k][j]=a[ik][j];

                    a[ik][j]=w;

                }

            else

            {

                sw=FALSE;

                return;

            }

        }

        ak=a[k][k];

        for(j=k;j<=N;j++)

            a[k][j]=a[k][j]/ak;

        for(i=0;i<N;i++)

            if(i!=k)

            {

                aik=a[i][k];

                for(j=k;j<=N;j++)

                    a[i][j]=a[i][j]-aik*a[k][j];

            }

    }

}

void op2()

{

    int i;

    if(sw==TRUE)

    {

        printf("\ng(x)=");

        for(i=0;i<N-1;i++)

            printf("(%f)x^%d +",a[i][N],i);

        printf("(%f)x^%d\n\n",a[N-1][N],N-1);

    }

    else

        printf("Impossible.\n\n");

}