source: trunk/test/math/remez-solver.h @ 1014

Last change on this file since 1014 was 1014, checked in by sam, 9 years ago

test: various improvements to the Remez exchange solver.

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