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 | |
---|
19 | using namespace lol; |
---|
20 | |
---|
21 | /* The order of the approximation we're looking for */ |
---|
22 | static int const ORDER = 3; |
---|
23 | |
---|
24 | /* The function we want to approximate */ |
---|
25 | double myfun(double x) |
---|
26 | { |
---|
27 | return exp(x); |
---|
28 | } |
---|
29 | |
---|
30 | real myfun(real const &x) |
---|
31 | { |
---|
32 | return exp(x); |
---|
33 | } |
---|
34 | |
---|
35 | /* Naive matrix inversion */ |
---|
36 | template<int N> class Matrix |
---|
37 | { |
---|
38 | public: |
---|
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 | |
---|
117 | static int cheby[ORDER + 1][ORDER + 1]; |
---|
118 | |
---|
119 | /* Fill the Chebyshev tables */ |
---|
120 | static 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 | |
---|
135 | int 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 | |
---|