Changeset 989 for trunk/test/math/remez.cpp
 Timestamp:
 Sep 27, 2011, 7:09:35 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test/math/remez.cpp
r985 r989 19 19 using namespace lol; 20 20 21 static int const ORDER = 12; 21 /* The order of the approximation we're looking for */ 22 static int const ORDER = 3; 22 23 23 static int cheby[ORDER][ORDER]; 24 /* The function we want to approximate */ 25 double myfun(double x) 26 { 27 return exp(x); 28 } 24 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 nonzero. 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 */ 25 120 static void make_table() 26 121 { 27 memset(cheby[0], 0, ORDER * sizeof(int)); 122 memset(cheby, 0, sizeof(cheby)); 123 28 124 cheby[0][0] = 1; 29 memset(cheby[1], 0, ORDER * sizeof(int));30 125 cheby[1][1] = 1; 31 126 32 for (int i = 2; i < ORDER ; i++)127 for (int i = 2; i < ORDER + 1; i++) 33 128 { 34 129 cheby[i][0] = cheby[i  2][0]; 35 for (int j = 1; j < ORDER ; j++)130 for (int j = 1; j < ORDER + 1; j++) 36 131 cheby[i][j] = 2 * cheby[i  1][j  1]  cheby[i  2][j]; 37 132 } … … 42 137 make_table(); 43 138 44 for (int i = 0; i < ORDER; i++) 45 { 46 for (int j = 0; j < ORDER; j++) 47 printf("% 5i ", cheby[i][j]); 48 printf("\n"); 49 } 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 } 50 182 51 183 return EXIT_SUCCESS;
Note: See TracChangeset
for help on using the changeset viewer.