样条之连分式插值函数

核心代码：

 1 //////////////////////////////////////////////////////////////////////
 2 // 连分式等距插值
 3 //////////////////////////////////////////////////////////////////////
 4 static float GetValuePqs(const void* valuesPtr, int stride, int n, float t)
 5 { 
 6     int i,j,k,m, l;
 7     float z,hh,xi,xj;
 8     float b[8];
 9 
10     // 初值
11     z = 0.0;
12 
13     // 特例处理
14     if (n < 1) 
15     {
16         return(z);
17     }
18     if (n == 1) 
19     { 
20         z = YfGetFloatValue(valuesPtr, stride, 0);  
21         return(z);
22     }
23     if (n == 2)
24     { 
25         float y0 = YfGetFloatValue(valuesPtr, stride, 0); 
26         float y1 = YfGetFloatValue(valuesPtr, stride, 1); 
27         z = y0 + (y1 - y0)*t;
28         return(z);
29     }
30 
31     float xStep = 1.0f/(n - 1);
32 
33     // 连分式插值
34     if (n <= 8) 
35     { 
36         k = 0; 
37         m = n;
38     }
39     else if (t < (4.0f*xStep)) 
40     { 
41         k=0; 
42         m=8;
43     }
44     else if (t > ((n-5)*xStep)) 
45     { 
46         k = n-8; 
47         m = 8;
48     }
49     else 
50     { 
51         k = (int)(t/xStep)-3; 
52         m = 8;
53     }
54 
55     b[0] = YfGetFloatValue(valuesPtr, stride, k); 
56     for (i = 2; i <= m; i++)
57     { 
58         hh = YfGetFloatValue(valuesPtr, stride, i+k-1); 
59         l = 0; 
60         j = 1;
61 
62         while ((l == 0) && (j <= i-1))
63         { 
64             if (fabs(hh-b[j-1])+1.0f == 1.0f) 
65             {
66                 l = 1;
67             }
68             else
69             { 
70                 xi = (i+k-1)*xStep;
71                 xj = (j+k-1)*xStep;
72                 hh = (xi-xj)/(hh-b[j-1]);
73             }
74 
75             j = j+1;
76         }
77 
78         b[i-1]=hh;
79         if (l != 0) 
80         {
81             b[i-1] = 1.0e+35F;
82         }
83     }
84 
85     z = b[m-1];
86     for (i = m-1; i >= 1; i--)
87     {
88         z = b[i-1]+(t-(i+k-1)*xStep)/z;
89     }
90 
91     return(z);
92 }