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

Last change on this file since 1012 was 1012, checked in by sam, 11 years ago

test: more Remez exchange experimentations.

  • Property svn:keywords set to Id
File size: 3.7 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
22static float adjustf(float f, int i)
23{
24    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 };
32    u.x += i;
33    return u.f;
34}
35
36static void inspect(float f)
37{
38    union { float f; uint32_t x; } u = { f };
39    printf("%08x %g  --  ", u.x, u.f);
40    int i = (u.x & 0x7fffffu) | 0x800000u;
41    int j = 23 - ((u.x >> 23) & 0xff) + ((1 << 7) - 1);
42    if (u.f <= 0)
43        i = -i;
44    printf("%i / 2^%i = %g\n", i, j, (float)i / (1LLu << j));
45}
46
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
57static float const a0 = adjust(1.0, 0);
58static float const a1 = adjust(-0.1666666666372165, 0);
59static float const a2 = adjust(0.008333332748323419, 0);
60static float const a3 = adjust(-0.0001984108245332497, 0);
61static float const a4 = adjust(2.753619853326498e-06, 0);
62static 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
73
74static float floatsin(float f)
75{
76    //return lol_sin(f);
77    //static float const k = (float)real::R_2_PI;
78
79    //f *= k;
80    float f2 = f * f;
81    float f4 = 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);
85#undef float
86}
87
88int main(void)
89{
90    int error[5] = { 0 };
91
92    inspect(a0);
93    inspect(a1);
94    inspect(a2);
95    inspect(a3);
96    inspect(a4);
97    inspect(a5);
98
99    for (float f = (float)real::R_PI_2; f > 1e-30; f = adjustf(f, -1))
100    {
101        union { float f; uint32_t x; } u = { f };
102        union { float f; uint32_t x; } s1 = { sinf(f) };
103        union { float f; uint32_t x; } s2 = { floatsin(f) };
104        int e = abs((int)(s1.x - s2.x));
105        switch (e)
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:
112        case 0:
113            error[e]++;
114            break;
115        default:
116            error[4]++;
117            break;
118        }
119    }
120
121    printf("exact: %i  off by 1: %i  by 2: %i  by 3: %i  error: %i\n",
122           error[0], error[1], error[2], error[3], error[4]);
123
124    return EXIT_SUCCESS;
125}
126
Note: See TracBrowser for help on using the repository browser.