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

核心代码：

  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