source: trunk/test/math/poly.cpp @ 2297

Last change on this file since 2297 was 2297, checked in by sam, 8 years ago

base: implement UNUSED() macro and put it here and there.

  • Property svn:keywords set to Id
File size: 4.1 KB
Line 
1//
2// Lol Engine - Sample math program: chek trigonometric functions
3//
4// Copyright: (c) 2005-2011 Sam Hocevar <sam@hocevar.net>
5//   This program is free software; you can redistribute it and/or
6//   modify it under the terms of the Do What The Fuck You Want To
7//   Public License, Version 2, as published by Sam Hocevar. See
8//   http://www.wtfpl.net/ for more details.
9//
10
11#if defined HAVE_CONFIG_H
12#   include "config.h"
13#endif
14
15#include <cstdio>
16
17#include "core.h"
18
19using namespace lol;
20using namespace std;
21
22float adjustf(float f, int i) __attribute__((noinline));
23float adjustf(float f, int i)
24{
25    union { float f; uint32_t x; } u = { f };
26    u.x += i;
27    return u.f;
28}
29
30double adjust(double f, int i) __attribute__((noinline));
31double adjust(double f, int i)
32{
33    union { double f; uint64_t x; } u = { f };
34    u.x += i;
35    return u.f;
36}
37
38static void inspect(float f)
39{
40    union { float f; uint32_t x; } u = { f };
41    printf("%08x %14.12g  --  ", u.x, u.f);
42    int i = (u.x & 0x7fffffu) | 0x800000u;
43    int j = 23 - ((u.x >> 23) & 0xff) + ((1 << 7) - 1);
44    if (u.f <= 0)
45        i = -i;
46    printf("%i / 2^%i = %14.12g\n", i, j, (float)i / (1LLu << j));
47}
48
49//#define float double
50#if 1
51static float const a0 = 1.0;
52static float const a1 = -0.1666666666663036;
53static float const a2 = 0.008333333325075758;
54static float const a3 = -0.0001984126372689299;
55static float const a4 = 2.755533925906394e-06;
56static float const a5 = -2.476042626296988e-08;
57static float const a6 = 0.0;
58#elif 0
59static float const a0 = adjust(0.9999999999999376, 0);
60static float const a1 = adjust(-0.1666666666643236, 0);
61static float const a2 = adjust(0.008333333318766562, 0);
62static float const a3 = adjust(-0.0001984126641174625, 0);
63static float const a4 = adjust(2.755693193297308e-006, 0);
64static float const a5 = adjust(-2.502951900290311e-008, 0);
65static float const a6 = adjust(1.540117123154927e-010, 0);
66#elif 0
67static float const a0 = adjust(1.0, 0);
68static float const a1 = adjust(-0.1666666666372165, 0);
69static float const a2 = adjust(0.008333332748323419, 0);
70static float const a3 = adjust(-0.0001984108245332497, 0);
71static float const a4 = adjust(2.753619853326498e-06, 0);
72static float const a5 = adjust(-2.407483949485896e-08, 0);
73static float const a6 = 0.0;
74#else
75static float const a0 = adjust(0.9999999946887117, 0);
76static float const a1 = adjust(-0.1666665668590824, 0);
77static float const a2 = adjust(0.008333025160523476, 0);
78static float const a3 = adjust(-0.0001980741944205014, 0);
79static float const a4 = adjust(2.60190356966559e-06, 0); // -900 in floats
80static float const a5 = 0.0;
81static float const a6 = 0.0;
82#endif
83
84static float floatsin(float f)
85{
86    return lol_sin(f);
87    //static float const k = (float)real::R_2_PI();
88
89    //f *= k;
90    float f2 = f * f;
91    float f4 = f2 * f2;
92    return f * (a0 + f4 * (a2 + f4 * (a4 + f4 * a6)) + f2 * (a1 + f4 * (a3 + f4 * a5)));
93    //return f * (a0 + f2 * (a1 + f2 * (a2 + f2 * (a3 + f2 * (a4 + f2 * (a5 + f2 * a6))))));
94    //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);
95#undef float
96}
97
98int main(int argc, char **argv)
99{
100    UNUSED(argc, argv);
101
102    typedef union { float f; uint32_t x; } flint;
103
104    int error[5] = { 0 };
105
106    inspect(a0);
107    inspect(a1);
108    inspect(a2);
109    inspect(a3);
110    inspect(a4);
111    inspect(a5);
112
113    for (flint u = { (float)real::R_PI_2() }; u.f > 1e-30; u.x -= 1)
114    {
115        union { float f; uint32_t x; } s1 = { sinf(adjustf(u.f, 0)) };
116        union { float f; uint32_t x; } s2 = { floatsin(adjustf(u.f, 0)) };
117        int e = lol::abs((int)(s1.x - s2.x));
118        switch (e)
119        {
120        case 3:
121        case 2:
122        case 1:
123            if (lol::abs((double)s1.f - (double)s2.f) > 1e-10 * lol::abs(s1.f))
124                printf("%15.13g %08x: %15.13g %15.13g %08x %08x\n", u.f, u.x, s1.f, s2.f, s1.x, s2.x);
125        case 0:
126            error[e]++;
127            break;
128        default:
129            error[4]++;
130            break;
131        }
132    }
133
134    printf("exact: %i  off by 1: %i  by 2: %i  by 3: %i  error: %i\n",
135           error[0], error[1], error[2], error[3], error[4]);
136
137    return EXIT_SUCCESS;
138}
139
Note: See TracBrowser for help on using the repository browser.