Changeset 1012


Ignore:
Timestamp:
Oct 6, 2011, 6:55:47 PM (10 years ago)
Author:
sam
Message:

test: more Remez exchange experimentations.

Location:
trunk/test/math
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/test/math/poly.cpp

    r1011 r1012  
    2020using namespace std;
    2121
    22 static float adjust(float f, int i)
     22static float adjustf(float f, int i)
    2323{
    2424    union { float f; uint32_t x; } u = { f };
     25    u.x += i;
     26    return u.f;
     27}
     28
     29static double adjust(double f, int i)
     30{
     31    union { double f; uint64_t x; } u = { f };
    2532    u.x += i;
    2633    return u.f;
    … …  
    3845}
    3946
     47//#define float double
     48#if 1
     49static float const a0 = adjust(0.9999999999999376, 0);
     50static float const a1 = adjust(-0.1666666666643236, 0);
     51static float const a2 = adjust(0.008333333318766562, 0);
     52static float const a3 = adjust(-0.0001984126641174625, 0);
     53static float const a4 = adjust(2.755693193297308e-006, 0);
     54static float const a5 = adjust(-2.502951900290311e-008, 0);
     55static float const a6 = adjust(1.540117123154927e-010, 0);
     56#elif 0
    4057static float const a0 = adjust(1.0, 0);
    4158static float const a1 = adjust(-0.1666666666372165, 0);
    … …  
    4461static float const a4 = adjust(2.753619853326498e-06, 0);
    4562static float const a5 = adjust(-2.407483949485896e-08, 0);
     63static float const a6 = 0.0;
     64#else
     65static float const a0 = adjust(0.9999999946887117, 0);
     66static float const a1 = adjust(-0.1666665668590824, 0);
     67static float const a2 = adjust(0.008333025160523476, 0);
     68static float const a3 = adjust(-0.0001980741944205014, 0);
     69static float const a4 = adjust(2.60190356966559e-06, 0); // -900 in floats
     70static float const a5 = 0.0;
     71static float const a6 = 0.0;
     72#endif
    4673
    4774static float floatsin(float f)
    4875{
    4976    //return lol_sin(f);
    50 //#define float double
    5177    //static float const k = (float)real::R_2_PI;
    5278
    … …  
    5480    float f2 = f * f;
    5581    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);
    5985#undef float
    6086}
    … …  
    7197    inspect(a5);
    7298
    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))
    74100    {
     101        union { float f; uint32_t x; } u = { f };
    75102        union { float f; uint32_t x; } s1 = { sinf(f) };
    76103        union { float f; uint32_t x; } s2 = { floatsin(f) };
    … …  
    78105        switch (e)
    79106        {
     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:
    80112        case 0:
    81         case 1:
    82         case 2:
    83         case 3:
    84113            error[e]++;
    85114            break;
    86115        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);
    89116            error[4]++;
    90117            break;

  • trunk/test/math/remez.cpp

    r1011 r1012  
    3232    static real const b2 = real(1) / real(120);
    3333    real y = sqrt(x);
    34     if (!y)
    35         return b2;
    36     return ((sin(y) / y - a0) / x - a1) / x;
     34    return sin(y) / y;
    3735}
    3836
    … …  
    4442int main(int argc, char **argv)
    4543{
    46     RemezSolver<3> solver;
     44    RemezSolver<6> solver;
    4745    //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);
    4947    //solver.Run(-1, 1, myfun, myfun, 15);
    5048    //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15);
Note: See TracChangeset for help on using the changeset viewer.