Changeset 1012
- Timestamp:
- Oct 6, 2011, 6:55:47 PM (11 years ago)
- Location:
- trunk/test/math
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/test/math/poly.cpp
r1011 r1012 20 20 using namespace std; 21 21 22 static float adjust (float f, int i)22 static float adjustf(float f, int i) 23 23 { 24 24 union { float f; uint32_t x; } u = { f }; 25 u.x += i; 26 return u.f; 27 } 28 29 static double adjust(double f, int i) 30 { 31 union { double f; uint64_t x; } u = { f }; 25 32 u.x += i; 26 33 return u.f; … … 38 45 } 39 46 47 //#define float double 48 #if 1 49 static float const a0 = adjust(0.9999999999999376, 0); 50 static float const a1 = adjust(-0.1666666666643236, 0); 51 static float const a2 = adjust(0.008333333318766562, 0); 52 static float const a3 = adjust(-0.0001984126641174625, 0); 53 static float const a4 = adjust(2.755693193297308e-006, 0); 54 static float const a5 = adjust(-2.502951900290311e-008, 0); 55 static float const a6 = adjust(1.540117123154927e-010, 0); 56 #elif 0 40 57 static float const a0 = adjust(1.0, 0); 41 58 static float const a1 = adjust(-0.1666666666372165, 0); … … 44 61 static float const a4 = adjust(2.753619853326498e-06, 0); 45 62 static float const a5 = adjust(-2.407483949485896e-08, 0); 63 static float const a6 = 0.0; 64 #else 65 static float const a0 = adjust(0.9999999946887117, 0); 66 static float const a1 = adjust(-0.1666665668590824, 0); 67 static float const a2 = adjust(0.008333025160523476, 0); 68 static float const a3 = adjust(-0.0001980741944205014, 0); 69 static float const a4 = adjust(2.60190356966559e-06, 0); // -900 in floats 70 static float const a5 = 0.0; 71 static float const a6 = 0.0; 72 #endif 46 73 47 74 static float floatsin(float f) 48 75 { 49 76 //return lol_sin(f); 50 //#define float double51 77 //static float const k = (float)real::R_2_PI; 52 78 … … 54 80 float f2 = f * f; 55 81 float f4 = f2 * f2; 56 //return f * (a0 + f4 * (a2 + f4 * a4) + f2 * (a1 + f4 * (a3 + f4 * a5)));57 //return f * (a0 + f2 * (a1 + f2 * (a2 + f2 * (a3 + f2 * (a4 + f2 * a5)))));58 return f * (a0 + a1 * f2 + a2 * f2 * f2 + a3 * f2 * f2 * f2 + a4 * f2 * f2 * f2 * f2 + a5* f2 * f2 * f2 * f2 * f2);82 return f * (a0 + f4 * (a2 + f4 * (a4 + f4 * a6)) + f2 * (a1 + f4 * (a3 + f4 * a5))); 83 //return f * (a0 + f2 * (a1 + f2 * (a2 + f2 * (a3 + f2 * (a4 + f2 * (a5 + f2 * a6)))))); 84 //return f * (a0 + a1 * f2 + a2 * f2 * f2 + a3 * f2 * f2 * f2 + a4 * f2 * f2 * f2 * f2 + a5 * f2 * f2 * f2 * f2 * f2 + a6 * f2 * f2 * f2 * f2 * f2 * f2); 59 85 #undef float 60 86 } … … 71 97 inspect(a5); 72 98 73 for (float f = (float)real::R_PI_2; f > 1e-30; f = adjust (f, -1))99 for (float f = (float)real::R_PI_2; f > 1e-30; f = adjustf(f, -1)) 74 100 { 101 union { float f; uint32_t x; } u = { f }; 75 102 union { float f; uint32_t x; } s1 = { sinf(f) }; 76 103 union { float f; uint32_t x; } s2 = { floatsin(f) }; … … 78 105 switch (e) 79 106 { 107 case 3: 108 if (fabs((double)s1.f - (double)s2.f) > 1e-10 * fabs(s1.f)) 109 printf("%15.13g %08x: %15.13g %15.13g %08x %08x\n", f, u.x, s1.f, s2.f, s1.x, s2.x); 110 case 2: 111 case 1: 80 112 case 0: 81 case 1:82 case 2:83 case 3:84 113 error[e]++; 85 114 break; 86 115 default: 87 if (fabs((double)s1.f - (double)s2.f) > 1e-10 * fabs(s1.f))88 printf("%16.14g: %16.14g %16.14g %08x %08x\n", f, s1.f, s2.f, s1.x, s2.x);89 116 error[4]++; 90 117 break; -
trunk/test/math/remez.cpp
r1011 r1012 32 32 static real const b2 = real(1) / real(120); 33 33 real y = sqrt(x); 34 if (!y) 35 return b2; 36 return ((sin(y) / y - a0) / x - a1) / x; 34 return sin(y) / y; 37 35 } 38 36 … … 44 42 int main(int argc, char **argv) 45 43 { 46 RemezSolver< 3> solver;44 RemezSolver<6> solver; 47 45 //solver.Run(0, 1, myfun, myfun, 15); 48 solver.Run( 0, real::R_PI * real::R_PI >> 2, myfun, myfun, 15);46 solver.Run(exp((real)-100), real::R_PI * real::R_PI >> 2, myfun, myfun, 15); 49 47 //solver.Run(-1, 1, myfun, myfun, 15); 50 48 //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15);
Note: See TracChangeset
for help on using the changeset viewer.