#include <iostream.h>
 #include <fstream.h>  
 #include <math.h>
 #include <stdlib.h>
 #define N 600
 #define M 7

 int brinv(double f[],int n);                                      //矩陣求逆
 int brmul(double a[],double b[],int m,int n,int k,double c[]);    //矩陣相乘

 void main()
 {
     double u[2*N],e[2*N],z[2*N],z1[N];
	 double x1[N*2*M],x1t[2*M*N],x2[N*2],x2t[2*N];
	 double g1[2*M*2*M],g2[2*M*N],g3[2*M],g4[2*2],g5[2*2*M],g6[2*N],g7[2*2];
	 double A[2*M*2],B[2*2],sita1[2*(M-1)],sita2[2],sita[2*M];
	 double f1[2*M*2],f2[2*M*N],f3[N],f4[N],f5[2*(M-1)],f6[2];
	 double w1[2*M*N],s[N],w[N],fpe,ww[1];
	 int i,k,j,k1;
     ofstream fop("data.txt");
     ifstream fip1("M序列.txt");
     ifstream fip2("Gauss.txt");
     for(i=0;i<2*N;i++)
	 {
	    fip1>>u[i];
	    fip2>>e[i];
	 }
	 z[0]=e[0];
     z[1]=1.185*z[0]+1.08*u[0]+e[1];
     z[2]=1.185*z[1]-0.814*z[0]+1.08*u[1]-0.745*u[0]+e[2];
     z[3]=1.185*z[2]-0.814*z[1]+0.518*z[0]+1.08*u[2]-0.745*u[1]+0.475*u[0]+e[3];
     z[4]=1.185*z[3]-0.814*z[2]+0.518*z[1]-0.349*z[0]+1.08*u[3]-0.745*u[2]+0.475*u[1]-0.253*u[0]+e[4];
     for(i=5;i<2*N;i++)
	 {
	    z[i]=1.185*z[i-1]-0.814*z[i-2]+0.518*z[i-3]-0.349*z[i-4]+0.117*z[i-5]+1.08*u[i-1]-0.745*u[i-2]+0.475*u[i-3]-0.253*u[i-4]+0.123*u[i-5]+e[i];
	 }
     
	 i=0;
	 for(k=0;k<N;k++)
	 {
		 x2[2*i]=-z[k];
		 x2[2*i+1]=u[k];
		 i=i+1;
	 }
	 
	 i=0;
     for(j=0;j<N;j++)
	 {
		x2t[i]=-z[j];
		i=i+1;
	 }
     for(j=0;j<N;j++)
	 {
		x2t[i]=u[j];
		i=i+1;
        
	 }
	 //遞推過程
     for(k1=1;k1<M;k1++)
	 {
        i=0;
	    for(k=0;k<N;k++)
		{
	        for(j=0;j<k1;j++)
			{
	           w1[i]=-z[k+k1-1-j];
               i=i+1;
			}
	       for(j=0;j<k1;j++)
		   {
		      w1[i]=u[k+k1-1-j];
		      i=i+1;
		   }
		}
        i=0;
	    for(k=0;k<N;k++)
		{
           z1[k]=z[k+k1];
	       for(j=0;j<k1-1;j++)
		   {
		      x1[2*i]=-z[k1-1+k-j];
		      x1[2*i+1]=u[k1-1+k-j];
              i=i+1;
		   }
		}
     		 
	    i=0;
        for(k=0;k<k1-1;k++)
		{
		    for(j=0;j<N;j++)
			{
			    x1t[i]=-z[j+k1-1-k];
			    i=i+1;
			}
	    for(j=0;j<N;j++)
		{
             x1t[i]=u[j+k1-1-k];
			 i=i+1;
		 }
	 }
	 
      brmul(x1t,x1,2*k1-2,N,2*k1-2,g1);
      brinv(g1,2*k1-2);
      brmul(g1,x1t,2*k1-2,2*k1-2,600,g2);
      brmul(g2,z1,2*k1-2,N,1,g3);
      brmul(x2t,x2,2,N,2,g4);
	  brmul(x2t,x1,2,N,2*k1-2,g5);
      brmul(g5,g2,2,2*k1-2,N,g6);
      brmul(g6,x2,2,N,2,g7);
	  for(i=0;i<4;i++)
		  B[i]=g4[i]-g7[i];
	  brinv(B,2);
      brmul(g2,x2,2*k1-2,N,2,f1);
	  brmul(f1,B,2*k1-2,2,2,A);
      brmul(A,x2t,2*k1-2,2,N,f2);
      brmul(x1,g3,N,2*k1-2,1,f3);
	  for(i=0;i<N;i++)
	  {
		  f4[i]=z1[i]-f3[i];
	  }
      brmul(f2,f4,2*k1-2,N,1,f5);
	  for(i=0;i<2*k1-2;i++)
	  {
		  sita1[i]=g3[i]-f5[i];
	  }
      brmul(x2t,f4,2,N,1,f6);
      brmul(B,f6,2,2,1,sita2);
	  j=0;
	  for(i=0;i<k1-1;i++)
	  {
		  sita[j]=sita1[2*i];
		  j=j+1;
	  }
	  j=j+1;
	  for(i=0;i<k1-1;i++)
	  {
		  sita[j]=sita1[2*i+1];
		  j=j+1;
	  }
	  sita[k1-1]=sita2[0];
	  sita[2*k1-1]=sita2[1];
	  
	brmul(w1,sita,N,2*k1,1,s);
	for(i=0;i<N;i++)
		 w[i]=z1[i]-s[i];
	brmul(w,w,1,N,1,ww);
	fpe=((N+2*k1)*ww[0])/((N-2*k1)*(N-2*k1));
	cout<<"n: "<<k1<<endl;
	fop<<"n: "<<k1<<endl;
	cout<<"判斷階次FPE的值: "<<fpe<<endl;
	fop<<"判斷階次FPE的值: "<<fpe<<endl;
	
	for(i=0;i<2*k1;i++)
	{
		cout<<sita[i]<<"  ";
		fop<<sita[i]<<"  ";
	}
	 cout<<endl;
	 fop<<endl;
     }
  }
		  
  int brinv(double f[],int n)
  { int *is,*js,i,j,k,l,u,v;
    double d,p;
    is=(int *)malloc(n*sizeof(int));
    js=(int *)malloc(n*sizeof(int));
    for (k=0; k<=n-1; k++)
      { d=0.0;
        for (i=k; i<=n-1; i++)
        for (j=k; j<=n-1; j++)
          { l=i*n+j; p=fabs(f[l]);
            if (p>d) { d=p; is[k]=i; js[k]=j;}
          }
        if (d+1.0==1.0)
          { free(is); free(js); cout<<"err**not inv\n";
            return(0);
          }
        if (is[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=is[k]*n+j;
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
        if (js[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+js[k];
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
        l=k*n+k;
        f[l]=1.0/f[l];
        for (j=0; j<=n-1; j++)
          if (j!=k)
            { u=k*n+j; f[u]=f[u]*f[l];}
        for (i=0; i<=n-1; i++)
          if (i!=k)
            for (j=0; j<=n-1; j++)
              if (j!=k)
                { u=i*n+j;
                  f[u]=f[u]-f[i*n+k]*f[k*n+j];
                }
        for (i=0; i<=n-1; i++)
          if (i!=k)
            { u=i*n+k; f[u]=-f[u]*f[l];}
      }
    for (k=n-1; k>=0; k--)
      { if (js[k]!=k)
          for (j=0; j<=n-1; j++)
            { u=k*n+j; v=js[k]*n+j;
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
        if (is[k]!=k)
          for (i=0; i<=n-1; i++)
            { u=i*n+k; v=i*n+is[k];
              p=f[u]; f[u]=f[v]; f[v]=p;
            }
      }
    free(is); free(js);
    return(1);
  }
 int brmul(double a[],double b[],int m,int n,int k,double c[])
    { int i,j,l,u;
    for (i=0; i<=m-1; i++)
    for (j=0; j<=k-1; j++)
      { u=i*k+j; c[u]=0.0;
        for (l=0; l<=n-1; l++)
          c[u]=c[u]+a[i*n+l]*b[l*k+j];
      }
    return 0;
  }