1 // 使用切比雪夫算法求多项式  2 void YcChebyshevFitSpline::CalculateMultinomialValues(const void* valuesPtr, int stride, int n, int m, float* a) const  3 {  4     memset(a, 0, sizeof(float)*m);  5   6     float xStep = 1.0f/(n - 1);  7   8     int m1,i,j,l,ii,k,im,ix[21];  9     float h[21],ha,hh,y1,y2,h1,h2,d,hm; 10  11     for (i=0; i<=m; i++)  12     { 13         a[i]=0.0; 14     } 15     if (m>=n)  16     { 17         m=n-1; 18     } 19     if (m>=20)  20     { 21         m=19; 22     } 23  24     m1=m+1; 25     ha=0.0f; 26     ix[0]=0;  27     ix[m]=n-1; 28     l=(n-1)/m;  29     j=l; 30     for (i=1; i<=m-1; i++) 31     {  32         ix[i]=j; j=j+l; 33     } 34  35     while (1==1) 36     {  37         hh=1.0f; 38         for (i=0; i<=m; i++) 39         {  40             a[i]=YfGetFloatValue(valuesPtr, stride, ix[i]);  41             h[i]=-hh;  42             hh=-hh; 43         } 44          45         for (j=1; j<=m; j++) 46         {  47             ii=m1;  48             y2=a[ii-1];  49             h2=h[ii-1]; 50             for (i=j; i<=m; i++) 51             {  52                 d=xStep*ix[ii-1] - xStep*ix[m1-i-1]; 53                 y1=a[m-i+j-1]; 54                 h1=h[m-i+j-1]; 55                 a[ii-1]=(y2-y1)/d; 56                 h[ii-1]=(h2-h1)/d; 57                 ii=m-i+j; y2=y1; h2=h1; 58             } 59         } 60  61         hh=-a[m]/h[m]; 62         for (i=0; i<=m; i++) 63         { 64             a[i]=a[i]+h[i]*hh; 65         } 66          67         for (j=1; j<=m-1; j++) 68         {  69             ii=m-j;  70             d=xStep*ix[ii-1]; 71             y2=a[ii-1]; 72             for (k=m1-j; k<=m; k++) 73             {  74                 y1=a[k-1];  75                 a[ii-1]=y2-d*y1; 76                 y2=y1; ii=k; 77             } 78         } 79  80         hm=fabs(hh); 81         if (hm<=ha)  82         {  83             a[m]=-hm;  84             return; 85         } 86  87         a[m]=hm; ha=hm; im=ix[0]; h1=hh; 88         j=0; 89         for (i=0; i<=n-1; i++) 90         {  91             if (i==ix[j]) 92             {  93                 if (j<m) j=j+1; 94             } 95             else 96             {  97                 h2=a[m-1]; 98                 for (k=m-2; k>=0; k--) 99                 {100                     h2=h2*xStep*i+a[k];101                 }102                 h2=h2-YfGetFloatValue(valuesPtr, stride, i);103                 if (fabs(h2)>hm)104                 { 105                     hm=fabs(h2); 106                     h1=h2; 107                     im=i;108                 }109             }110         }111         if (im==ix[0]) 112         {113             return;114         }115 116         i=0;117         l=1;118         while (l==1)119         { 120             l=0;121             if (im>=ix[i])122             { 123                 i=i+1;124                 if (i<=m) 125                 {126                     l=1;127                 }128             }129         }130 131         if (i>m) 132         {133             i=m;134         }135         if (i==(i/2)*2) 136         {137             h2=-hh;138         }139         else 140         {141             h2=hh;142         }143 144         if (h1*h2>=0.0f) 145         {146             ix[i]=im;147         }148         else149         { 150             if (im<ix[0])151             { 152                 for (j=m-1; j>=0; j--)153                 {154                     ix[j+1]=ix[j];155                 }156                 ix[0]=im;157             }158             else159             { 160                 if (im>ix[m])161                 { 162                     for (j=1; j<=m; j++)163                     {164                         ix[j-1]=ix[j];165                     }166                     ix[m]=im;167                 }168                 else 169                 {170                     ix[i-1]=im;171                 }172             }173         }174     }175 }