source: trunk/src/lol/math/remez.h @ 1139

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

math: rename matrix.h to vector.h and simplify some stuff, especially in
the matrix code itself.

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