Changeset 991
 Timestamp:
 Sep 28, 2011, 2:59:42 AM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test/math/remez.cpp
r989 r991 20 20 21 21 /* The order of the approximation we're looking for */ 22 static int const ORDER = 3;22 static int const ORDER = 4; 23 23 24 24 /* The function we want to approximate */ 25 double myfun(double x)26 {27 return exp(x);28 }29 30 25 real myfun(real const &x) 31 26 { 32 return exp(x); 27 static real const one = 1.0; 28 return cos(x)  one; 29 //return exp(x); 33 30 } 34 31 35 32 /* Naive matrix inversion */ 36 template<int N> class Matrix 37 { 38 public: 33 template<int N> struct Matrix 34 { 39 35 inline Matrix() {} 40 36 … … 54 50 55 51 /* Inversion method: iterate through all columns and make sure 56 * all the terms are zero except on the diagonal where it is one */52 * all the terms are 1 on the diagonal and 0 everywhere else */ 57 53 for (int i = 0; i < N; i++) 58 54 { … … 118 114 119 115 /* Fill the Chebyshev tables */ 120 static void make_table()116 static void cheby_init() 121 117 { 122 118 memset(cheby, 0, sizeof(cheby)); … … 133 129 } 134 130 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]; 131 static void cheby_coeff(real *coeff, real *bn) 132 { 133 for (int i = 0; i < ORDER + 1; i++) 134 { 135 bn[i] = 0; 136 for (int j = 0; j < ORDER + 1; j++) 137 if (cheby[j][i]) 138 bn[i] += coeff[j] * (real)cheby[j][i]; 139 } 140 } 141 142 static real cheby_eval(real *coeff, real const &x) 143 { 144 real ret = 0.0, xn = 1.0; 145 146 for (int i = 0; i < ORDER + 1; i++) 147 { 148 real mul = 0; 149 for (int j = 0; j < ORDER + 1; j++) 150 if (cheby[j][i]) 151 mul += coeff[j] * (real)cheby[j][i]; 152 ret += mul * xn; 153 xn *= x; 154 } 155 156 return ret; 157 } 158 159 static void remez_init(real *coeff, real *zeroes) 160 { 161 /* Pick up x_i where error will be 0 and compute f(x_i) */ 141 162 real fxn[ORDER + 1]; 142 163 for (int i = 0; i < ORDER + 1; i++) 143 164 { 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. */ 165 zeroes[i] = real(2 * i  ORDER) / real(ORDER + 1); 166 fxn[i] = myfun(zeroes[i]); 167 } 168 169 /* We build a matrix of Chebishev evaluations: row i contains the 170 * evaluations of x_i for polynomial order n = 0, 1, ... */ 152 171 Matrix<ORDER + 1> mat; 153 for (int j = 0; j < ORDER + 1; j++)154 { 155 /* Compute the powers of x_ j*/172 for (int i = 0; i < ORDER + 1; i++) 173 { 174 /* Compute the powers of x_i */ 156 175 real powers[ORDER + 1]; 157 176 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++)177 for (int n = 1; n < ORDER + 1; n++) 178 powers[n] = powers[n  1] * zeroes[i]; 179 180 /* Compute the Chebishev evaluations at x_i */ 181 for (int n = 0; n < ORDER + 1; n++) 163 182 { 164 183 real sum = 0.0; 165 184 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*/185 if (cheby[n][k]) 186 sum += real(cheby[n][k]) * powers[k]; 187 mat.m[i][n] = sum; 188 } 189 } 190 191 /* Solve the system */ 173 192 mat = mat.inv(); 174 real an[ORDER + 1]; 193 194 /* Compute interpolation coefficients */ 175 195 for (int j = 0; j < ORDER + 1; j++) 176 196 { 177 an[j] = 0;197 coeff[j] = 0; 178 198 for (int i = 0; i < ORDER + 1; i++) 179 an[j] += mat.m[j][i] * fxn[i]; 180 an[j].print(10); 181 } 199 coeff[j] += mat.m[j][i] * fxn[i]; 200 } 201 } 202 203 static void remez_findzeroes(real *coeff, real *zeroes, real *control) 204 { 205 /* FIXME: this is fake for now */ 206 for (int i = 0; i < ORDER + 1; i++) 207 { 208 real a = control[i]; 209 real ea = cheby_eval(coeff, a)  myfun(a); 210 real b = control[i + 1]; 211 real eb = cheby_eval(coeff, b)  myfun(b); 212 213 while (fabs(a  b) > (real)1e140) 214 { 215 real c = (a + b) * (real)0.5; 216 real ec = cheby_eval(coeff, c)  myfun(c); 217 218 if ((ea < (real)0 && ec < (real)0) 219  (ea > (real)0 && ec > (real)0)) 220 { 221 a = c; 222 ea = ec; 223 } 224 else 225 { 226 b = c; 227 eb = ec; 228 } 229 } 230 231 zeroes[i] = a; 232 } 233 } 234 235 static void remez_finderror(real *coeff, real *zeroes, real *control) 236 { 237 real final = 0; 238 239 for (int i = 0; i < ORDER + 2; i++) 240 { 241 real a = 1, b = 1; 242 if (i > 0) 243 a = zeroes[i  1]; 244 if (i < ORDER + 1) 245 b = zeroes[i]; 246 247 printf("Error for [%g..%g]: ", (double)a, (double)b); 248 for (;;) 249 { 250 real c = a, delta = (b  a) / (real)10.0; 251 real maxerror = 0; 252 int best = 1; 253 for (int k = 0; k <= 10; k++) 254 { 255 real e = fabs(cheby_eval(coeff, c)  myfun(c)); 256 if (e > maxerror) 257 { 258 maxerror = e; 259 best = k; 260 } 261 c += delta; 262 } 263 264 if (best == 0) 265 best = 1; 266 if (best == 10) 267 best = 9; 268 269 b = a + (real)(best + 1) * delta; 270 a = a + (real)(best  1) * delta; 271 272 if (b  a < (real)1e15) 273 { 274 if (maxerror > final) 275 final = maxerror; 276 control[i] = (a + b) * (real)0.5; 277 printf("%g (in %g)\n", (double)maxerror, (double)control[i]); 278 break; 279 } 280 } 281 } 282 283 printf("Final error: %g\n", (double)final); 284 } 285 286 static void remez_step(real *coeff, real *control) 287 { 288 /* Pick up x_i where error will be 0 and compute f(x_i) */ 289 real fxn[ORDER + 2]; 290 for (int i = 0; i < ORDER + 2; i++) 291 fxn[i] = myfun(control[i]); 292 293 /* We build a matrix of Chebishev evaluations: row i contains the 294 * evaluations of x_i for polynomial order n = 0, 1, ... */ 295 Matrix<ORDER + 2> mat; 296 for (int i = 0; i < ORDER + 2; i++) 297 { 298 /* Compute the powers of x_i */ 299 real powers[ORDER + 1]; 300 powers[0] = 1.0; 301 for (int n = 1; n < ORDER + 1; n++) 302 powers[n] = powers[n  1] * control[i]; 303 304 /* Compute the Chebishev evaluations at x_i */ 305 for (int n = 0; n < ORDER + 1; n++) 306 { 307 real sum = 0.0; 308 for (int k = 0; k < ORDER + 1; k++) 309 if (cheby[n][k]) 310 sum += real(cheby[n][k]) * powers[k]; 311 mat.m[i][n] = sum; 312 } 313 mat.m[i][ORDER + 1] = (real)(1 + (i & 1) * 2); 314 } 315 316 /* Solve the system */ 317 mat = mat.inv(); 318 319 /* Compute interpolation coefficients */ 320 for (int j = 0; j < ORDER + 1; j++) 321 { 322 coeff[j] = 0; 323 for (int i = 0; i < ORDER + 2; i++) 324 coeff[j] += mat.m[j][i] * fxn[i]; 325 } 326 327 /* Compute the error */ 328 real error = 0; 329 for (int i = 0; i < ORDER + 2; i++) 330 error += mat.m[ORDER + 1][i] * fxn[i]; 331 } 332 333 int main(int argc, char **argv) 334 { 335 cheby_init(); 336 337 /* ORDER + 1 chebyshev coefficients and 1 error value */ 338 real coeff[ORDER + 2]; 339 /* ORDER + 1 zeroes of the error function */ 340 real zeroes[ORDER + 1]; 341 /* ORDER + 2 control points */ 342 real control[ORDER + 2]; 343 344 real bn[ORDER + 1]; 345 346 remez_init(coeff, zeroes); 347 348 cheby_coeff(coeff, bn); 349 for (int j = 0; j < ORDER + 1; j++) 350 printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j); 351 printf("\n"); 352 353 for (int n = 0; n < 200; n++) 354 { 355 remez_finderror(coeff, zeroes, control); 356 remez_step(coeff, control); 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 remez_findzeroes(coeff, zeroes, control); 364 } 365 366 remez_finderror(coeff, zeroes, control); 367 remez_step(coeff, control); 368 369 cheby_coeff(coeff, bn); 370 for (int j = 0; j < ORDER + 1; j++) 371 printf("%s%12.10gx^%i", j ? "+" : "", (double)bn[j], j); 372 printf("\n"); 182 373 183 374 return EXIT_SUCCESS;
Note: See TracChangeset
for help on using the changeset viewer.