Changeset 1014
 Timestamp:
 Oct 10, 2011, 3:33:15 AM (9 years ago)
 Location:
 trunk/test/math
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test/math/remezsolver.h
r1010 r1014 21 21 } 22 22 23 void Run(real a, real b, RealFunc *func, RealFunc * error, int steps)23 void Run(real a, real b, RealFunc *func, RealFunc *weight, int steps) 24 24 { 25 25 m_func = func; 26 m_ error = error;26 m_weight = weight; 27 27 m_k1 = (b + a) >> 1; 28 28 m_k2 = (b  a) >> 1; … … 36 36 for (int n = 0; n < steps; n++) 37 37 { 38 FindE rror();38 FindExtrema(); 39 39 Step(); 40 40 … … 44 44 } 45 45 46 FindE rror();46 FindExtrema(); 47 47 Step(); 48 48 … … 111 111 void FindZeroes() 112 112 { 113 for (int i = 0; i < ORDER + 1; i++) 114 { 115 real a = control[i]; 116 real ea = ChebyEval(a)  Value(a); 117 real b = control[i + 1]; 118 real eb = ChebyEval(b)  Value(b); 119 120 while (fabs(a  b) > (real)1e140) 113 /* Find ORDER + 1 zeroes of the error function. No need to 114 * compute the relative error: its zeroes are at the same 115 * place as the absolute error! */ 116 for (int i = 0; i < ORDER + 1; i++) 117 { 118 struct { real value, error; } left, right, mid; 119 120 left.value = control[i]; 121 left.error = ChebyEval(left.value)  Value(left.value); 122 right.value = control[i + 1]; 123 right.error = ChebyEval(right.value)  Value(right.value); 124 125 static real limit = real::R_1 >> 500; 126 while (fabs(left.value  right.value) > limit) 121 127 { 122 real c = (a + b) * (real)0.5; 123 real ec = ChebyEval(c)  Value(c); 124 125 if ((ea < (real)0 && ec < (real)0) 126  (ea > (real)0 && ec > (real)0)) 127 { 128 a = c; 129 ea = ec; 130 } 128 mid.value = (left.value + right.value) >> 1; 129 mid.error = ChebyEval(mid.value)  Value(mid.value); 130 131 if ((left.error < real::R_0 && mid.error < real::R_0) 132  (left.error > real::R_0 && mid.error > real::R_0)) 133 left = mid; 131 134 else 132 { 133 b = c; 134 eb = ec; 135 } 135 right = mid; 136 136 } 137 137 138 zeroes[i] = a; 139 } 140 } 141 142 void FindError() 143 { 138 zeroes[i] = mid.value; 139 } 140 } 141 142 void FindExtrema() 143 { 144 /* Find ORDER + 2 extrema of the error function. We need to 145 * compute the relative error, since its extrema are at slightly 146 * different locations than the absolute error’s. */ 144 147 real final = 0; 145 148 … … 155 158 for (;;) 156 159 { 157 real c = a, delta = (b  a) / (real)10.0;160 real c = a, delta = (b  a) >> 3; 158 161 real maxerror = 0; 162 real maxweight = 0; 159 163 int best = 1; 160 for (int k = 0; k <= 10; k++)164 for (int k = 1; k <= 7; k++) 161 165 { 162 real e = fabs(ChebyEval(c)  Value(c)); 163 if (e > maxerror) 166 real error = ChebyEval(c)  Value(c); 167 real weight = Weight(c); 168 if (fabs(error * maxweight) >= fabs(maxerror * weight)) 164 169 { 165 maxerror = e; 170 maxerror = error; 171 maxweight = weight; 166 172 best = k; 167 173 } … … 169 175 } 170 176 171 if (best == 0)172 best = 1;173 if (best == 10)174 best = 9;175 176 177 b = a + (real)(best + 1) * delta; 177 178 a = a + (real)(best  1) * delta; 178 179 179 if (b  a < (real)1e1 5)180 if (b  a < (real)1e18) 180 181 { 181 if (maxerror > final) 182 final = maxerror; 183 control[i] = (a + b) * (real)0.5; 184 printf("%g (at %g)\n", (double)maxerror, (double)control[i]); 182 real e = maxerror / maxweight; 183 if (e > final) 184 final = e; 185 control[i] = (a + b) >> 1; 186 printf("%g (at %g)\n", (double)e, (double)control[i]); 185 187 break; 186 188 } … … 219 221 } 220 222 if (i & 1) 221 mat.m[i][ORDER + 1] = fabs( Error(control[i]));223 mat.m[i][ORDER + 1] = fabs(Weight(control[i])); 222 224 else 223 mat.m[i][ORDER + 1] = fabs( Error(control[i]));225 mat.m[i][ORDER + 1] = fabs(Weight(control[i])); 224 226 } 225 227 … … 289 291 } 290 292 293 printf("Final polynomial:\n"); 291 294 for (int j = 0; j < ORDER + 1; j++) 292 printf("%s%18.16gx^%i", j && (an[j] >= real::R_0) ? "+" : "", (double)an[j], j); 295 { 296 if (j) 297 printf("+"); 298 printf("x^%i*", j); 299 an[j].print(40); 300 } 293 301 printf("\n"); 294 302 } … … 299 307 } 300 308 301 real Error(real const &x)302 { 303 return m_ error(x * m_k2 + m_k1);309 real Weight(real const &x) 310 { 311 return m_weight(x * m_k2 + m_k1); 304 312 } 305 313 … … 312 320 313 321 private: 314 RealFunc *m_func, *m_ error;322 RealFunc *m_func, *m_weight; 315 323 real m_k1, m_k2, m_invk1, m_invk2; 316 324 }; 
trunk/test/math/remez.cpp
r1012 r1014 27 27 static real myfun(real const &x) 28 28 { 29 static real const a0 = real::R_1;30 static real const a1 = real(11184811) >> 26;31 static real const b1 = real(1) / real(6);32 static real const b2 = real(1) / real(120);33 29 real y = sqrt(x); 34 return sin(y) / y;30 return (sin(y)  y) / (x * y); 35 31 } 36 32 37 33 static real myerr(real const &x) 38 34 { 39 return real::R_1; 35 real y = sqrt(x); 36 return re(x * y); 40 37 } 41 38 42 int main( int argc, char **argv)39 int main(void) 43 40 { 44 RemezSolver<6> solver; 45 //solver.Run(0, 1, myfun, myfun, 15); 46 solver.Run(exp((real)100), real::R_PI * real::R_PI >> 2, myfun, myfun, 15); 47 //solver.Run(1, 1, myfun, myfun, 15); 48 //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15); 41 RemezSolver<4> solver; 42 solver.Run(real::R_1 >> 400, real::R_PI_2 * real::R_PI_2, myfun, myerr, 15); 49 43 50 44 return EXIT_SUCCESS;
Note: See TracChangeset
for help on using the changeset viewer.