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

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

test: use namespace "std" to avoid PS3 build errors.

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