Changeset 1117
- Timestamp:
- Dec 29, 2011, 7:36:48 PM (11 years ago)
- Location:
- trunk
- Files:
-
- 1 added
- 1 deleted
- 10 edited
- 3 moved
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/Makefile.am
r1106 r1117 3 3 4 4 liblol_a_SOURCES = \ 5 core.h matrix.cpp matrix.htiler.cpp tiler.h dict.cpp dict.h \5 core.h matrix.cpp real.cpp tiler.cpp tiler.h dict.cpp dict.h \ 6 6 audio.cpp audio.h scene.cpp scene.h font.cpp font.h layer.cpp layer.h \ 7 7 map.cpp map.h entity.cpp entity.h ticker.cpp ticker.h lolgl.h \ … … 12 12 worldentity.cpp worldentity.h gradient.cpp gradient.h half.cpp half.h \ 13 13 platform.cpp platform.h sprite.cpp sprite.h trig.cpp trig.h \ 14 real.cpp real.h \15 14 \ 16 15 lol/unit.h \ 16 lol/math/matrix.h lol/math/real.h lol/math/remez.h \ 17 17 \ 18 18 application/application.cpp application/application.h \ -
trunk/src/core.h
r1106 r1117 65 65 #include "trig.h" 66 66 #include "half.h" 67 #include " real.h"68 #include " matrix.h"67 #include "lol/math/real.h" 68 #include "lol/math/matrix.h" 69 69 #include "numeric.h" 70 70 #include "timer.h" -
trunk/src/eglapp.h
r863 r1117 17 17 #define __LOL_EGLAPP_H__ 18 18 19 #include " matrix.h"19 #include "lol/math/matrix.h" 20 20 21 21 namespace lol -
trunk/src/image/image.h
r864 r1117 17 17 #define __LOL_IMAGE_H__ 18 18 19 #include " matrix.h"19 #include "lol/math/matrix.h" 20 20 21 21 namespace lol -
trunk/src/input.h
r866 r1117 17 17 #define __LOL_INPUT_H__ 18 18 19 #include " matrix.h"19 #include "lol/math/matrix.h" 20 20 21 21 namespace lol -
trunk/src/lol/math/matrix.h
r1116 r1117 14 14 // 15 15 16 #if !defined __LOL_MAT RIX_H__17 #define __LOL_MAT RIX_H__16 #if !defined __LOL_MATH_MATRIX_H__ 17 #define __LOL_MATH_MATRIX_H__ 18 18 19 19 #include <cmath> … … 24 24 namespace lol 25 25 { 26 27 class half; 28 class real; 26 29 27 30 #define VECTOR_TYPES(tname, suffix) \ … … 588 591 }; 589 592 593 /* 594 * Arbitrarily-sized square matrices; for now this only supports 595 * naive inversion and is used for the Remez inversion method. 596 */ 597 598 template<int N, typename T> struct Mat 599 { 600 inline Mat<N, T>() {} 601 602 Mat(T x) 603 { 604 for (int j = 0; j < N; j++) 605 for (int i = 0; i < N; i++) 606 if (i == j) 607 m[i][j] = x; 608 else 609 m[i][j] = 0; 610 } 611 612 /* Naive matrix inversion */ 613 Mat<N, T> inv() const 614 { 615 Mat a = *this, b((T)1); 616 617 /* Inversion method: iterate through all columns and make sure 618 * all the terms are 1 on the diagonal and 0 everywhere else */ 619 for (int i = 0; i < N; i++) 620 { 621 /* If the expected coefficient is zero, add one of 622 * the other lines. The first we meet will do. */ 623 if (!a.m[i][i]) 624 { 625 for (int j = i + 1; j < N; j++) 626 { 627 if (!a.m[i][j]) 628 continue; 629 /* Add row j to row i */ 630 for (int n = 0; n < N; n++) 631 { 632 a.m[n][i] += a.m[n][j]; 633 b.m[n][i] += b.m[n][j]; 634 } 635 break; 636 } 637 } 638 639 /* Now we know the diagonal term is non-zero. Get its inverse 640 * and use that to nullify all other terms in the column */ 641 T x = (T)1 / a.m[i][i]; 642 for (int j = 0; j < N; j++) 643 { 644 if (j == i) 645 continue; 646 T mul = x * a.m[i][j]; 647 for (int n = 0; n < N; n++) 648 { 649 a.m[n][j] -= mul * a.m[n][i]; 650 b.m[n][j] -= mul * b.m[n][i]; 651 } 652 } 653 654 /* Finally, ensure the diagonal term is 1 */ 655 for (int n = 0; n < N; n++) 656 { 657 a.m[n][i] *= x; 658 b.m[n][i] *= x; 659 } 660 } 661 662 return b; 663 } 664 665 T m[N][N]; 666 }; 667 590 668 } /* namespace lol */ 591 669 592 #endif // __LOL_MAT RIX_H__593 670 #endif // __LOL_MATH_MATRIX_H__ 671 -
trunk/src/lol/math/real.h
-
Property
svn:keywords
set to
Id
r1116 r1117 14 14 // 15 15 16 #if !defined __LOL_ REAL_H__17 #define __LOL_ REAL_H__16 #if !defined __LOL_MATH_REAL_H__ 17 #define __LOL_MATH_REAL_H__ 18 18 19 19 #include <stdint.h> … … 148 148 } /* namespace lol */ 149 149 150 #endif // __LOL_ REAL_H__150 #endif // __LOL_MATH_REAL_H__ 151 151 -
Property
svn:keywords
set to
-
trunk/src/lol/math/remez.h
r1116 r1117 1 1 // 2 // Lol Engine - Sample math program: Chebyshev polynomials3 // 4 // Copyright: (c) 20 05-2011 Sam Hocevar <sam@hocevar.net>2 // Lol Engine 3 // 4 // Copyright: (c) 2010-2011 Sam Hocevar <sam@hocevar.net> 5 5 // This program is free software; you can redistribute it and/or 6 6 // modify it under the terms of the Do What The Fuck You Want To … … 9 9 // 10 10 11 #if !defined __REMEZ_SOLVER_H__ 12 #define __REMEZ_SOLVER_H__ 13 14 template<int ORDER> class RemezSolver 11 // 12 // The Remez class 13 // --------------- 14 // 15 16 #if !defined __LOL_MATH_REMEZ_H__ 17 #define __LOL_MATH_REMEZ_H__ 18 19 namespace lol 20 { 21 22 template<int ORDER, typename T> class RemezSolver 15 23 { 16 24 public: 17 typedef real RealFunc(realconst &x);25 typedef T RealFunc(T const &x); 18 26 19 27 RemezSolver() … … 21 29 } 22 30 23 void Run( real a, realb, RealFunc *func, RealFunc *weight, int steps)31 void Run(T a, T b, RealFunc *func, RealFunc *weight, int steps) 24 32 { 25 33 m_func = func; … … 50 58 } 51 59 52 real ChebyEval(realconst &x)53 { 54 realret = 0.0, xn = 1.0;55 56 for (int i = 0; i < ORDER + 1; i++) 57 { 58 realmul = 0;60 T ChebyEval(T const &x) 61 { 62 T ret = 0.0, xn = 1.0; 63 64 for (int i = 0; i < ORDER + 1; i++) 65 { 66 T mul = 0; 59 67 for (int j = 0; j < ORDER + 1; j++) 60 mul += coeff[j] * ( real)Cheby(j, i);68 mul += coeff[j] * (T)Cheby(j, i); 61 69 ret += mul * xn; 62 70 xn *= x; … … 69 77 { 70 78 /* Pick up x_i where error will be 0 and compute f(x_i) */ 71 realfxn[ORDER + 1];72 for (int i = 0; i < ORDER + 1; i++) 73 { 74 zeroes[i] = ( real)(2 * i - ORDER) / (real)(ORDER + 1);79 T fxn[ORDER + 1]; 80 for (int i = 0; i < ORDER + 1; i++) 81 { 82 zeroes[i] = (T)(2 * i - ORDER) / (T)(ORDER + 1); 75 83 fxn[i] = Value(zeroes[i]); 76 84 } … … 78 86 /* We build a matrix of Chebishev evaluations: row i contains the 79 87 * evaluations of x_i for polynomial order n = 0, 1, ... */ 80 Matrix<ORDER + 1> mat;88 lol::Mat<ORDER + 1, T> mat; 81 89 for (int i = 0; i < ORDER + 1; i++) 82 90 { 83 91 /* Compute the powers of x_i */ 84 realpowers[ORDER + 1];92 T powers[ORDER + 1]; 85 93 powers[0] = 1.0; 86 94 for (int n = 1; n < ORDER + 1; n++) … … 90 98 for (int n = 0; n < ORDER + 1; n++) 91 99 { 92 realsum = 0.0;100 T sum = 0.0; 93 101 for (int k = 0; k < ORDER + 1; k++) 94 sum += ( real)Cheby(n, k) * powers[k];102 sum += (T)Cheby(n, k) * powers[k]; 95 103 mat.m[i][n] = sum; 96 104 } … … 116 124 for (int i = 0; i < ORDER + 1; i++) 117 125 { 118 struct { realvalue, error; } left, right, mid;126 struct { T value, error; } left, right, mid; 119 127 120 128 left.value = control[i]; … … 123 131 right.error = ChebyEval(right.value) - Value(right.value); 124 132 125 static real limit = real::R_1 >> 500; 133 static T limit = (T)1 >> 500; 134 static T zero = (T)0; 126 135 while (fabs(left.value - right.value) > limit) 127 136 { … … 129 138 mid.error = ChebyEval(mid.value) - Value(mid.value); 130 139 131 if ((left.error < real::R_0 && mid.error < real::R_0)132 || (left.error > real::R_0 && mid.error > real::R_0))140 if ((left.error < zero && mid.error < zero) 141 || (left.error > zero && mid.error > zero)) 133 142 left = mid; 134 143 else … … 147 156 * compute the relative error, since its extrema are at slightly 148 157 * different locations than the absolute error’s. */ 149 realfinal = 0;158 T final = 0; 150 159 151 160 for (int i = 0; i < ORDER + 2; i++) 152 161 { 153 reala = -1, b = 1;162 T a = -1, b = 1; 154 163 if (i > 0) 155 164 a = zeroes[i - 1]; … … 159 168 for (;;) 160 169 { 161 realc = a, delta = (b - a) >> 3;162 realmaxerror = 0;163 realmaxweight = 0;170 T c = a, delta = (b - a) >> 3; 171 T maxerror = 0; 172 T maxweight = 0; 164 173 int best = -1; 165 174 for (int k = 1; k <= 7; k++) 166 175 { 167 realerror = ChebyEval(c) - Value(c);168 realweight = Weight(c);176 T error = ChebyEval(c) - Value(c); 177 T weight = Weight(c); 169 178 if (fabs(error * maxweight) >= fabs(maxerror * weight)) 170 179 { … … 176 185 } 177 186 178 b = a + ( real)(best + 1) * delta;179 a = a + ( real)(best - 1) * delta;180 181 if (b - a < ( real)1e-18)187 b = a + (T)(best + 1) * delta; 188 a = a + (T)(best - 1) * delta; 189 190 if (b - a < (T)1e-18) 182 191 { 183 reale = maxerror / maxweight;192 T e = maxerror / maxweight; 184 193 if (e > final) 185 194 final = e; … … 197 206 { 198 207 /* Pick up x_i where error will be 0 and compute f(x_i) */ 199 realfxn[ORDER + 2];208 T fxn[ORDER + 2]; 200 209 for (int i = 0; i < ORDER + 2; i++) 201 210 fxn[i] = Value(control[i]); … … 203 212 /* We build a matrix of Chebishev evaluations: row i contains the 204 213 * evaluations of x_i for polynomial order n = 0, 1, ... */ 205 Matrix<ORDER + 2> mat;214 lol::Mat<ORDER + 2, T> mat; 206 215 for (int i = 0; i < ORDER + 2; i++) 207 216 { 208 217 /* Compute the powers of x_i */ 209 realpowers[ORDER + 1];218 T powers[ORDER + 1]; 210 219 powers[0] = 1.0; 211 220 for (int n = 1; n < ORDER + 1; n++) … … 215 224 for (int n = 0; n < ORDER + 1; n++) 216 225 { 217 realsum = 0.0;226 T sum = 0.0; 218 227 for (int k = 0; k < ORDER + 1; k++) 219 sum += ( real)Cheby(n, k) * powers[k];228 sum += (T)Cheby(n, k) * powers[k]; 220 229 mat.m[i][n] = sum; 221 230 } … … 238 247 239 248 /* Compute the error */ 240 realerror = 0;249 T error = 0; 241 250 for (int i = 0; i < ORDER + 2; i++) 242 251 error += mat.m[ORDER + 1][i] * fxn[i]; … … 265 274 /* Transform Chebyshev polynomial weights into powers of X^i 266 275 * in the [-1..1] range. */ 267 realbn[ORDER + 1];276 T bn[ORDER + 1]; 268 277 269 278 for (int i = 0; i < ORDER + 1; i++) … … 271 280 bn[i] = 0; 272 281 for (int j = 0; j < ORDER + 1; j++) 273 bn[i] += coeff[j] * ( real)Cheby(j, i);282 bn[i] += coeff[j] * (T)Cheby(j, i); 274 283 } 275 284 276 285 /* Transform a polynomial in the [-1..1] range into a polynomial 277 286 * in the [a..b] range. */ 278 realk1p[ORDER + 1], k2p[ORDER + 1];279 realan[ORDER + 1];280 281 for (int i = 0; i < ORDER + 1; i++) 282 { 283 k1p[i] = i ? k1p[i - 1] * m_invk1 : real::R_1;284 k2p[i] = i ? k2p[i - 1] * m_invk2 : real::R_1;287 T k1p[ORDER + 1], k2p[ORDER + 1]; 288 T an[ORDER + 1]; 289 290 for (int i = 0; i < ORDER + 1; i++) 291 { 292 k1p[i] = i ? k1p[i - 1] * m_invk1 : (T)1; 293 k2p[i] = i ? k2p[i - 1] * m_invk2 : (T)1; 285 294 } 286 295 … … 289 298 an[i] = 0; 290 299 for (int j = i; j < ORDER + 1; j++) 291 an[i] += ( real)Comb(j, i) * k1p[j - i] * bn[j];300 an[i] += (T)Comb(j, i) * k1p[j - i] * bn[j]; 292 301 an[i] *= k2p[i]; 293 302 } … … 298 307 if (j) 299 308 printf("+"); 300 printf("x ^%i*", j);309 printf("x**%i*", j); 301 310 an[j].print(40); 302 311 } … … 304 313 } 305 314 306 real Value(realconst &x)315 T Value(T const &x) 307 316 { 308 317 return m_func(x * m_k2 + m_k1); 309 318 } 310 319 311 real Weight(realconst &x)320 T Weight(T const &x) 312 321 { 313 322 return m_weight(x * m_k2 + m_k1); … … 315 324 316 325 /* ORDER + 1 Chebyshev coefficients and 1 error value */ 317 realcoeff[ORDER + 2];326 T coeff[ORDER + 2]; 318 327 /* ORDER + 1 zeroes of the error function */ 319 realzeroes[ORDER + 1];328 T zeroes[ORDER + 1]; 320 329 /* ORDER + 2 control points */ 321 realcontrol[ORDER + 2];330 T control[ORDER + 2]; 322 331 323 332 private: 324 333 RealFunc *m_func, *m_weight; 325 realm_k1, m_k2, m_invk1, m_invk2;334 T m_k1, m_k2, m_invk1, m_invk2; 326 335 }; 327 336 328 #endif /* __REMEZ_SOLVER_H__ */ 329 337 } /* namespace lol */ 338 339 #endif /* __LOL_MATH_REMEZ_H__ */ 340 -
trunk/src/platform/nacl/naclapp.h
r1082 r1117 17 17 #define __LOL_NACLAPP_H__ 18 18 19 #include " matrix.h"19 #include "lol/math/matrix.h" 20 20 21 21 namespace lol -
trunk/src/platform/ps3/ps3app.h
r1063 r1117 17 17 #define __LOL_PS3APP_H__ 18 18 19 #include " matrix.h"19 #include "lol/math/matrix.h" 20 20 21 21 namespace lol -
trunk/src/platform/sdl/sdlapp.h
r1063 r1117 17 17 #define __LOL_SDLAPP_H__ 18 18 19 #include " matrix.h"19 #include "lol/math/matrix.h" 20 20 21 21 namespace lol -
trunk/test/math/Makefile.am
r1112 r1117 23 23 poly_DEPENDENCIES = $(top_builddir)/src/liblol.a 24 24 25 remez_SOURCES = remez.cpp remez-matrix.h remez-solver.h25 remez_SOURCES = remez.cpp 26 26 remez_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@ 27 27 remez_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@ -
trunk/test/math/remez.cpp
r1091 r1117 13 13 #endif 14 14 15 #include <cstring>16 15 #include <cstdio> 16 #include <cstdlib> 17 17 18 18 #if USE_SDL && defined __APPLE__ … … 20 20 #endif 21 21 22 #include "core.h" 22 #include "lol/math/real.h" 23 #include "lol/math/matrix.h" 24 #include "lol/math/remez.h" 23 25 24 26 using lol::real; 25 26 #include "remez-matrix.h" 27 #include "remez-solver.h" 27 using lol::RemezSolver; 28 28 29 29 /* The function we want to approximate */ … … 42 42 int main(int argc, char **argv) 43 43 { 44 RemezSolver<2 > solver;44 RemezSolver<2, real> solver; 45 45 solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 40); 46 46
Note: See TracChangeset
for help on using the changeset viewer.