# source:trunk/src/lol/math/remez.h@1117Tweet

Last change on this file since 1117 was 1117, checked in by sam, 11 years ago

math: move the Remez algorithm implementation to the core.

• Property svn:keywords set to `Id`
File size: 8.9 KB
Line
1//
2// Lol Engine
3//
4// Copyright: (c) 2010-2011 Sam Hocevar <sam@hocevar.net>
5//   This program is free software; you can redistribute it and/or
6//   modify it under the terms of the Do What The Fuck You Want To
8//   http://sam.zoy.org/projects/COPYING.WTFPL for more details.
9//
10
11//
12// The Remez class
13// ---------------
14//
15
16#if !defined __LOL_MATH_REMEZ_H__
17#define __LOL_MATH_REMEZ_H__
18
19namespace lol
20{
21
22template<int ORDER, typename T> class RemezSolver
23{
24public:
25    typedef T RealFunc(T const &x);
26
27    RemezSolver()
28    {
29    }
30
31    void Run(T a, T b, RealFunc *func, RealFunc *weight, int steps)
32    {
33        m_func = func;
34        m_weight = weight;
35        m_k1 = (b + a) >> 1;
36        m_k2 = (b - a) >> 1;
37        m_invk2 = re(m_k2);
38        m_invk1 = -m_k1 * m_invk2;
39
40        Init();
41
42        PrintPoly();
43
44        for (int n = 0; n < steps; n++)
45        {
46            FindExtrema();
47            Step();
48
49            PrintPoly();
50
51            FindZeroes();
52        }
53
54        FindExtrema();
55        Step();
56
57        PrintPoly();
58    }
59
60    T ChebyEval(T const &x)
61    {
62        T ret = 0.0, xn = 1.0;
63
64        for (int i = 0; i < ORDER + 1; i++)
65        {
66            T mul = 0;
67            for (int j = 0; j < ORDER + 1; j++)
68                mul += coeff[j] * (T)Cheby(j, i);
69            ret += mul * xn;
70            xn *= x;
71        }
72
73        return ret;
74    }
75
76    void Init()
77    {
78        /* Pick up x_i where error will be 0 and compute f(x_i) */
79        T fxn[ORDER + 1];
80        for (int i = 0; i < ORDER + 1; i++)
81        {
82            zeroes[i] = (T)(2 * i - ORDER) / (T)(ORDER + 1);
83            fxn[i] = Value(zeroes[i]);
84        }
85
86        /* We build a matrix of Chebishev evaluations: row i contains the
87         * evaluations of x_i for polynomial order n = 0, 1, ... */
88        lol::Mat<ORDER + 1, T> mat;
89        for (int i = 0; i < ORDER + 1; i++)
90        {
91            /* Compute the powers of x_i */
92            T powers[ORDER + 1];
93            powers[0] = 1.0;
94            for (int n = 1; n < ORDER + 1; n++)
95                 powers[n] = powers[n - 1] * zeroes[i];
96
97            /* Compute the Chebishev evaluations at x_i */
98            for (int n = 0; n < ORDER + 1; n++)
99            {
100                T sum = 0.0;
101                for (int k = 0; k < ORDER + 1; k++)
102                    sum += (T)Cheby(n, k) * powers[k];
103                mat.m[i][n] = sum;
104            }
105        }
106
107        /* Solve the system */
108        mat = mat.inv();
109
110        /* Compute interpolation coefficients */
111        for (int j = 0; j < ORDER + 1; j++)
112        {
113            coeff[j] = 0;
114            for (int i = 0; i < ORDER + 1; i++)
115                coeff[j] += mat.m[j][i] * fxn[i];
116        }
117    }
118
119    void FindZeroes()
120    {
121        /* Find ORDER + 1 zeroes of the error function. No need to
122         * compute the relative error: its zeroes are at the same
123         * place as the absolute error! */
124        for (int i = 0; i < ORDER + 1; i++)
125        {
126            struct { T value, error; } left, right, mid;
127
128            left.value = control[i];
129            left.error = ChebyEval(left.value) - Value(left.value);
130            right.value = control[i + 1];
131            right.error = ChebyEval(right.value) - Value(right.value);
132
133            static T limit = (T)1 >> 500;
134            static T zero = (T)0;
135            while (fabs(left.value - right.value) > limit)
136            {
137                mid.value = (left.value + right.value) >> 1;
138                mid.error = ChebyEval(mid.value) - Value(mid.value);
139
140                if ((left.error < zero && mid.error < zero)
141                     || (left.error > zero && mid.error > zero))
142                    left = mid;
143                else
144                    right = mid;
145            }
146
147            zeroes[i] = mid.value;
148        }
149    }
150
151    void FindExtrema()
152    {
153        using std::printf;
154
155        /* Find ORDER + 2 extrema of the error function. We need to
156         * compute the relative error, since its extrema are at slightly
157         * different locations than the absolute error’s. */
158        T final = 0;
159
160        for (int i = 0; i < ORDER + 2; i++)
161        {
162            T a = -1, b = 1;
163            if (i > 0)
164                a = zeroes[i - 1];
165            if (i < ORDER + 1)
166                b = zeroes[i];
167
168            for (;;)
169            {
170                T c = a, delta = (b - a) >> 3;
171                T maxerror = 0;
172                T maxweight = 0;
173                int best = -1;
174                for (int k = 1; k <= 7; k++)
175                {
176                    T error = ChebyEval(c) - Value(c);
177                    T weight = Weight(c);
178                    if (fabs(error * maxweight) >= fabs(maxerror * weight))
179                    {
180                        maxerror = error;
181                        maxweight = weight;
182                        best = k;
183                    }
184                    c += delta;
185                }
186
187                b = a + (T)(best + 1) * delta;
188                a = a + (T)(best - 1) * delta;
189
190                if (b - a < (T)1e-18)
191                {
192                    T e = maxerror / maxweight;
193                    if (e > final)
194                        final = e;
195                    control[i] = (a + b) >> 1;
196                    break;
197                }
198            }
199        }
200
201        printf("Final error: ");
202        final.print(40);
203    }
204
205    void Step()
206    {
207        /* Pick up x_i where error will be 0 and compute f(x_i) */
208        T fxn[ORDER + 2];
209        for (int i = 0; i < ORDER + 2; i++)
210            fxn[i] = Value(control[i]);
211
212        /* We build a matrix of Chebishev evaluations: row i contains the
213         * evaluations of x_i for polynomial order n = 0, 1, ... */
214        lol::Mat<ORDER + 2, T> mat;
215        for (int i = 0; i < ORDER + 2; i++)
216        {
217            /* Compute the powers of x_i */
218            T powers[ORDER + 1];
219            powers[0] = 1.0;
220            for (int n = 1; n < ORDER + 1; n++)
221                 powers[n] = powers[n - 1] * control[i];
222
223            /* Compute the Chebishev evaluations at x_i */
224            for (int n = 0; n < ORDER + 1; n++)
225            {
226                T sum = 0.0;
227                for (int k = 0; k < ORDER + 1; k++)
228                    sum += (T)Cheby(n, k) * powers[k];
229                mat.m[i][n] = sum;
230            }
231            if (i & 1)
232                mat.m[i][ORDER + 1] = fabs(Weight(control[i]));
233            else
234                mat.m[i][ORDER + 1] = -fabs(Weight(control[i]));
235        }
236
237        /* Solve the system */
238        mat = mat.inv();
239
240        /* Compute interpolation coefficients */
241        for (int j = 0; j < ORDER + 1; j++)
242        {
243            coeff[j] = 0;
244            for (int i = 0; i < ORDER + 2; i++)
245                coeff[j] += mat.m[j][i] * fxn[i];
246        }
247
248        /* Compute the error */
249        T error = 0;
250        for (int i = 0; i < ORDER + 2; i++)
251            error += mat.m[ORDER + 1][i] * fxn[i];
252    }
253
254    int Cheby(int n, int k)
255    {
256        if (k > n || k < 0)
257            return 0;
258        if (n <= 1)
259            return (n ^ k ^ 1) & 1;
260        return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k);
261    }
262
263    int Comb(int n, int k)
264    {
265        if (k == 0 || k == n)
266            return 1;
267        return Comb(n - 1, k - 1) + Comb(n - 1, k);
268    }
269
270    void PrintPoly()
271    {
272        using std::printf;
273
274        /* Transform Chebyshev polynomial weights into powers of X^i
275         * in the [-1..1] range. */
276        T bn[ORDER + 1];
277
278        for (int i = 0; i < ORDER + 1; i++)
279        {
280            bn[i] = 0;
281            for (int j = 0; j < ORDER + 1; j++)
282                bn[i] += coeff[j] * (T)Cheby(j, i);
283        }
284
285        /* Transform a polynomial in the [-1..1] range into a polynomial
286         * in the [a..b] range. */
287        T k1p[ORDER + 1], k2p[ORDER + 1];
288        T an[ORDER + 1];
289
290        for (int i = 0; i < ORDER + 1; i++)
291        {
292            k1p[i] = i ? k1p[i - 1] * m_invk1 : (T)1;
293            k2p[i] = i ? k2p[i - 1] * m_invk2 : (T)1;
294        }
295
296        for (int i = 0; i < ORDER + 1; i++)
297        {
298            an[i] = 0;
299            for (int j = i; j < ORDER + 1; j++)
300                an[i] += (T)Comb(j, i) * k1p[j - i] * bn[j];
301            an[i] *= k2p[i];
302        }
303
304        printf("Polynomial estimate:\n");
305        for (int j = 0; j < ORDER + 1; j++)
306        {
307            if (j)
308                printf("+");
309            printf("x**%i*", j);
310            an[j].print(40);
311        }
312        printf("\n");
313    }
314
315    T Value(T const &x)
316    {
317        return m_func(x * m_k2 + m_k1);
318    }
319
320    T Weight(T const &x)
321    {
322        return m_weight(x * m_k2 + m_k1);
323    }
324
325    /* ORDER + 1 Chebyshev coefficients and 1 error value */
326    T coeff[ORDER + 2];
327    /* ORDER + 1 zeroes of the error function */
328    T zeroes[ORDER + 1];
329    /* ORDER + 2 control points */
330    T control[ORDER + 2];
331
332private:
333    RealFunc *m_func, *m_weight;
334    T m_k1, m_k2, m_invk1, m_invk2;
335};
336
337} /* namespace lol */
338
339#endif /* __LOL_MATH_REMEZ_H__ */
340
Note: See TracBrowser for help on using the repository browser.