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

Last change on this file since 1050 was 1027, checked in by sam, 10 years ago

test: minor Remez algorithm tests.

  • Property svn:keywords set to Id
File size: 8.9 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            for (;;)
158            {
159                real c = a, delta = (b - a) >> 3;
160                real maxerror = 0;
161                real maxweight = 0;
162                int best = -1;
163                for (int k = 1; k <= 7; k++)
164                {
165                    real error = ChebyEval(c) - Value(c);
166                    real weight = Weight(c);
167                    if (fabs(error * maxweight) >= fabs(maxerror * weight))
168                    {
169                        maxerror = error;
170                        maxweight = weight;
171                        best = k;
172                    }
173                    c += delta;
174                }
175
176                b = a + (real)(best + 1) * delta;
177                a = a + (real)(best - 1) * delta;
178
179                if (b - a < (real)1e-18)
180                {
181                    real e = maxerror / maxweight;
182                    if (e > final)
183                        final = e;
184                    control[i] = (a + b) >> 1;
185                    break;
186                }
187            }
188        }
189
190        printf("Final error: ");
191        final.print(40);
192    }
193
194    void Step()
195    {
196        /* Pick up x_i where error will be 0 and compute f(x_i) */
197        real fxn[ORDER + 2];
198        for (int i = 0; i < ORDER + 2; i++)
199            fxn[i] = Value(control[i]);
200
201        /* We build a matrix of Chebishev evaluations: row i contains the
202         * evaluations of x_i for polynomial order n = 0, 1, ... */
203        Matrix<ORDER + 2> mat;
204        for (int i = 0; i < ORDER + 2; i++)
205        {
206            /* Compute the powers of x_i */
207            real powers[ORDER + 1];
208            powers[0] = 1.0;
209            for (int n = 1; n < ORDER + 1; n++)
210                 powers[n] = powers[n - 1] * control[i];
211
212            /* Compute the Chebishev evaluations at x_i */
213            for (int n = 0; n < ORDER + 1; n++)
214            {
215                real sum = 0.0;
216                for (int k = 0; k < ORDER + 1; k++)
217                    sum += (real)Cheby(n, k) * powers[k];
218                mat.m[i][n] = sum;
219            }
220            if (i & 1)
221                mat.m[i][ORDER + 1] = fabs(Weight(control[i]));
222            else
223                mat.m[i][ORDER + 1] = -fabs(Weight(control[i]));
224        }
225
226        /* Solve the system */
227        mat = mat.inv();
228
229        /* Compute interpolation coefficients */
230        for (int j = 0; j < ORDER + 1; j++)
231        {
232            coeff[j] = 0;
233            for (int i = 0; i < ORDER + 2; i++)
234                coeff[j] += mat.m[j][i] * fxn[i];
235        }
236
237        /* Compute the error */
238        real error = 0;
239        for (int i = 0; i < ORDER + 2; i++)
240            error += mat.m[ORDER + 1][i] * fxn[i];
241    }
242
243    int Cheby(int n, int k)
244    {
245        if (k > n || k < 0)
246            return 0;
247        if (n <= 1)
248            return (n ^ k ^ 1) & 1;
249        return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k);
250    }
251
252    int Comb(int n, int k)
253    {
254        if (k == 0 || k == n)
255            return 1;
256        return Comb(n - 1, k - 1) + Comb(n - 1, k);
257    }
258
259    void PrintPoly()
260    {
261        /* Transform Chebyshev polynomial weights into powers of X^i
262         * in the [-1..1] range. */
263        real bn[ORDER + 1];
264
265        for (int i = 0; i < ORDER + 1; i++)
266        {
267            bn[i] = 0;
268            for (int j = 0; j < ORDER + 1; j++)
269                bn[i] += coeff[j] * (real)Cheby(j, i);
270        }
271
272        /* Transform a polynomial in the [-1..1] range into a polynomial
273         * in the [a..b] range. */
274        real k1p[ORDER + 1], k2p[ORDER + 1];
275        real an[ORDER + 1];
276
277        for (int i = 0; i < ORDER + 1; i++)
278        {
279            k1p[i] = i ? k1p[i - 1] * m_invk1 : real::R_1;
280            k2p[i] = i ? k2p[i - 1] * m_invk2 : real::R_1;
281        }
282
283        for (int i = 0; i < ORDER + 1; i++)
284        {
285            an[i] = 0;
286            for (int j = i; j < ORDER + 1; j++)
287                an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j];
288            an[i] *= k2p[i];
289        }
290
291        printf("Polynomial estimate:\n");
292        for (int j = 0; j < ORDER + 1; j++)
293        {
294            if (j)
295                printf("+");
296            printf("x^%i*", j);
297            an[j].print(40);
298        }
299        printf("\n");
300    }
301
302    real Value(real const &x)
303    {
304        return m_func(x * m_k2 + m_k1);
305    }
306
307    real Weight(real const &x)
308    {
309        return m_weight(x * m_k2 + m_k1);
310    }
311
312    /* ORDER + 1 Chebyshev coefficients and 1 error value */
313    real coeff[ORDER + 2];
314    /* ORDER + 1 zeroes of the error function */
315    real zeroes[ORDER + 1];
316    /* ORDER + 2 control points */
317    real control[ORDER + 2];
318
319private:
320    RealFunc *m_func, *m_weight;
321    real m_k1, m_k2, m_invk1, m_invk2;
322};
323
324#endif /* __REMEZ_SOLVER_H__ */
325
Note: See TracBrowser for help on using the repository browser.