# source:trunk/test/math/poly.cpp@1011Tweet

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

test: more Remez exchange experimentations.

• Property svn:keywords set to `Id`
File size: 2.6 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
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 adjust(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 void inspect(float f)
30{
31    union { float f; uint32_t x; } u = { f };
32    printf("%08x %g  --  ", u.x, u.f);
33    int i = (u.x & 0x7fffffu) | 0x800000u;
34    int j = 23 - ((u.x >> 23) & 0xff) + ((1 << 7) - 1);
35    if (u.f <= 0)
36        i = -i;
37    printf("%i / 2^%i = %g\n", i, j, (float)i / (1LLu << j));
38}
39
40static float const a0 = adjust(1.0, 0);
41static float const a1 = adjust(-0.1666666666372165, 0);
42static float const a2 = adjust(0.008333332748323419, 0);
43static float const a3 = adjust(-0.0001984108245332497, 0);
44static float const a4 = adjust(2.753619853326498e-06, 0);
45static float const a5 = adjust(-2.407483949485896e-08, 0);
46
47static float floatsin(float f)
48{
49    //return lol_sin(f);
50//#define float double
51    //static float const k = (float)real::R_2_PI;
52
53    //f *= k;
54    float f2 = f * f;
55    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);
59#undef float
60}
61
62int main(void)
63{
64    int error[5] = { 0 };
65
66    inspect(a0);
67    inspect(a1);
68    inspect(a2);
69    inspect(a3);
70    inspect(a4);
71    inspect(a5);
72
73    for (float f = (float)real::R_PI_2; f > 1e-30; f = adjust(f, -1))
74    {
75        union { float f; uint32_t x; } s1 = { sinf(f) };
76        union { float f; uint32_t x; } s2 = { floatsin(f) };
77        int e = abs((int)(s1.x - s2.x));
78        switch (e)
79        {
80        case 0:
81        case 1:
82        case 2:
83        case 3:
84            error[e]++;
85            break;
86        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);
89            error[4]++;
90            break;
91        }
92    }
93
94    printf("exact: %i  off by 1: %i  by 2: %i  by 3: %i  error: %i\n",
95           error[0], error[1], error[2], error[3], error[4]);
96
97    return EXIT_SUCCESS;
98}
99
Note: See TracBrowser for help on using the repository browser.