#define TestNow 0
//一维多项式计算
double plyv(double a[],int n,double x)
{ int i;
double u;
u=a[n-1];
for (i=n-2; i>=0; i--)
u=u*x+a[i];
return(u);
}
//二维多项式计算
double bply(double a[6][6],int n,double x,double y)
{ int i,j;
double u,v;
v = a[0][5];
for(j=n-2;j>=0;j--)
{
u=a[n-1-j][j];
for (i=n-2-j; i>=0; i--)
{
u=u*x+a[i][j];
}
v = v*y+u;
}
return(v);
}
//针对二维5阶多项式的优化函数
double bply55(double xin,double yin/*,double a[21]*/)
{
double x=xin;
double y=yin;
double z = 0;
/*
double P00 = a[0];
double P10 = a[1];
double P01 = a[2];
double P20 = a[3];
double P11 = a[4];
double P02 = a[5];
double P30 = a[6];
double P21 = a[7];
double P12 = a[8];
double P03 = a[9];
double P40 = a[10];
double P31 = a[11];
double P22 = a[12];
double P13 = a[13];
double P04 = a[14];
double P50 = a[15];
double P41 = a[16];
double P32 = a[17];
double P23 = a[18];
double P14 = a[19];
double P05 = a[20];
*/
/*
double P00 = -157.6646;
double P10 = 107.4548;
double P01 = -137.6695;
double P20 = -49.1798;
double P11 = 61.5766;
double P02 = -15.4219;
double P30 = 0.039778;
double P21 = -13.2057;
double P12 = 33.0038;
double P03 = -7.0569;
double P40 = 5.0611;
double P31 = 5.0977;
double P22 = -16.3451;
double P13 = 1.6553;
double P04 = 0.0029376;
double P50 = -0.91358;
double P41 = -1.6222;
double P32 = 2.6572;
double P23 = -0.079645;
double P14 = 0.053844;
double P05 = 0.022143;
*/
/********************************************************************************
觉得用宏这种方式替代系数最好
**********************************************************************************/
#define P00 -157.6646#define P10 107.4548
#define P01 -137.6695
#define P20 -49.1798
#define P11 61.5766
#define P02 -15.4219
#define P30 0.039778
#define P21 -13.2057
#define P12 33.0038
#define P03 -7.0569
#define P40 5.0611
#define P31 5.0977
#define P22 -16.3451
#define P13 1.6553
#define P04 0.0029376
#define P50 -0.91358
#define P41 -1.6222
#define P32 2.6572
#define P23 -0.079645
#define P14 0.053844
#define P05 0.022143
z = P00+x*(P10+x*(P20+x*(P30+x*(P40+x*P50))))
+ y*( (P01+x*(P11+x*(P21+x*(P31+x*P41))))
+ y*( (P02+x*(P12+x*(P22+x*P32)))
+ y*( (P03+x*(P13+x*P23))
+ y*( (P04+x*P14)
+ y*P05))));
return z;
}
#if TestNow
void plyv0()
{
int i;
static double a[7]={-20.0,7.0,-7.0,1.0,3.0,-5.0,2.0};
static double x[6]={0.9,-0.9,1.1,-1.1,1.3,-1.3};
printf("\n");
for (i=0; i<=5; i++)
printf("x(%d)=%5.2lf p(%d)=%13.7e\n",
i,x[i],i,plyv(a,7,x[i]));
printf("\n");
}
void bply0(double xin,double yin,double a[21])
{
double x=xin;
double y=yin;
double z = 0;
// double a[21]={-157.6646,107.4548,-137.6695,-49.1798,61.5766,-15.4219,0.039778,-13.2057,33.0038,-7.0569,5.0611,5.0977,-16.3451,1.6553,0.0029376,-0.91358,-1.6222,2.6572,-0.079645,0.053844,0.022143};
double b[6][6]={0};
printf("x=%5.2lf, y=%5.2lf z=%10.4lf\n",x,y,z);
b[0][0] = a[0];
b[1][0] = a[1];
b[0][1] = a[2];
b[2][0] = a[3];
b[1][1] = a[4];
b[0][2] = a[5];
b[3][0] = a[6];
b[2][1] = a[7];
b[1][2] = a[8];
b[0][3] = a[9];
b[4][0] = a[10];
b[3][1] = a[11];
b[2][2] = a[12];
b[1][3] = a[13];
b[0][4] = a[14];
b[5][0] = a[15];
b[4][1] = a[16];
b[3][2] = a[17];
b[2][3] = a[18];
b[1][4] = a[19];
b[0][5] = a[20];
z = bply(b,6,x,y);
printf("\n");
printf("x=%5.2lf, y=%5.2lf z=%10.4lf\n",x,y,z);
printf("\n");
}
void main()
{
double x=0.233333333;
double y=500;
double a[21]={-157.6646,107.4548,-137.6695,-49.1798,61.5766,-15.4219,0.039778,-13.2057,33.0038,-7.0569,5.0611,5.0977,-16.3451,1.6553,0.0029376,-0.91358,-1.6222,2.6572,-0.079645,0.053844,0.022143};
x=(x-(0.076924))/(0.034819);
y=(y-(97.2882))/(85.3981);
bply0(x,y,a);
bply55(x,y,a);
}
#endif