test: more Remez exchange experimentations.

.gitignore

@@ -60,6 +60,7 @@ test/benchsuite
 test/quad
 test/sandbox
 test/math/pi
+test/math/poly
 test/math/remez
 tools/make-font
 # Our data

test/math/Makefile.am

@@ -3,19 +3,25 @@ AM_CPPFLAGS = -I$(top_srcdir)/src
 
 all-local: $(noinst_PROGRAMS)
 	test x$(MAKE_FSELF) = xno || make_fself pi$(EXEEXT) pi.self
+	test x$(MAKE_FSELF) = xno || make_fself poly$(EXEEXT) pi.self
 	test x$(MAKE_FSELF) = xno || make_fself remez$(EXEEXT) remez.self
 
 CLEANFILES = $(noinst_PROGRAMS:%$(EXEEXT)=%.self) \
              $(noinst_PROGRAMS:%$(EXEEXT)=%.elf) \
              $(noinst_PROGRAMS:%$(EXEEXT)=%.exe)
 
-noinst_PROGRAMS = pi remez
+noinst_PROGRAMS = pi poly remez
 
 pi_SOURCES = pi.cpp
 pi_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
 pi_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@
 pi_DEPENDENCIES = $(top_builddir)/src/liblol.a
 
+poly_SOURCES = poly.cpp
+poly_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
+poly_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@
+poly_DEPENDENCIES = $(top_builddir)/src/liblol.a
+
 remez_SOURCES = remez.cpp remez-matrix.h remez-solver.h
 remez_CPPFLAGS = @LOL_CFLAGS@ @PIPI_CFLAGS@
 remez_LDFLAGS = $(top_builddir)/src/liblol.a @LOL_LIBS@ @PIPI_LIBS@

test/math/poly.cpp

@@ -0,0 +1,99 @@
+//
+// Lol Engine - Sample math program: chek trigonometric functions
+//
+// Copyright: (c) 2005-2011 Sam Hocevar <sam@hocevar.net>
+//   This program is free software; you can redistribute it and/or
+//   modify it under the terms of the Do What The Fuck You Want To
+//   Public License, Version 2, as published by Sam Hocevar. See
+//   http://sam.zoy.org/projects/COPYING.WTFPL for more details.
+//
+
+#if defined HAVE_CONFIG_H
+#   include "config.h"
+#endif
+
+#include <cstdio>
+
+#include "core.h"
+
+using namespace lol;
+using namespace std;
+
+static float adjust(float f, int i)
+{
+    union { float f; uint32_t x; } u = { f };
+    u.x += i;
+    return u.f;
+}
+
+static void inspect(float f)
+{
+    union { float f; uint32_t x; } u = { f };
+    printf("%08x %g  --  ", u.x, u.f);
+    int i = (u.x & 0x7fffffu) | 0x800000u;
+    int j = 23 - ((u.x >> 23) & 0xff) + ((1 << 7) - 1);
+    if (u.f <= 0)
+        i = -i;
+    printf("%i / 2^%i = %g\n", i, j, (float)i / (1LLu << j));
+}
+
+static float const a0 = adjust(1.0, 0);
+static float const a1 = adjust(-0.1666666666372165, 0);
+static float const a2 = adjust(0.008333332748323419, 0);
+static float const a3 = adjust(-0.0001984108245332497, 0);
+static float const a4 = adjust(2.753619853326498e-06, 0);
+static float const a5 = adjust(-2.407483949485896e-08, 0);
+
+static float floatsin(float f)
+{
+    //return lol_sin(f);
+//#define float double
+    //static float const k = (float)real::R_2_PI;
+
+    //f *= k;
+    float f2 = f * f;
+    float f4 = f2 * f2;
+    //return f * (a0 + f4 * (a2 + f4 * a4) + f2 * (a1 + f4 * (a3 + f4 * a5)));
+    //return f * (a0 + f2 * (a1 + f2 * (a2 + f2 * (a3 + f2 * (a4 + f2 * a5)))));
+    return f * (a0 + a1 * f2 + a2 * f2 * f2 + a3 * f2 * f2 * f2 + a4 * f2 * f2 * f2 * f2 + a5 * f2 * f2 * f2 * f2 * f2);
+#undef float
+}
+
+int main(void)
+{
+    int error[5] = { 0 };
+
+    inspect(a0);
+    inspect(a1);
+    inspect(a2);
+    inspect(a3);
+    inspect(a4);
+    inspect(a5);
+
+    for (float f = (float)real::R_PI_2; f > 1e-30; f = adjust(f, -1))
+    {
+        union { float f; uint32_t x; } s1 = { sinf(f) };
+        union { float f; uint32_t x; } s2 = { floatsin(f) };
+        int e = abs((int)(s1.x - s2.x));
+        switch (e)
+        {
+        case 0:
+        case 1:
+        case 2:
+        case 3:
+            error[e]++;
+            break;
+        default:
+            if (fabs((double)s1.f - (double)s2.f) > 1e-10 * fabs(s1.f))
+                printf("%16.14g: %16.14g %16.14g %08x %08x\n", f, s1.f, s2.f, s1.x, s2.x);
+            error[4]++;
+            break;
+        }
+    }
+
+    printf("exact: %i  off by 1: %i  by 2: %i  by 3: %i  error: %i\n",
+           error[0], error[1], error[2], error[3], error[4]);
+
+    return EXIT_SUCCESS;
+}
+

test/math/remez.cpp

@@ -26,10 +26,14 @@ using namespace std;
 /* The function we want to approximate */
 static real myfun(real const &x)
 {
+    static real const a0 = real::R_1;
+    static real const a1 = real(-11184811) >> 26;
+    static real const b1 = real(-1) / real(6);
+    static real const b2 = real(1) / real(120);
     real y = sqrt(x);
     if (!y)
-        return real::R_PI_2;
-    return sin(real::R_PI_2 * y) / y;
+        return b2;
+    return ((sin(y) / y - a0) / x - a1) / x;
 }
 
 static real myerr(real const &x)
@@ -39,8 +43,9 @@ static real myerr(real const &x)
 
 int main(int argc, char **argv)
 {
-    RemezSolver<4> solver;
-    solver.Run(0, 1, myfun, myfun, 15);
+    RemezSolver<3> solver;
+    //solver.Run(0, 1, myfun, myfun, 15);
+    solver.Run(0, real::R_PI * real::R_PI >> 2, myfun, myfun, 15);
     //solver.Run(-1, 1, myfun, myfun, 15);
     //solver.Run(0, real::R_PI * real::R_PI >> 4, myfun, myfun, 15);