// 最佳一致逼近多項式.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "iostream.h"

  
void bua(double a,double b,double p[],int n,double e,double (*f)(double))
{ 
	int i,j,k,m;
    double x[21],g[21],d,t,u,s,xx,x0,h,yy;
    if (n>20) n=20;		//最高逼近19次
    m=n+1; d=1.0e+35;
    for (k=0; k<=n; k++)	//初始點采集n+1階切比雪夫多項式的交錯點組
      { 
		t=cos((n-k)*3.1415926/(1.0*n));
        x[k]=(b+a+(b-a)*t)/2.0;
      }
    while (1==1)	//迭代
      { u=1.0;
        for (i=0; i<=m-1; i++)		//計算f(x)及(-1)u
          { p[i]=(*f)(x[i]);
            g[i]=-u;
			u=-u;
          }
        for (j=0; j<=n-1; j++)		//確定參考偏差u
          { 
			k=m;
			s=p[k-1]; 
			xx=g[k-1];
            for (i=j; i<=n-1; i++)
              { 
				t=p[n-i+j-1];
				x0=g[n-i+j-1];
                p[k-1]=(s-t)/(x[k-1]-x[m-i-2]);
                g[k-1]=(xx-x0)/(x[k-1]-x[m-i-2]);
                k=n-i+j;
				s=t;
				xx=x0;
              }
          }
        u=-p[m-1]/g[m-1];
        for (i=0; i<=m-1; i++)		//牛頓插值公式確定P系數
          p[i]=p[i]+g[i]*u;
        for (j=1; j<=n-1; j++)
          { 
			k=n-j; 
			h=x[k-1]; 
			s=p[k-1];
            for (i=m-j; i<=n; i++)
              { 
				t=p[i-1]; 
				p[k-1]=s-h*t;
                s=t; 
				k=i;
              }
          }
        p[m-1]=fabs(u); 
		u=p[m-1];
        if (fabs(u-d)<=e) return;		//若精度滿足，返回
        d=u; 
		h=0.1*(b-a)/(1.0*n);
        xx=a; 
		x0=a;
        while (x0<=b)	//掃描最大值點xx
          {
			s=(*f)(x0);
			t=p[n-1];
            for (i=n-2; i>=0; i--)
              t=t*x0+p[i];
            s=fabs(s-t);
            if (s>u) 
				{ 
					u=s; 
					xx=x0;
				}
            x0=x0+h;
          }
        s=(*f)(xx);		//計算yy=f(x)-P(x)
		t=p[n-1];
        for (i=n-2; i>=0; i--)
          t=t*xx+p[i];
        yy=s-t;
		i=1;
		j=n+1;
        while ((j-i)!=1)	//尋找最大值點應處在點集中的位置
          { 
			k=(i+j)/2;
            if (xx<x[k-1]) j=k;
            else i=k;
          }
        if (xx<x[0])	//最大值點xx應在a與x0之間
          { 
			s=(*f)(x[0]);
			t=p[n-1];
            for (k=n-2; k>=0; k--)	//計算f(x0)-P(x0)
              t=t*x[0]+p[k];
            s=s-t;
            if (s*yy>0.0) x[0]=xx;	//若f(x)-P(x)與f(x0)-P(x0)同號,用xx代替x0
            else
              {
				for (k=n-1; k>=0; k--)
                  x[k+1]=x[k];
                x[0]=xx;
              }
          }
        else
          { 
			if (xx>x[n])	//最大值點xx應在xn+1與b之間
              {
				s=(*f)(x[n]); t=p[n-1];
                for (k=n-2; k>=0; k--)	//計算f(xn+1)-P(xn+1)
                  t=t*x[n]+p[k];
                s=s-t;
                if (s*yy>0.0) x[n]=xx;	//若f(x)-P(x)與f(x0)-P(x0)同號,用xx代替xn+1
                else
                  { 
					for (k=0; k<=n-1; k++)
                      x[k]=x[k+1];
                    x[n]=xx;
                  }
              }
            else	//最大值點xx應在xk+1與xk之間
              { 
				i=i-1;
				j=j-1;
                s=(*f)(x[i]);
				t=p[n-1];
                for (k=n-2; k>=0; k--)	//計算f(xk)-P(xk)
                  t=t*x[i]+p[k];
                s=s-t;
                if (s*yy>0.0) x[i]=xx;	//若f(x)-P(x)與f(x0)-P(x0)同號,用xx代替xk,否則用xx代替xk+1
                else x[j]=xx;
              }
          }
      }
}

double f(double x)
{ 
	double y;
    y=sqrt(x);	//待逼函數
    return(y);
}

void main()
{ 
	cout << "  *******************************************" << endl;
	cout << "  **                                       **" << endl;
	cout << "  **          最佳一致逼近多項式           **" << endl;
	cout << "  **                                       **" << endl;
	cout << "  *******************************************" << endl << endl;
	double a,b,e,p[5],f(double);
    a=0.25;		//逼近區間下限
	b=1.0;		//逼近區間上限
	e=1.0e-10;	//逼近精度
    bua(a,b,p,2,e,f);
    cout << "所求一次最佳一致逼近多項式為:" << endl;
	cout << "P1(x)=" << p[1] << "x+" << p[0] << endl;
 }