# Changeset 989 for trunk/test/mathTweet

Ignore:
Timestamp:
Sep 27, 2011, 7:09:35 PM (10 years ago)
Message:

test: more work on the Remez exchange algorithm.

Location:
trunk/test/math
Files:
2 edited

Unmodified
Added
Removed
• ## trunk/test/math/pi.cpp

 r984 sum.print(); /* Bonus: compute e for fun. */ sum = 0.0; x0 = 1.0; for (int i = 1; i < 100; i++) { sum += fres(x0); x0 *= (real)i; } sum.print(); return EXIT_SUCCESS; }
• ## trunk/test/math/remez.cpp

 r985 using namespace lol; static int const ORDER = 12; /* The order of the approximation we're looking for */ static int const ORDER = 3; static int cheby[ORDER][ORDER]; /* The function we want to approximate */ double myfun(double x) { return exp(x); } real myfun(real const &x) { return exp(x); } /* Naive matrix inversion */ template class Matrix { public: inline Matrix() {} Matrix(real x) { for (int j = 0; j < N; j++) for (int i = 0; i < N; i++) if (i == j) m[i][j] = x; else m[i][j] = 0; } Matrix inv() const { Matrix a = *this, b(real(1.0)); /* Inversion method: iterate through all columns and make sure * all the terms are zero except on the diagonal where it is one */ for (int i = 0; i < N; i++) { /* If the expected coefficient is zero, add one of * the other lines. The first we meet will do. */ if ((double)a.m[i][i] == 0.0) { for (int j = i + 1; j < N; j++) { if ((double)a.m[i][j] == 0.0) continue; /* Add row j to row i */ for (int n = 0; n < N; n++) { a.m[n][i] += a.m[n][j]; b.m[n][i] += b.m[n][j]; } break; } } /* Now we know the diagonal term is non-zero. Get its inverse * and use that to nullify all other terms in the column */ real x = fres(a.m[i][i]); for (int j = 0; j < N; j++) { if (j == i) continue; real mul = x * a.m[i][j]; for (int n = 0; n < N; n++) { a.m[n][j] -= mul * a.m[n][i]; b.m[n][j] -= mul * b.m[n][i]; } } /* Finally, ensure the diagonal term is 1 */ for (int n = 0; n < N; n++) { a.m[n][i] *= x; b.m[n][i] *= x; } } return b; } void print() const { for (int j = 0; j < N; j++) { for (int i = 0; i < N; i++) printf("%9.5f ", (double)m[j][i]); printf("\n"); } } real m[N][N]; }; static int cheby[ORDER + 1][ORDER + 1]; /* Fill the Chebyshev tables */ static void make_table() { memset(cheby[0], 0, ORDER * sizeof(int)); memset(cheby, 0, sizeof(cheby)); cheby[0][0] = 1; memset(cheby[1], 0, ORDER * sizeof(int)); cheby[1][1] = 1; for (int i = 2; i < ORDER; i++) for (int i = 2; i < ORDER + 1; i++) { cheby[i][0] = -cheby[i - 2][0]; for (int j = 1; j < ORDER; j++) for (int j = 1; j < ORDER + 1; j++) cheby[i][j] = 2 * cheby[i - 1][j - 1] - cheby[i - 2][j]; } make_table(); for (int i = 0; i < ORDER; i++) { for (int j = 0; j < ORDER; j++) printf("% 5i ", cheby[i][j]); printf("\n"); } /* We start with ORDER+1 points and their images through myfun() */ real xn[ORDER + 1]; real fxn[ORDER + 1]; for (int i = 0; i < ORDER + 1; i++) { //xn[i] = real(2 * i - ORDER) / real(ORDER + 1); xn[i] = real(2 * i - ORDER + 1) / real(ORDER - 1); fxn[i] = myfun(xn[i]); } /* We build a matrix of Chebishev evaluations: one row per point * in our array, and column i is the evaluation of the ith order * polynomial. */ Matrix mat; for (int j = 0; j < ORDER + 1; j++) { /* Compute the powers of x_j */ real powers[ORDER + 1]; powers[0] = 1.0; for (int i = 1; i < ORDER + 1; i++) powers[i] = powers[i - 1] * xn[j]; /* Compute the Chebishev evaluations at x_j */ for (int i = 0; i < ORDER + 1; i++) { real sum = 0.0; for (int k = 0; k < ORDER + 1; k++) if (cheby[i][k]) sum += real(cheby[i][k]) * powers[k]; mat.m[j][i] = sum; } } /* Invert the matrix and build interpolation coefficients */ mat = mat.inv(); real an[ORDER + 1]; for (int j = 0; j < ORDER + 1; j++) { an[j] = 0; for (int i = 0; i < ORDER + 1; i++) an[j] += mat.m[j][i] * fxn[i]; an[j].print(10); } return EXIT_SUCCESS;
Note: See TracChangeset for help on using the changeset viewer.