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

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

test: add support for relative error in the Remez exchange test program.

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