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

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

test: more work on the Remez exchange algorithm.

  • Property svn:keywords set to Id
File size: 4.6 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 = 3;
23
24/* The function we want to approximate */
25double myfun(double x)
26{
27    return exp(x);
28}
29
30real myfun(real const &x)
31{
32    return exp(x);
33}
34
35/* Naive matrix inversion */
36template<int N> class Matrix
37{
38public:
39    inline Matrix() {}
40
41    Matrix(real x)
42    {
43        for (int j = 0; j < N; j++)
44            for (int i = 0; i < N; i++)
45                if (i == j)
46                    m[i][j] = x;
47                else
48                    m[i][j] = 0;
49    }
50
51    Matrix<N> inv() const
52    {
53        Matrix a = *this, b(real(1.0));
54
55        /* Inversion method: iterate through all columns and make sure
56         * all the terms are zero except on the diagonal where it is one */
57        for (int i = 0; i < N; i++)
58        {
59            /* If the expected coefficient is zero, add one of
60             * the other lines. The first we meet will do. */
61            if ((double)a.m[i][i] == 0.0)
62            {
63                for (int j = i + 1; j < N; j++)
64                {
65                    if ((double)a.m[i][j] == 0.0)
66                        continue;
67                    /* Add row j to row i */
68                    for (int n = 0; n < N; n++)
69                    {
70                        a.m[n][i] += a.m[n][j];
71                        b.m[n][i] += b.m[n][j];
72                    }
73                    break;
74                }
75            }
76
77            /* Now we know the diagonal term is non-zero. Get its inverse
78             * and use that to nullify all other terms in the column */
79            real x = fres(a.m[i][i]);
80            for (int j = 0; j < N; j++)
81            {
82                if (j == i)
83                    continue;
84                real mul = x * a.m[i][j];
85                for (int n = 0; n < N; n++)
86                {
87                    a.m[n][j] -= mul * a.m[n][i];
88                    b.m[n][j] -= mul * b.m[n][i];
89                }
90            }
91
92            /* Finally, ensure the diagonal term is 1 */
93            for (int n = 0; n < N; n++)
94            {
95                a.m[n][i] *= x;
96                b.m[n][i] *= x;
97            }
98        }
99
100        return b;
101    }
102
103    void print() const
104    {
105        for (int j = 0; j < N; j++)
106        {
107            for (int i = 0; i < N; i++)
108                printf("%9.5f ", (double)m[j][i]);
109            printf("\n");
110        }
111    }
112
113    real m[N][N];
114};
115
116
117static int cheby[ORDER + 1][ORDER + 1];
118
119/* Fill the Chebyshev tables */
120static void make_table()
121{
122    memset(cheby, 0, sizeof(cheby));
123
124    cheby[0][0] = 1;
125    cheby[1][1] = 1;
126
127    for (int i = 2; i < ORDER + 1; i++)
128    {
129        cheby[i][0] = -cheby[i - 2][0];
130        for (int j = 1; j < ORDER + 1; j++)
131            cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j];
132    }
133}
134
135int main(int argc, char **argv)
136{
137    make_table();
138
139    /* We start with ORDER+1 points and their images through myfun() */
140    real xn[ORDER + 1];
141    real fxn[ORDER + 1];
142    for (int i = 0; i < ORDER + 1; i++)
143    {
144        //xn[i] = real(2 * i - ORDER) / real(ORDER + 1);
145        xn[i] = real(2 * i - ORDER + 1) / real(ORDER - 1);
146        fxn[i] = myfun(xn[i]);
147    }
148
149    /* We build a matrix of Chebishev evaluations: one row per point
150     * in our array, and column i is the evaluation of the ith order
151     * polynomial. */
152    Matrix<ORDER + 1> mat;
153    for (int j = 0; j < ORDER + 1; j++)
154    {
155        /* Compute the powers of x_j */
156        real powers[ORDER + 1];
157        powers[0] = 1.0;
158        for (int i = 1; i < ORDER + 1; i++)
159             powers[i] = powers[i - 1] * xn[j];
160
161        /* Compute the Chebishev evaluations at x_j */
162        for (int i = 0; i < ORDER + 1; i++)
163        {
164            real sum = 0.0;
165            for (int k = 0; k < ORDER + 1; k++)
166                if (cheby[i][k])
167                    sum += real(cheby[i][k]) * powers[k];
168            mat.m[j][i] = sum;
169        }
170    }
171
172    /* Invert the matrix and build interpolation coefficients */
173    mat = mat.inv();
174    real an[ORDER + 1];
175    for (int j = 0; j < ORDER + 1; j++)
176    {
177        an[j] = 0;
178        for (int i = 0; i < ORDER + 1; i++)
179            an[j] += mat.m[j][i] * fxn[i];
180        an[j].print(10);
181    }
182
183    return EXIT_SUCCESS;
184}
185
Note: See TracBrowser for help on using the repository browser.