样条之切比雪夫算法求多项式

核心代码：
  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         else
149         { 
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             else
159             { 
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 }
切图：
相关软件的下载地址为：http://files.cnblogs.com/WhyEngine/TestSpline.zip