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

Last change on this file since 1309 was 1050, checked in by sam, 9 years ago

test: fix OS X compilation; we still need SDLmain.a on that platform.

  • Property svn:keywords set to Id
File size: 4.5 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://sam.zoy.org/projects/COPYING.WTFPL for more details.
9//
10
11#if defined HAVE_CONFIG_H
12#   include "config.h"
13#endif
14
15#include <cstdio>
16
17#if USE_SDL && defined __APPLE__
18#   include <SDL_main.h>
19#endif
20
21#include "core.h"
22
23using namespace lol;
24using namespace std;
25
26float adjustf(float f, int i) __attribute__((noinline));
27float adjustf(float f, int i)
28{
29    union { float f; uint32_t x; } u = { f };
30    u.x += i;
31    return u.f;
32}
33
34double adjust(double f, int i) __attribute__((noinline));
35double adjust(double f, int i)
36{
37    union { double f; uint64_t x; } u = { f };
38    u.x += i;
39    return u.f;
40}
41
42static void inspect(float f)
43{
44    union { float f; uint32_t x; } u = { f };
45    printf("%08x %14.12g  --  ", u.x, u.f);
46    int i = (u.x & 0x7fffffu) | 0x800000u;
47    int j = 23 - ((u.x >> 23) & 0xff) + ((1 << 7) - 1);
48    if (u.f <= 0)
49        i = -i;
50    printf("%i / 2^%i = %14.12g\n", i, j, (float)i / (1LLu << j));
51}
52
53//#define float double
54#if 1
55static float const a0 = 1.0;
56static float const a1 = -0.1666666666663036;
57static float const a2 = 0.008333333325075758;
58static float const a3 = -0.0001984126372689299;
59static float const a4 = 2.755533925906394e-06;
60static float const a5 = -2.476042626296988e-08;
61static float const a6 = 0.0;
62#elif 0
63static float const a0 = adjust(0.9999999999999376, 0);
64static float const a1 = adjust(-0.1666666666643236, 0);
65static float const a2 = adjust(0.008333333318766562, 0);
66static float const a3 = adjust(-0.0001984126641174625, 0);
67static float const a4 = adjust(2.755693193297308e-006, 0);
68static float const a5 = adjust(-2.502951900290311e-008, 0);
69static float const a6 = adjust(1.540117123154927e-010, 0);
70#elif 0
71static float const a0 = adjust(1.0, 0);
72static float const a1 = adjust(-0.1666666666372165, 0);
73static float const a2 = adjust(0.008333332748323419, 0);
74static float const a3 = adjust(-0.0001984108245332497, 0);
75static float const a4 = adjust(2.753619853326498e-06, 0);
76static float const a5 = adjust(-2.407483949485896e-08, 0);
77static float const a6 = 0.0;
78#else
79static float const a0 = adjust(0.9999999946887117, 0);
80static float const a1 = adjust(-0.1666665668590824, 0);
81static float const a2 = adjust(0.008333025160523476, 0);
82static float const a3 = adjust(-0.0001980741944205014, 0);
83static float const a4 = adjust(2.60190356966559e-06, 0); // -900 in floats
84static float const a5 = 0.0;
85static float const a6 = 0.0;
86#endif
87
88static float floatsin(float f)
89{
90    return lol_sin(f);
91    //static float const k = (float)real::R_2_PI;
92
93    //f *= k;
94    float f2 = f * f;
95    float f4 = f2 * f2;
96    return f * (a0 + f4 * (a2 + f4 * (a4 + f4 * a6)) + f2 * (a1 + f4 * (a3 + f4 * a5)));
97    //return f * (a0 + f2 * (a1 + f2 * (a2 + f2 * (a3 + f2 * (a4 + f2 * (a5 + f2 * a6))))));
98    //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);
99#undef float
100}
101
102int main(int argc, char **argv)
103{
104    typedef union { float f; uint32_t x; } flint;
105
106    int error[5] = { 0 };
107
108    inspect(a0);
109    inspect(a1);
110    inspect(a2);
111    inspect(a3);
112    inspect(a4);
113    inspect(a5);
114
115//flint v = { 1.433971524239 };
116flint v = { 1.555388212204 };
117inspect(v.f);
118printf("sinf: ");
119flint w = { sinf(adjustf(v.f, 0)) };
120inspect(w.f);
121printf("lols: ");
122flint z = { lol_sin(adjustf(v.f, 0)) };
123inspect(z.f);
124
125printf("-- START --\n");
126    for (flint u = { (float)real::R_PI_2 }; u.f > 1e-30; u.x -= 1)
127    {
128        union { float f; uint32_t x; } s1 = { sinf(adjustf(u.f, 0)) };
129        union { float f; uint32_t x; } s2 = { floatsin(adjustf(u.f, 0)) };
130        int e = abs((int)(s1.x - s2.x));
131        switch (e)
132        {
133        case 3:
134        case 2:
135        case 1:
136inspect(u.f);
137printf("sinf: ");
138inspect(sinf(u.f));
139            if (fabs((double)s1.f - (double)s2.f) > 1e-10 * fabs(s1.f))
140                printf("%15.13g %08x: %15.13g %15.13g %08x %08x\n", u.f, u.x, s1.f, s2.f, s1.x, s2.x);
141        case 0:
142            error[e]++;
143            break;
144        default:
145            error[4]++;
146            break;
147        }
148    }
149
150    printf("exact: %i  off by 1: %i  by 2: %i  by 3: %i  error: %i\n",
151           error[0], error[1], error[2], error[3], error[4]);
152
153    return EXIT_SUCCESS;
154}
155
Note: See TracBrowser for help on using the repository browser.