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

Last change on this file since 1164 was 1164, checked in by sam, 8 years ago

ps3: fix PS3 build after the Queue refactoring.

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