1 // 使用最小二乘算法求多项式  2 void YcLeastSquaresFitSpline::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 i,j,k;  9     float z,p,c,g,q,d1,d2,s[20],t[20],b[20]; 10     for (i=0; i<=m-1; i++)  11     { 12         a[i]=0.0f; 13     } 14  15     if (m>n)  16     { 17         m=n; 18     } 19     if (m>20)  20     { 21         m=20; 22     } 23  24     z=0.0f; 25     for (i=0; i<=n-1; i++)  26     { 27         z=z+xStep*i/(1.0f*n); 28     } 29  30     b[0]=1.0f;  31     d1=1.0f*n;  32     p=0.0f;  33     c=0.0f; 34     for (i=0; i<=n-1; i++) 35     {  36         p=p+(xStep*i-z);  37         c=c+YfGetFloatValue(valuesPtr, stride, i);  38     } 39     c=c/d1;  40     p=p/d1; 41     a[0]=c*b[0]; 42     if (m>1) 43     {  44         t[1]=1.0f;  45         t[0]=-p; 46         d2=0.0f;  47         c=0.0f;  48         g=0.0f; 49  50         for (i=0; i<=n-1; i++) 51         {  52             q=xStep*i-z-p;  53             d2=d2+q*q; 54             c=c+YfGetFloatValue(valuesPtr, stride, i)*q; 55             g=g+(xStep*i-z)*q*q; 56         } 57         c=c/d2;  58         p=g/d2;  59         q=d2/d1; 60         d1=d2; 61         a[1]=c*t[1];  62         a[0]=c*t[0]+a[0]; 63     } 64  65     for (j=2; j<=m-1; j++) 66     {  67         s[j]=t[j-1]; 68         s[j-1]=-p*t[j-1]+t[j-2]; 69         if (j>=3) 70         { 71             for (k=j-2; k>=1; k--) 72             { 73                 s[k]=-p*t[k]+t[k-1]-q*b[k]; 74             } 75         } 76  77         s[0]=-p*t[0]-q*b[0]; 78         d2=0.0f;  79         c=0.0f;  80         g=0.0f; 81         for (i=0; i<=n-1; i++) 82         {  83             q=s[j]; 84             for (k=j-1; k>=0; k--) 85             { 86                 q=q*(xStep*i-z)+s[k]; 87             } 88             d2=d2+q*q;  89             c=c+YfGetFloatValue(valuesPtr, stride, i)*q; 90             g=g+(xStep*i-z)*q*q; 91         } 92  93         c=c/d2;  94         p=g/d2;  95         q=d2/d1; 96         d1=d2; 97         a[j]=c*s[j];  98         t[j]=s[j]; 99         for (k=j-1; k>=0; k--)100         { 101             a[k]=c*s[k]+a[k];102             b[k]=t[k]; 103             t[k]=s[k];104         }105     }106 }