Changeset 1127
 Timestamp:
 Jan 9, 2012, 11:49:28 PM (7 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/lol/math/remez.h
r1122 r1127 33 33 } 34 34 35 void Run(T a, T b, RealFunc *func, RealFunc *weight, int steps)35 void Run(T a, T b, RealFunc *func, RealFunc *weight, int decimals) 36 36 { 37 37 m_func = func; … … 41 41 m_invk2 = re(m_k2); 42 42 m_invk1 = m_k1 * m_invk2; 43 m_decimals = decimals; 44 m_epsilon = pow((T)10, (T)(decimals + 2)); 43 45 44 46 Init(); … … 46 48 PrintPoly(); 47 49 48 for (int n = 0; n < steps; n++) 49 { 50 FindExtrema(); 50 T error = 1; 51 52 for (int n = 0; ; n++) 53 { 54 T newerror = FindExtrema(); 55 printf("Step %i error: ", n); 56 newerror.print(m_decimals); 57 51 58 Step(); 52 59 60 if (error >= (T)0 && fabs(newerror  error) < error * m_epsilon) 61 break; 62 error = newerror; 63 53 64 PrintPoly(); 54 65 … … 56 67 } 57 68 58 FindExtrema();59 Step();60 61 69 PrintPoly(); 62 70 } 63 71 64 inline void Run(T a, T b, RealFunc *func, int steps)65 { 66 Run(a, b, func, NULL, steps);72 inline void Run(T a, T b, RealFunc *func, int decimals) 73 { 74 Run(a, b, func, NULL, decimals); 67 75 } 68 76 … … 158 166 } 159 167 160 voidFindExtrema()168 real FindExtrema() 161 169 { 162 170 using std::printf; … … 175 183 b = zeroes[i]; 176 184 177 for (;;) 185 T maxerror = 0, maxweight = 0; 186 int best = 1; 187 188 for (int round = 0; ; round++) 178 189 { 179 T c = a, delta = (b  a) / 8; 180 T maxerror = 0; 181 T maxweight = 0; 182 int best = 1; 183 for (int k = 1; k <= 7; k++) 190 T c = a, delta = (b  a) / 4; 191 for (int k = 0; k <= 4; k++) 184 192 { 185 T error = ChebyEval(c)  Value(c); 186 T weight = Weight(c); 187 if (fabs(error * maxweight) >= fabs(maxerror * weight)) 193 if (round == 0  (k & 1)) 188 194 { 189 maxerror = error; 190 maxweight = weight; 191 best = k; 195 T error = ChebyEval(c)  Value(c); 196 T weight = Weight(c); 197 if (fabs(error * maxweight) >= fabs(maxerror * weight)) 198 { 199 maxerror = error; 200 maxweight = weight; 201 best = k; 202 } 192 203 } 193 204 c += delta; 194 205 } 195 206 196 b = a + (T)(best + 1) * delta; 197 a = a + (T)(best  1) * delta; 198 199 if (b  a < (T)1e18) 207 switch (best) 208 { 209 case 0: 210 b = a + delta * 2; 211 break; 212 case 4: 213 a = b  delta * 2; 214 break; 215 default: 216 b = a + delta * (best + 1); 217 a = a + delta * (best  1); 218 best = 2; 219 break; 220 } 221 222 if (delta < m_epsilon) 200 223 { 201 224 T e = maxerror / maxweight; … … 208 231 } 209 232 210 printf("Final error: "); 211 final.print(40); 233 return final; 212 234 } 213 235 … … 317 339 printf("+"); 318 340 printf("x**%i*", j); 319 an[j].print( 40);341 an[j].print(m_decimals); 320 342 } 321 343 printf("\n"); … … 343 365 private: 344 366 RealFunc *m_func, *m_weight; 345 T m_k1, m_k2, m_invk1, m_invk2; 367 T m_k1, m_k2, m_invk1, m_invk2, m_epsilon; 368 int m_decimals; 346 369 }; 347 370 
trunk/test/math/remez.cpp
r1126 r1127 30 30 { 31 31 RemezSolver<4, real> solver; 32 solver.Run(1, 1, f, g, 30);32 solver.Run(1, 1, f, g, 40); 33 33 return 0; 34 34 }
Note: See TracChangeset
for help on using the changeset viewer.