#include "math.h"//a[m+1] 前m個元素返回m-1次擬合多項式的m個系數，最后一個a[m]返回
//擬合多項式pm-1（x）的最大偏差，a[m]為負則說明在跌代過程中參考偏差不再增大其絕對值
//為當前選擇的參考偏差  m為擬合多項式項數，最高次為m-1要求m《n（給定數據點個數）且
//m〈20，若m〉n或m》20，則本函數按m=min[n-1,20]處理
  #include "stdio.h"//
  void chir(x,y,n,a,m)//
  int n,m;
  double x[],y[],a[];
  { int m1,i,j,l,ii,k,im,ix[21];
    double h[21],ha,hh,y1,y2,h1,h2,d,hm;
    for (i=0; i<=m; i++) a[i]=0.0;
    if (m>=n) m=n-1;
    if (m>=20) m=19;
    m1=m+1;
    ha=0.0;
    ix[0]=0; ix[m]=n-1;
    l=(n-1)/m; j=l;
    for (i=1; i<=m-1; i++)
      { ix[i]=j; j=j+l;}
    while (1==1)
      { hh=1.0;
        for (i=0; i<=m; i++)
          { a[i]=y[ix[i]]; h[i]=-hh; hh=-hh;}
        for (j=1; j<=m; j++)
          { ii=m1; y2=a[ii-1]; h2=h[ii-1];
            for (i=j; i<=m; i++)
              { d=x[ix[ii-1]]-x[ix[m1-i-1]];
                y1=a[m-i+j-1];
                h1=h[m-i+j-1];
                a[ii-1]=(y2-y1)/d;
                h[ii-1]=(h2-h1)/d;
                ii=m-i+j; y2=y1; h2=h1;
              }
          }
        hh=-a[m]/h[m];
        for (i=0; i<=m; i++)
          a[i]=a[i]+h[i]*hh;
        for (j=1; j<=m-1; j++)
          { ii=m-j; d=x[ix[ii-1]];
            y2=a[ii-1];
            for (k=m1-j; k<=m; k++)
              { y1=a[k-1]; a[ii-1]=y2-d*y1;
                y2=y1; ii=k;
              }
          }
        hm=fabs(hh);
        if (hm<=ha) { a[m]=-hm; return;}
        a[m]=hm; ha=hm; im=ix[0]; h1=hh;
        j=0;
        for (i=0; i<=n-1; i++)
          { if (i==ix[j])
              { if (j<m) j=j+1;}
            else
              { h2=a[m-1];
                for (k=m-2; k>=0; k--)
                  h2=h2*x[i]+a[k];
                h2=h2-y[i];
                if (fabs(h2)>hm)
                  { hm=fabs(h2); h1=h2; im=i;}
              }
          }
        if (im==ix[0]) return;
        i=0;l=1;
        while (l==1)
          { l=0;
            if (im>=ix[i])
              { i=i+1;
                if (i<=m) l=1;
              }
          }
        if (i>m) i=m;
        if (i==(i/2)*2) h2=-hh;
        else h2=hh;
        if (h1*h2>=0.0) ix[i]=im;
        else
          { if (im<ix[0])
              { for (j=m-1; j>=0; j--)
                  ix[j+1]=ix[j];
                ix[0]=im;
              }
            else
              { if (im>ix[m])
                  { for (j=1; j<=m; j++)
                      ix[j-1]=ix[j];
                    ix[m]=im;
                  }
                else ix[i-1]=im;
              }
          }
      }
  }


  main()
  { int i;
    double x[101],y[101],a[7];
    for (i=0; i<=100; i++)
      { x[i]=-1.0+0.02*i;
        y[i]=atan(x[i]);
      }
    chir(x,y,101,a,6);
    printf("\n");
    for (i=0; i<=5; i++)
      printf("a(%2d)=%e\n",i,a[i]);
    printf("\n");
    printf("MAX(p-f)=%e\n",a[6]);
    printf("\n");
  }