From 3391a4c904d34558ded8675192540421ef700f96 Mon Sep 17 00:00:00 2001 From: Sam Hocevar Date: Thu, 6 Oct 2011 07:18:17 +0000 Subject: [PATCH] test: more Remez exchange experimentations. --- .gitignore | 1 + test/math/Makefile.am | 8 +++- test/math/poly.cpp | 99 +++++++++++++++++++++++++++++++++++++++++++ test/math/remez.cpp | 13 ++++-- 4 files changed, 116 insertions(+), 5 deletions(-) create mode 100644 test/math/poly.cpp diff --git a/.gitignore b/.gitignore index fca8d554..f0d32f94 100644 --- a/.gitignore +++ b/.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 diff --git a/test/math/Makefile.am b/test/math/Makefile.am index 735389b8..c8b6ea11 100644 --- a/test/math/Makefile.am +++ b/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@ diff --git a/test/math/poly.cpp b/test/math/poly.cpp new file mode 100644 index 00000000..1a2766bc --- /dev/null +++ b/test/math/poly.cpp @@ -0,0 +1,99 @@ +// +// Lol Engine - Sample math program: chek trigonometric functions +// +// Copyright: (c) 2005-2011 Sam Hocevar +// 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 + +#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; +} + diff --git a/test/math/remez.cpp b/test/math/remez.cpp index 70c1d358..03f12052 100644 --- a/test/math/remez.cpp +++ b/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);