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 }