source: trunk/test/math/remez.cpp @ 997

Last change on this file since 997 was 997, checked in by sam, 11 years ago

core: add boolean operators on real numbers, add unit tests for that,
and simplify the Remez code accordingly.

  • Property svn:keywords set to Id
File size: 9.7 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 HAVE_CONFIG_H
12#   include "config.h"
13#endif
14
15#include <cstring>
16
17#include "core.h"
18
19using namespace lol;
20
21/* The order of the approximation we're looking for */
22static int const ORDER = 8;
23
24/* The function we want to approximate */
25static real myfun(real const &x)
26{
27    static real const one = 1.0;
28    if (!x)
29        return real::R_PI_2;
30    return sin(x * real::R_PI_2) / x;
31    //return cos(x) - one;
32    //return exp(x);
33}
34
35/* Naive matrix inversion */
36template<int N> struct Matrix
37{
38    inline Matrix() {}
39
40    Matrix(real x)
41    {
42        for (int j = 0; j < N; j++)
43            for (int i = 0; i < N; i++)
44                if (i == j)
45                    m[i][j] = x;
46                else
47                    m[i][j] = 0;
48    }
49
50    Matrix<N> inv() const
51    {
52        Matrix a = *this, b((real)1.0);
53
54        /* Inversion method: iterate through all columns and make sure
55         * all the terms are 1 on the diagonal and 0 everywhere else */
56        for (int i = 0; i < N; i++)
57        {
58            /* If the expected coefficient is zero, add one of
59             * the other lines. The first we meet will do. */
60            if ((double)a.m[i][i] == 0.0)
61            {
62                for (int j = i + 1; j < N; j++)
63                {
64                    if ((double)a.m[i][j] == 0.0)
65                        continue;
66                    /* Add row j to row i */
67                    for (int n = 0; n < N; n++)
68                    {
69                        a.m[n][i] += a.m[n][j];
70                        b.m[n][i] += b.m[n][j];
71                    }
72                    break;
73                }
74            }
75
76            /* Now we know the diagonal term is non-zero. Get its inverse
77             * and use that to nullify all other terms in the column */
78            real x = (real)1.0 / a.m[i][i];
79            for (int j = 0; j < N; j++)
80            {
81                if (j == i)
82                    continue;
83                real mul = x * a.m[i][j];
84                for (int n = 0; n < N; n++)
85                {
86                    a.m[n][j] -= mul * a.m[n][i];
87                    b.m[n][j] -= mul * b.m[n][i];
88                }
89            }
90
91            /* Finally, ensure the diagonal term is 1 */
92            for (int n = 0; n < N; n++)
93            {
94                a.m[n][i] *= x;
95                b.m[n][i] *= x;
96            }
97        }
98
99        return b;
100    }
101
102    void print() const
103    {
104        for (int j = 0; j < N; j++)
105        {
106            for (int i = 0; i < N; i++)
107                printf("%9.5f ", (double)m[j][i]);
108            printf("\n");
109        }
110    }
111
112    real m[N][N];
113};
114
115
116static int cheby[ORDER + 1][ORDER + 1];
117
118/* Fill the Chebyshev tables */
119static void cheby_init()
120{
121    memset(cheby, 0, sizeof(cheby));
122
123    cheby[0][0] = 1;
124    cheby[1][1] = 1;
125
126    for (int i = 2; i < ORDER + 1; i++)
127    {
128        cheby[i][0] = -cheby[i - 2][0];
129        for (int j = 1; j < ORDER + 1; j++)
130            cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j];
131    }
132}
133
134static void cheby_coeff(real *coeff, real *bn)
135{
136    for (int i = 0; i < ORDER + 1; i++)
137    {
138        bn[i] = 0;
139        for (int j = 0; j < ORDER + 1; j++)
140            if (cheby[j][i])
141                bn[i] += coeff[j] * (real)cheby[j][i];
142    }
143}
144
145static real cheby_eval(real *coeff, real const &x)
146{
147    real ret = 0.0, xn = 1.0;
148
149    for (int i = 0; i < ORDER + 1; i++)
150    {
151        real mul = 0;
152        for (int j = 0; j < ORDER + 1; j++)
153            if (cheby[j][i])
154                mul += coeff[j] * (real)cheby[j][i];
155        ret += mul * xn;
156        xn *= x;
157    }
158
159    return ret;
160}
161
162static void remez_init(real *coeff, real *zeroes)
163{
164    /* Pick up x_i where error will be 0 and compute f(x_i) */
165    real fxn[ORDER + 1];
166    for (int i = 0; i < ORDER + 1; i++)
167    {
168        zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
169        fxn[i] = myfun(zeroes[i]);
170    }
171
172    /* We build a matrix of Chebishev evaluations: row i contains the
173     * evaluations of x_i for polynomial order n = 0, 1, ... */
174    Matrix<ORDER + 1> mat;
175    for (int i = 0; i < ORDER + 1; i++)
176    {
177        /* Compute the powers of x_i */
178        real powers[ORDER + 1];
179        powers[0] = 1.0;
180        for (int n = 1; n < ORDER + 1; n++)
181             powers[n] = powers[n - 1] * zeroes[i];
182
183        /* Compute the Chebishev evaluations at x_i */
184        for (int n = 0; n < ORDER + 1; n++)
185        {
186            real sum = 0.0;
187            for (int k = 0; k < ORDER + 1; k++)
188                if (cheby[n][k])
189                    sum += (real)cheby[n][k] * powers[k];
190            mat.m[i][n] = sum;
191        }
192    }
193
194    /* Solve the system */
195    mat = mat.inv();
196
197    /* Compute interpolation coefficients */
198    for (int j = 0; j < ORDER + 1; j++)
199    {
200        coeff[j] = 0;
201        for (int i = 0; i < ORDER + 1; i++)
202            coeff[j] += mat.m[j][i] * fxn[i];
203    }
204}
205
206static void remez_findzeroes(real *coeff, real *zeroes, real *control)
207{
208    /* FIXME: this is fake for now */
209    for (int i = 0; i < ORDER + 1; i++)
210    {
211        real a = control[i];
212        real ea = cheby_eval(coeff, a) - myfun(a);
213        real b = control[i + 1];
214        real eb = cheby_eval(coeff, b) - myfun(b);
215
216        while (fabs(a - b) > (real)1e-140)
217        {
218            real c = (a + b) * (real)0.5;
219            real ec = cheby_eval(coeff, c) - myfun(c);
220
221            if ((ea < (real)0 && ec < (real)0)
222                 || (ea > (real)0 && ec > (real)0))
223            {
224                a = c;
225                ea = ec;
226            }
227            else
228            {
229                b = c;
230                eb = ec;
231            }
232        }
233
234        zeroes[i] = a;
235    }
236}
237
238static void remez_finderror(real *coeff, real *zeroes, real *control)
239{
240    real final = 0;
241
242    for (int i = 0; i < ORDER + 2; i++)
243    {
244        real a = -1, b = 1;
245        if (i > 0)
246            a = zeroes[i - 1];
247        if (i < ORDER + 1)
248            b = zeroes[i];
249
250        printf("Error for [%g..%g]: ", (double)a, (double)b);
251        for (;;)
252        {
253            real c = a, delta = (b - a) / (real)10.0;
254            real maxerror = 0;
255            int best = -1;
256            for (int k = 0; k <= 10; k++)
257            {
258                real e = fabs(cheby_eval(coeff, c) - myfun(c));
259                if (e > maxerror)
260                {
261                    maxerror = e;
262                    best = k;
263                }
264                c += delta;
265            }
266
267            if (best == 0)
268                best = 1;
269            if (best == 10)
270                best = 9;
271
272            b = a + (real)(best + 1) * delta;
273            a = a + (real)(best - 1) * delta;
274
275            if (b - a < (real)1e-15)
276            {
277                if (maxerror > final)
278                    final = maxerror;
279                control[i] = (a + b) * (real)0.5;
280                printf("%g (in %g)\n", (double)maxerror, (double)control[i]);
281                break;
282            }
283        }
284    }
285
286    printf("Final error: %g\n", (double)final);
287}
288
289static void remez_step(real *coeff, real *control)
290{
291    /* Pick up x_i where error will be 0 and compute f(x_i) */
292    real fxn[ORDER + 2];
293    for (int i = 0; i < ORDER + 2; i++)
294        fxn[i] = myfun(control[i]);
295
296    /* We build a matrix of Chebishev evaluations: row i contains the
297     * evaluations of x_i for polynomial order n = 0, 1, ... */
298    Matrix<ORDER + 2> mat;
299    for (int i = 0; i < ORDER + 2; i++)
300    {
301        /* Compute the powers of x_i */
302        real powers[ORDER + 1];
303        powers[0] = 1.0;
304        for (int n = 1; n < ORDER + 1; n++)
305             powers[n] = powers[n - 1] * control[i];
306
307        /* Compute the Chebishev evaluations at x_i */
308        for (int n = 0; n < ORDER + 1; n++)
309        {
310            real sum = 0.0;
311            for (int k = 0; k < ORDER + 1; k++)
312                if (cheby[n][k])
313                    sum += (real)cheby[n][k] * powers[k];
314            mat.m[i][n] = sum;
315        }
316        mat.m[i][ORDER + 1] = (real)(-1 + (i & 1) * 2);
317    }
318
319    /* Solve the system */
320    mat = mat.inv();
321
322    /* Compute interpolation coefficients */
323    for (int j = 0; j < ORDER + 1; j++)
324    {
325        coeff[j] = 0;
326        for (int i = 0; i < ORDER + 2; i++)
327            coeff[j] += mat.m[j][i] * fxn[i];
328    }
329
330    /* Compute the error */
331    real error = 0;
332    for (int i = 0; i < ORDER + 2; i++)
333        error += mat.m[ORDER + 1][i] * fxn[i];
334}
335
336int main(int argc, char **argv)
337{
338    cheby_init();
339
340    /* ORDER + 1 chebyshev coefficients and 1 error value */
341    real coeff[ORDER + 2];
342    /* ORDER + 1 zeroes of the error function */
343    real zeroes[ORDER + 1];
344    /* ORDER + 2 control points */
345    real control[ORDER + 2];
346
347    real bn[ORDER + 1];
348
349    remez_init(coeff, zeroes);
350
351    cheby_coeff(coeff, bn);
352    for (int j = 0; j < ORDER + 1; j++)
353        printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
354    printf("\n");
355
356    for (int n = 0; n < 200; n++)
357    {
358        remez_finderror(coeff, zeroes, control);
359        remez_step(coeff, control);
360
361        cheby_coeff(coeff, bn);
362        for (int j = 0; j < ORDER + 1; j++)
363            printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
364        printf("\n");
365
366        remez_findzeroes(coeff, zeroes, control);
367    }
368
369    remez_finderror(coeff, zeroes, control);
370    remez_step(coeff, control);
371
372    cheby_coeff(coeff, bn);
373    for (int j = 0; j < ORDER + 1; j++)
374        printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j);
375    printf("\n");
376
377    return EXIT_SUCCESS;
378}
379
Note: See TracBrowser for help on using the repository browser.