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

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

test: allow to perform Remez solving on an arbitrary range.

  • Property svn:keywords set to Id
File size: 8.4 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 *error, int steps)
24    {
25        m_func = func;
26        m_error = error;
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            FindError();
39            Step();
40
41            PrintPoly();
42
43            FindZeroes();
44        }
45
46        FindError();
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        for (int i = 0; i < ORDER + 1; i++)
114        {
115            real a = control[i];
116            real ea = ChebyEval(a) - Value(a);
117            real b = control[i + 1];
118            real eb = ChebyEval(b) - Value(b);
119
120            while (fabs(a - b) > (real)1e-140)
121            {
122                real c = (a + b) * (real)0.5;
123                real ec = ChebyEval(c) - Value(c);
124
125                if ((ea < (real)0 && ec < (real)0)
126                     || (ea > (real)0 && ec > (real)0))
127                {
128                    a = c;
129                    ea = ec;
130                }
131                else
132                {
133                    b = c;
134                    eb = ec;
135                }
136            }
137
138            zeroes[i] = a;
139        }
140    }
141
142    void FindError()
143    {
144        real final = 0;
145
146        for (int i = 0; i < ORDER + 2; i++)
147        {
148            real a = -1, b = 1;
149            if (i > 0)
150                a = zeroes[i - 1];
151            if (i < ORDER + 1)
152                b = zeroes[i];
153
154            printf("Error for [%g..%g]: ", (double)a, (double)b);
155            for (;;)
156            {
157                real c = a, delta = (b - a) / (real)10.0;
158                real maxerror = 0;
159                int best = -1;
160                for (int k = 0; k <= 10; k++)
161                {
162                    real e = fabs(ChebyEval(c) - Value(c));
163                    if (e > maxerror)
164                    {
165                        maxerror = e;
166                        best = k;
167                    }
168                    c += delta;
169                }
170
171                if (best == 0)
172                    best = 1;
173                if (best == 10)
174                    best = 9;
175
176                b = a + (real)(best + 1) * delta;
177                a = a + (real)(best - 1) * delta;
178
179                if (b - a < (real)1e-15)
180                {
181                    if (maxerror > final)
182                        final = maxerror;
183                    control[i] = (a + b) * (real)0.5;
184                    printf("%g (at %g)\n", (double)maxerror, (double)control[i]);
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(Error(control[i]));
222            else
223                mat.m[i][ORDER + 1] = -fabs(Error(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        for (int j = 0; j < ORDER + 1; j++)
292            printf("%s%18.16gx^%i", j && (an[j] >= real::R_0) ? "+" : "", (double)an[j], j);
293        printf("\n");
294    }
295
296    real Value(real const &x)
297    {
298        return m_func(x * m_k2 + m_k1);
299    }
300
301    real Error(real const &x)
302    {
303        return m_error(x * m_k2 + m_k1);
304    }
305
306    /* ORDER + 1 Chebyshev coefficients and 1 error value */
307    real coeff[ORDER + 2];
308    /* ORDER + 1 zeroes of the error function */
309    real zeroes[ORDER + 1];
310    /* ORDER + 2 control points */
311    real control[ORDER + 2];
312
313private:
314    RealFunc *m_func, *m_error;
315    real m_k1, m_k2, m_invk1, m_invk2;
316};
317
318#endif /* __REMEZ_SOLVER_H__ */
319
Note: See TracBrowser for help on using the repository browser.