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

Last change on this file since 2506 was 2506, checked in by sam, 7 years ago

base: start removing occurrences of NULL on our long journey to nullptr.

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