Changeset 1010
 Timestamp:
 Oct 5, 2011, 7:01:33 PM (9 years ago)
 Location:
 trunk/test/math
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test/math/remezsolver.h
r1009 r1010 19 19 RemezSolver() 20 20 { 21 ChebyInit(); 22 } 23 24 void Run(RealFunc *func, RealFunc *error, int steps) 21 } 22 23 void Run(real a, real b, RealFunc *func, RealFunc *error, int steps) 25 24 { 26 25 m_func = func; 27 26 m_error = error; 27 m_k1 = (b + a) >> 1; 28 m_k2 = (b  a) >> 1; 29 m_invk2 = re(m_k2); 30 m_invk1 = m_k1 * m_invk2; 28 31 29 32 Init(); 30 33 31 ChebyCoeff(); 32 for (int j = 0; j < ORDER + 1; j++) 33 printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j); 34 printf("\n"); 34 PrintPoly(); 35 35 36 36 for (int n = 0; n < steps; n++) … … 39 39 Step(); 40 40 41 ChebyCoeff(); 42 for (int j = 0; j < ORDER + 1; j++) 43 printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j); 44 printf("\n"); 41 PrintPoly(); 45 42 46 43 FindZeroes(); … … 50 47 Step(); 51 48 52 ChebyCoeff(); 53 for (int j = 0; j < ORDER + 1; j++) 54 printf("%s%14.12gx^%i", j && (bn[j] >= real::R_0) ? "+" : "", (double)bn[j], j); 55 printf("\n"); 56 } 57 58 /* Fill the Chebyshev tables */ 59 void ChebyInit() 60 { 61 memset(cheby, 0, sizeof(cheby)); 62 63 cheby[0][0] = 1; 64 cheby[1][1] = 1; 65 66 for (int i = 2; i < ORDER + 1; i++) 67 { 68 cheby[i][0] = cheby[i  2][0]; 69 for (int j = 1; j < ORDER + 1; j++) 70 cheby[i][j] = 2 * cheby[i  1][j  1]  cheby[i  2][j]; 71 } 72 } 73 74 void ChebyCoeff() 75 { 76 for (int i = 0; i < ORDER + 1; i++) 77 { 78 bn[i] = 0; 79 for (int j = 0; j < ORDER + 1; j++) 80 if (cheby[j][i]) 81 bn[i] += coeff[j] * (real)cheby[j][i]; 82 } 49 PrintPoly(); 83 50 } 84 51 … … 91 58 real mul = 0; 92 59 for (int j = 0; j < ORDER + 1; j++) 93 if (cheby[j][i]) 94 mul += coeff[j] * (real)cheby[j][i]; 60 mul += coeff[j] * (real)Cheby(j, i); 95 61 ret += mul * xn; 96 62 xn *= x; … … 107 73 { 108 74 zeroes[i] = (real)(2 * i  ORDER) / (real)(ORDER + 1); 109 fxn[i] = m_func(zeroes[i]);75 fxn[i] = Value(zeroes[i]); 110 76 } 111 77 … … 126 92 real sum = 0.0; 127 93 for (int k = 0; k < ORDER + 1; k++) 128 if (cheby[n][k]) 129 sum += (real)cheby[n][k] * powers[k]; 94 sum += (real)Cheby(n, k) * powers[k]; 130 95 mat.m[i][n] = sum; 131 96 } … … 149 114 { 150 115 real a = control[i]; 151 real ea = ChebyEval(a)  m_func(a);116 real ea = ChebyEval(a)  Value(a); 152 117 real b = control[i + 1]; 153 real eb = ChebyEval(b)  m_func(b);118 real eb = ChebyEval(b)  Value(b); 154 119 155 120 while (fabs(a  b) > (real)1e140) 156 121 { 157 122 real c = (a + b) * (real)0.5; 158 real ec = ChebyEval(c)  m_func(c);123 real ec = ChebyEval(c)  Value(c); 159 124 160 125 if ((ea < (real)0 && ec < (real)0) … … 195 160 for (int k = 0; k <= 10; k++) 196 161 { 197 real e = fabs(ChebyEval(c)  m_func(c));162 real e = fabs(ChebyEval(c)  Value(c)); 198 163 if (e > maxerror) 199 164 { … … 223 188 } 224 189 225 printf("Final error: %g\n", (double)final); 190 printf("Final error: "); 191 final.print(40); 226 192 } 227 193 … … 231 197 real fxn[ORDER + 2]; 232 198 for (int i = 0; i < ORDER + 2; i++) 233 fxn[i] = m_func(control[i]);199 fxn[i] = Value(control[i]); 234 200 235 201 /* We build a matrix of Chebishev evaluations: row i contains the … … 249 215 real sum = 0.0; 250 216 for (int k = 0; k < ORDER + 1; k++) 251 if (cheby[n][k]) 252 sum += (real)cheby[n][k] * powers[k]; 217 sum += (real)Cheby(n, k) * powers[k]; 253 218 mat.m[i][n] = sum; 254 219 } 255 220 if (i & 1) 256 mat.m[i][ORDER + 1] = fabs( m_error(control[i]));221 mat.m[i][ORDER + 1] = fabs(Error(control[i])); 257 222 else 258 mat.m[i][ORDER + 1] = fabs( m_error(control[i]));223 mat.m[i][ORDER + 1] = fabs(Error(control[i])); 259 224 } 260 225 … … 276 241 } 277 242 278 int cheby[ORDER + 1][ORDER + 1]; 279 280 /* ORDER + 1 chebyshev coefficients and 1 error value */ 243 int Cheby(int n, int k) 244 { 245 if (k > n  k < 0) 246 return 0; 247 if (n <= 1) 248 return (n ^ k ^ 1) & 1; 249 return 2 * Cheby(n  1, k  1)  Cheby(n  2, k); 250 } 251 252 int Comb(int n, int k) 253 { 254 if (k == 0  k == n) 255 return 1; 256 return Comb(n  1, k  1) + Comb(n  1, k); 257 } 258 259 void PrintPoly() 260 { 261 /* Transform Chebyshev polynomial weights into powers of X^i 262 * in the [1..1] range. */ 263 real bn[ORDER + 1]; 264 265 for (int i = 0; i < ORDER + 1; i++) 266 { 267 bn[i] = 0; 268 for (int j = 0; j < ORDER + 1; j++) 269 bn[i] += coeff[j] * (real)Cheby(j, i); 270 } 271 272 /* Transform a polynomial in the [1..1] range into a polynomial 273 * in the [a..b] range. */ 274 real k1p[ORDER + 1], k2p[ORDER + 1]; 275 real an[ORDER + 1]; 276 277 for (int i = 0; i < ORDER + 1; i++) 278 { 279 k1p[i] = i ? k1p[i  1] * m_invk1 : real::R_1; 280 k2p[i] = i ? k2p[i  1] * m_invk2 : real::R_1; 281 } 282 283 for (int i = 0; i < ORDER + 1; i++) 284 { 285 an[i] = 0; 286 for (int j = i; j < ORDER + 1; j++) 287 an[i] += (real)Comb(j, i) * k1p[j  i] * bn[j]; 288 an[i] *= k2p[i]; 289 } 290 291 for (int j = 0; j < ORDER + 1; j++) 292 printf("%s%18.16gx^%i", j && (an[j] >= real::R_0) ? "+" : "", (double)an[j], j); 293 printf("\n"); 294 } 295 296 real Value(real const &x) 297 { 298 return m_func(x * m_k2 + m_k1); 299 } 300 301 real Error(real const &x) 302 { 303 return m_error(x * m_k2 + m_k1); 304 } 305 306 /* ORDER + 1 Chebyshev coefficients and 1 error value */ 281 307 real coeff[ORDER + 2]; 282 308 /* ORDER + 1 zeroes of the error function */ … … 285 311 real control[ORDER + 2]; 286 312 287 real bn[ORDER + 1];288 289 313 private: 290 RealFunc *m_func ;291 RealFunc *m_error;314 RealFunc *m_func, *m_error; 315 real m_k1, m_k2, m_invk1, m_invk2; 292 316 }; 293 317 
trunk/test/math/remez.cpp
r1009 r1010 27 27 static real myfun(real const &x) 28 28 { 29 return exp(x); 29 real y = sqrt(x); 30 if (!y) 31 return real::R_PI_2; 32 return sin(real::R_PI_2 * y) / y; 30 33 } 31 34 32 static real myerr or(real const &x)35 static real myerr(real const &x) 33 36 { 34 return myfun(x);37 return real::R_1; 35 38 } 36 39 … … 38 41 { 39 42 RemezSolver<4> solver; 40 41 solver.Run(myfun, myerror, 10); 43 solver.Run(0, 1, myfun, myfun, 15); 44 //solver.Run(1, 1, myfun, myfun, 15); 45 //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15); 42 46 43 47 return EXIT_SUCCESS;
Note: See TracChangeset
for help on using the changeset viewer.