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

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

test: minor updates to the Pi program (now almost deprecated) and the
Remez exchange program.

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