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

Last change on this file since 2056 was 2056, checked in by sam, 7 years ago

math: reimplement min(), max(), abs() and fmod() in the lol:: namespace to
avoid conflicts with the C++ stdlib.

  • Property svn:keywords set to Id
File size: 4.4 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#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    typedef union { float f; uint32_t x; } flint;
101
102    int error[5] = { 0 };
103
104    inspect(a0);
105    inspect(a1);
106    inspect(a2);
107    inspect(a3);
108    inspect(a4);
109    inspect(a5);
110
111//flint v = { 1.433971524239 };
112flint v = { 1.555388212204 };
113inspect(v.f);
114printf("sinf: ");
115flint w = { sinf(adjustf(v.f, 0)) };
116inspect(w.f);
117printf("lols: ");
118flint z = { lol_sin(adjustf(v.f, 0)) };
119inspect(z.f);
120
121printf("-- START --\n");
122    for (flint u = { (float)real::R_PI_2() }; u.f > 1e-30; u.x -= 1)
123    {
124        union { float f; uint32_t x; } s1 = { sinf(adjustf(u.f, 0)) };
125        union { float f; uint32_t x; } s2 = { floatsin(adjustf(u.f, 0)) };
126        int e = lol::abs((int)(s1.x - s2.x));
127        switch (e)
128        {
129        case 3:
130        case 2:
131        case 1:
132inspect(u.f);
133printf("sinf: ");
134inspect(sinf(u.f));
135            if (lol::abs((double)s1.f - (double)s2.f) > 1e-10 * lol::abs(s1.f))
136                printf("%15.13g %08x: %15.13g %15.13g %08x %08x\n", u.f, u.x, s1.f, s2.f, s1.x, s2.x);
137        case 0:
138            error[e]++;
139            break;
140        default:
141            error[4]++;
142            break;
143        }
144    }
145
146    printf("exact: %i  off by 1: %i  by 2: %i  by 3: %i  error: %i\n",
147           error[0], error[1], error[2], error[3], error[4]);
148
149    return EXIT_SUCCESS;
150}
151
Note: See TracBrowser for help on using the repository browser.