diff --git a/build/Lol (vs2015).sln b/build/Lol (vs2015).sln index d24d0a1d..7ecd0c1a 100644 --- a/build/Lol (vs2015).sln +++ b/build/Lol (vs2015).sln @@ -58,10 +58,6 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Test", "Test", "{E4DFEBF9-C EndProject Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "13_shader_builder", "..\doc\tutorial\13_shader_builder.vcxproj", "{F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}" EndProject -Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "LolRemez", "LolRemez", "{4C4BD478-3767-4C27-BD91-DAAFE7CD03A2}" -EndProject -Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "lolremez", "..\tools\lolremez\lolremez.vcxproj", "{73F1A804-1116-46C3-922A-9C0ADEB33F52}" -EndProject Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Samples", "Samples", "{B6297FF2-63D0-41EE-BE13-EFF720C9B0FA}" EndProject Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "15_lolimgui", "..\doc\tutorial\15_lolimgui.vcxproj", "{81C83B42-D00A-4FA3-9A3D-80F9D46524BF}" @@ -286,18 +282,6 @@ Global {F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}.Release|Win32.Build.0 = Release|Win32 {F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}.Release|x64.ActiveCfg = Release|x64 {F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}.Release|x64.Build.0 = Release|x64 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Debug|ORBIS.ActiveCfg = Debug|ORBIS - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Debug|ORBIS.Build.0 = Debug|ORBIS - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Debug|Win32.ActiveCfg = Debug|Win32 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Debug|Win32.Build.0 = Debug|Win32 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Debug|x64.ActiveCfg = Debug|x64 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Debug|x64.Build.0 = Debug|x64 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Release|ORBIS.ActiveCfg = Release|ORBIS - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Release|ORBIS.Build.0 = Release|ORBIS - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Release|Win32.ActiveCfg = Release|Win32 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Release|Win32.Build.0 = Release|Win32 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Release|x64.ActiveCfg = Release|x64 - {73F1A804-1116-46C3-922A-9C0ADEB33F52}.Release|x64.Build.0 = Release|x64 {81C83B42-D00A-4FA3-9A3D-80F9D46524BF}.Debug|ORBIS.ActiveCfg = Debug|ORBIS {81C83B42-D00A-4FA3-9A3D-80F9D46524BF}.Debug|ORBIS.Build.0 = Debug|ORBIS {81C83B42-D00A-4FA3-9A3D-80F9D46524BF}.Debug|Win32.ActiveCfg = Debug|Win32 @@ -351,8 +335,6 @@ Global {E05E23A5-67DE-42B5-98A3-E63CCE0CC0AF} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} {E4DFEBF9-C310-462F-9876-7EB59C1E4D4E} = {1AFD580B-98B8-4689-B661-38C41132C60E} {F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} - {4C4BD478-3767-4C27-BD91-DAAFE7CD03A2} = {3D341D8A-E400-4B1D-BC05-B5C35487D9B5} - {73F1A804-1116-46C3-922A-9C0ADEB33F52} = {4C4BD478-3767-4C27-BD91-DAAFE7CD03A2} {B6297FF2-63D0-41EE-BE13-EFF720C9B0FA} = {1AFD580B-98B8-4689-B661-38C41132C60E} {81C83B42-D00A-4FA3-9A3D-80F9D46524BF} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} {31B96262-1C41-43B9-BA38-27AA385B05DB} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} diff --git a/configure.ac b/configure.ac index 891738fd..4ffd0fab 100644 --- a/configure.ac +++ b/configure.ac @@ -21,8 +21,6 @@ AM_DEFAULT_VERBOSITY=0 dnl Versioning of the separate software we ship LOLUNIT_VERSION=0.1 AC_SUBST(LOLUNIT_VERSION) -LOLREMEZ_VERSION=0.2 -AC_SUBST(LOLREMEZ_VERSION) AC_SUBST(lol_srcdir, '${top_srcdir}') AC_SUBST(lol_builddir, '${top_builddir}') @@ -262,7 +260,6 @@ AC_CONFIG_FILES( doc/samples/sandbox/Makefile doc/tutorial/Makefile tools/Makefile - tools/lolremez/Makefile tools/lolunit/Makefile tools/vimlol/Makefile tools/vslol/Makefile diff --git a/tools/Makefile.am b/tools/Makefile.am index 80b16881..296c66f5 100644 --- a/tools/Makefile.am +++ b/tools/Makefile.am @@ -2,7 +2,6 @@ include $(top_srcdir)/build/autotools/common.am SUBDIRS = -SUBDIRS += lolremez SUBDIRS += lolunit SUBDIRS += vimlol SUBDIRS += vslol diff --git a/tools/lolremez/.gitignore b/tools/lolremez/.gitignore deleted file mode 100644 index cabca74d..00000000 --- a/tools/lolremez/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -# Our binaries -lolremez diff --git a/tools/lolremez/Makefile.am b/tools/lolremez/Makefile.am deleted file mode 100644 index e07526b6..00000000 --- a/tools/lolremez/Makefile.am +++ /dev/null @@ -1,14 +0,0 @@ - -include $(top_srcdir)/build/autotools/common.am - -EXTRA_DIST += NEWS.txt lolremez.sln lolremez.vcxproj lolremez.vcxproj.filters - -if BUILD_TOOLS -noinst_PROGRAMS = lolremez -endif - -lolremez_SOURCES = \ - lolremez.cpp solver.cpp solver.h matrix.h expression.h -lolremez_CPPFLAGS = $(AM_CPPFLAGS) -lolremez_DEPENDENCIES = @LOL_DEPS@ - diff --git a/tools/lolremez/NEWS.txt b/tools/lolremez/NEWS.txt deleted file mode 100644 index 35f7cbdc..00000000 --- a/tools/lolremez/NEWS.txt +++ /dev/null @@ -1,13 +0,0 @@ - -News for LolRemez 0.2: - - - significant performance and accuracy improvements thanks to various - bugfixes and a better extrema finder for the error function. - - user can now define accuracy of the final result. - - exp(), sin(), cos() and tan() are now about 20% faster. - - multiplying a real number by an integer power of two is now a virtually - free operation. - - fixed a rounding bug in the real number printing routine. - -Initial release: LolRemez 0.1 - diff --git a/tools/lolremez/expression.h b/tools/lolremez/expression.h deleted file mode 100644 index 25fe9741..00000000 --- a/tools/lolremez/expression.h +++ /dev/null @@ -1,349 +0,0 @@ -// -// LolRemez — Remez algorithm implementation -// -// Copyright © 2005—2015 Sam Hocevar -// -// This program is free software. It comes without any warranty, to -// the extent permitted by applicable law. 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 the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// Powerful arithmetic expression parser/evaluator -// -// Usage: -// expression e; -// e.parse(" 2*x^3 + 3 * sin(x - atan(x))"); -// auto y = e.eval("1.5"); -// - -#include "pegtl.hh" - -namespace grammar -{ - -using namespace pegtl; - -enum class id : uint8_t -{ - /* Variables and constants */ - x, - constant, - /* Unary functions/operators */ - plus, minus, abs, - sqrt, cbrt, - exp, exp2, log, log2, log10, - sin, cos, tan, - asin, acos, atan, - sinh, cosh, tanh, - /* Binary functions/operators */ - add, sub, mul, div, - atan2, pow, - min, max, -}; - -struct expression -{ - /* - * Evaluate expression at x - */ - lol::real eval(lol::real const &x) const - { - /* Use a stack */ - lol::array stack; - - for (int i = 0; i < m_ops.count(); ++i) - { - /* Rules that do not consume stack elements */ - if (m_ops[i].m1 == id::x) - { - stack.push(x); - continue; - } - else if (m_ops[i].m1 == id::constant) - { - stack.push(m_constants[m_ops[i].m2]); - continue; - } - - /* All other rules consume at least the head of the stack */ - lol::real head = stack.pop(); - - switch (m_ops[i].m1) - { - case id::plus: stack.push(head); break; - case id::minus: stack.push(-head); break; - - case id::abs: stack.push(fabs(head)); break; - case id::sqrt: stack.push(sqrt(head)); break; - case id::cbrt: stack.push(cbrt(head)); break; - case id::exp: stack.push(exp(head)); break; - case id::exp2: stack.push(exp2(head)); break; - case id::log: stack.push(log(head)); break; - case id::log2: stack.push(log2(head)); break; - case id::log10: stack.push(log10(head)); break; - case id::sin: stack.push(sin(head)); break; - case id::cos: stack.push(cos(head)); break; - case id::tan: stack.push(tan(head)); break; - case id::asin: stack.push(asin(head)); break; - case id::acos: stack.push(acos(head)); break; - case id::atan: stack.push(atan(head)); break; - case id::sinh: stack.push(sinh(head)); break; - case id::cosh: stack.push(cosh(head)); break; - case id::tanh: stack.push(tanh(head)); break; - - case id::add: stack.push(stack.pop() + head); break; - case id::sub: stack.push(stack.pop() - head); break; - case id::mul: stack.push(stack.pop() * head); break; - case id::div: stack.push(stack.pop() / head); break; - - case id::atan2: stack.push(atan2(stack.pop(), head)); break; - case id::pow: stack.push(pow(stack.pop(), head)); break; - case id::min: stack.push(min(stack.pop(), head)); break; - case id::max: stack.push(max(stack.pop(), head)); break; - - case id::x: - case id::constant: - /* Already handled above */ - break; - } - } - - ASSERT(stack.count() == 1); - return stack.pop(); - } - -private: - lol::array m_ops; - lol::array m_constants; - -private: - struct r_expr; - - // r_ <- * - struct _ : star {}; - - // r_constant <- + ( "." * ) ? ( [eE] [+-] ? + ) ? - struct r_constant : seq, - opt, - star>>, - opt, - opt>, - plus>>> {}; - - // r_var <- "x" - struct r_var : pegtl_string_t("x") {}; - - // r_binary_call <- "(" r_expr "," r_expr ")" - struct r_binary_fun : sor {}; - - struct r_binary_call : seq, - _, r_expr, - _, one<','>, - _, r_expr, - _, one<')'>> {}; - - // r_unary_call <- "(" r_expr ")" - struct r_unary_fun : sor {}; - - struct r_unary_call : seq, - _, r_expr, - _, one<')'>> {}; - - // r_call <- r_binary_call / r_unary_call - struct r_call : sor {}; - - // r_parentheses <- "(" r_expr ")" - struct r_parentheses : seq, - pad, - one<')'>> {}; - - // r_terminal <- r_call / r_var / r_constant / r_parentheses - struct r_terminal : sor {}; - - // r_signed <- "-" r_signed / "+" r_signed / r_terminal - struct r_signed; - struct r_minus : seq, _, r_signed> {}; - struct r_signed : sor, _, r_signed>, - r_terminal> {}; - - // r_exponent <- ( "^" / "**" ) r_signed - struct r_exponent : seq, - string<'*', '*'>>, space>, - r_signed> {}; - - // r_factor <- r_signed ( r_exponent ) * - struct r_factor : seq> {}; - - // r_mul <- "*" r_factor - // r_div <- "/" r_factor - // r_term <- r_factor ( r_mul / r_div ) * - struct r_mul : seq<_, one<'*'>, _, r_factor> {}; - struct r_div : seq<_, one<'/'>, _, r_factor> {}; - struct r_term : seq>> {}; - - // r_add <- "+" r_term - // r_sub <- "-" r_term - // r_expr <- r_term ( r_add / r_sub ) * - struct r_add : seq<_, one<'+'>, _, r_term> {}; - struct r_sub : seq<_, one<'-'>, _, r_term> {}; - struct r_expr : seq>> {}; - - // r_stmt <- r_expr - struct r_stmt : must, pegtl::eof> {}; - - // - // Default actions - // - - template - struct action : nothing {}; - - template - struct generic_action - { - static void apply(action_input const &in, expression *that) - { - UNUSED(in); - that->m_ops.push(OP, -1); - } - }; - -public: - /* - * Parse arithmetic expression in x, e.g. 2*x+3 - */ - void parse(std::string const &str) - { - m_ops.empty(); - m_constants.empty(); - - pegtl::parse_string(str, "expression", this); - } -}; - -// -// Rule specialisations for simple operators -// - -template<> struct expression::action : generic_action {}; -template<> struct expression::action : generic_action {}; -template<> struct expression::action : generic_action {}; -template<> struct expression::action : generic_action {}; -template<> struct expression::action : generic_action {}; -template<> struct expression::action : generic_action {}; -template<> struct expression::action : generic_action {}; - -// -// Rule specialisations for unary and binary function calls -// - -template<> -struct expression::action -{ - static void apply(action_input const &in, expression *that) - { - struct { id ret; char const *name; } lut[] = - { - { id::atan2, "atan2" }, - { id::pow, "pow" }, - { id::min, "min" }, - { id::max, "max" }, - }; - - for (auto pair : lut) - { - if (strncmp(in.string().c_str(), pair.name, strlen(pair.name)) != 0) - continue; - - that->m_ops.push(pair.ret, -1); - return; - } - } -}; - -template<> -struct expression::action -{ - static void apply(action_input const &in, expression *that) - { - struct { id ret; char const *name; } lut[] = - { - { id::abs, "abs" }, - { id::sqrt, "sqrt" }, - { id::cbrt, "cbrt" }, - { id::exp2, "exp2" }, - { id::exp, "exp" }, - { id::log10, "log10" }, - { id::log2, "log2" }, - { id::log, "log" }, - { id::sinh, "sinh" }, - { id::cosh, "cosh" }, - { id::tanh, "tanh" }, - { id::sin, "sin" }, - { id::cos, "cos" }, - { id::tan, "tan" }, - { id::asin, "asin" }, - { id::acos, "acos" }, - { id::atan, "atan" }, - }; - - for (auto pair : lut) - { - if (strncmp(in.string().c_str(), pair.name, strlen(pair.name)) != 0) - continue; - - that->m_ops.push(pair.ret, -1); - return; - } - } -}; - -template<> -struct expression::action -{ - static void apply(action_input const &in, expression *that) - { - /* FIXME: check if the constant is already in the list */ - that->m_ops.push(id::constant, that->m_constants.count()); - that->m_constants.push(lol::real(in.string().c_str())); - } -}; - -} /* namespace grammar */ - -using grammar::expression; - diff --git a/tools/lolremez/lolremez.cpp b/tools/lolremez/lolremez.cpp deleted file mode 100644 index 03d49bc3..00000000 --- a/tools/lolremez/lolremez.cpp +++ /dev/null @@ -1,129 +0,0 @@ -// -// LolRemez - Remez algorithm implementation -// -// Copyright © 2005—2016 Sam Hocevar -// -// This program is free software. It comes without any warranty, to -// the extent permitted by applicable law. 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 the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#if HAVE_CONFIG_H -# include "config.h" -#endif - -#include - -#include - -#include "solver.h" -#include "expression.h" - -using lol::array; -using lol::real; -using lol::String; - -static void version(void) -{ - printf("lolremez %s\n", PACKAGE_VERSION); - printf("Copyright © 2005—2016 Sam Hocevar \n"); - printf("This program is free software. It comes without any warranty, to the extent\n"); - printf("permitted by applicable law. You can redistribute it and/or modify it under\n"); - printf("the terms of the Do What the Fuck You Want to Public License, Version 2, as\n"); - printf("published by the WTFPL Task Force. See http://www.wtfpl.net/ for more details.\n"); - printf("\n"); - printf("Written by Sam Hocevar. Report bugs to .\n"); -} - -static void usage() -{ - printf("Usage: lolremez [-d degree] [-r xmin:xmax] x-expression [x-error]\n"); - printf(" lolremez -h | --help\n"); - printf(" lolremez -V | --version\n"); - printf("Find a polynomial approximation for x-expression.\n"); - printf("\n"); - printf("Mandatory arguments to long options are mandatory for short options too.\n"); - printf(" -d, --degree degree of final polynomial\n"); - printf(" -r, --range : range over which to approximate\n"); - printf(" -h, --help display this help and exit\n"); - printf(" -V, --version output version information and exit\n"); - printf("\n"); - printf("Examples:\n"); - printf(" lolremez -d 4 -r -1:1 \"atan(exp(1+x))\"\n"); - printf(" lolremez -d 4 -r -1:1 \"atan(exp(1+x))\" \"exp(1+x)\"\n"); - printf("\n"); - printf("Written by Sam Hocevar. Report bugs to .\n"); -} - -static void FAIL(char const *message = nullptr) -{ - if (message) - printf("Error: %s\n", message); - printf("Try 'lolremez --help' for more information.\n"); - exit(EXIT_FAILURE); -} - -/* See the tutorial at http://lolengine.net/wiki/doc/maths/remez */ -int main(int argc, char **argv) -{ - String xmin("-1"), xmax("1"); - char const *f = nullptr, *g = nullptr; - int degree = 4; - - lol::getopt opt(argc, argv); - opt.add_opt('h', "help", false); - opt.add_opt('v', "version", false); - opt.add_opt('d', "degree", true); - opt.add_opt('r', "range", true); - - for (;;) - { - int c = opt.parse(); - if (c == -1) - break; - - switch (c) - { - case 'd': /* --degree */ - degree = atoi(opt.arg); - break; - case 'r': { /* --range */ - array arg = String(opt.arg).split(':'); - if (arg.count() != 2) - FAIL("invalid range"); - xmin = arg[0]; - xmax = arg[1]; - } break; - case 'h': /* --help */ - usage(); - return EXIT_SUCCESS; - case 'v': /* --version */ - version(); - return EXIT_SUCCESS; - default: - FAIL(); - } - } - - if (opt.index < argc) - f = argv[opt.index++]; - - if (opt.index < argc) - g = argv[opt.index++]; - - if (!f) - FAIL("no function specified"); - else if (opt.index < argc) - FAIL("too many arguments"); - - if (real(xmin.C()) >= real(xmax.C())) - FAIL("invalid range"); - - remez_solver solver(degree, 20); - solver.run(xmin.C(), xmax.C(), f, g); - - return 0; -} - diff --git a/tools/lolremez/lolremez.sln b/tools/lolremez/lolremez.sln deleted file mode 100644 index ba6b3637..00000000 --- a/tools/lolremez/lolremez.sln +++ /dev/null @@ -1,31 +0,0 @@ - -Microsoft Visual Studio Solution File, Format Version 11.00 -# Visual Studio 2010 -Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "lolremez", "lolremez.vcxproj", "{A2767731-6FBF-442C-900F-301D095516E7}" -EndProject -Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution Items", "{663F5D32-D212-4A88-A681-84B54D3188A3}" - ProjectSection(SolutionItems) = preProject - Performance1.psess = Performance1.psess - EndProjectSection -EndProject -Global - GlobalSection(SolutionConfigurationPlatforms) = preSolution - Debug|Win32 = Debug|Win32 - Debug|x64 = Debug|x64 - Release|Win32 = Release|Win32 - Release|x64 = Release|x64 - EndGlobalSection - GlobalSection(ProjectConfigurationPlatforms) = postSolution - {A2767731-6FBF-442C-900F-301D095516E7}.Debug|Win32.ActiveCfg = Debug|Win32 - {A2767731-6FBF-442C-900F-301D095516E7}.Debug|Win32.Build.0 = Debug|Win32 - {A2767731-6FBF-442C-900F-301D095516E7}.Debug|x64.ActiveCfg = Debug|x64 - {A2767731-6FBF-442C-900F-301D095516E7}.Debug|x64.Build.0 = Debug|x64 - {A2767731-6FBF-442C-900F-301D095516E7}.Release|Win32.ActiveCfg = Release|Win32 - {A2767731-6FBF-442C-900F-301D095516E7}.Release|Win32.Build.0 = Release|Win32 - {A2767731-6FBF-442C-900F-301D095516E7}.Release|x64.ActiveCfg = Release|x64 - {A2767731-6FBF-442C-900F-301D095516E7}.Release|x64.Build.0 = Release|x64 - EndGlobalSection - GlobalSection(SolutionProperties) = preSolution - HideSolutionNode = FALSE - EndGlobalSection -EndGlobal diff --git a/tools/lolremez/lolremez.vcxproj b/tools/lolremez/lolremez.vcxproj deleted file mode 100644 index 0cfd3598..00000000 --- a/tools/lolremez/lolremez.vcxproj +++ /dev/null @@ -1,73 +0,0 @@ - - - - $(SolutionDir)\lol - $(SolutionDir)\.. - - - - Debug - ORBIS - - - Debug - Win32 - - - Debug - x64 - - - Release - ORBIS - - - Release - Win32 - - - Release - x64 - - - - - - - - - - - - - - {9e62f2fe-3408-4eae-8238-fd84238ceeda} - - - {83d3b207-c601-4025-8f41-01dedc354661} - - - {d84021ca-b233-4e0f-8a52-071b83bbccc4} - - - - {73F1A804-1116-46C3-922A-9C0ADEB33F52} - Application - Win32Proj - - - - - - - - - - - - - - - - - diff --git a/tools/lolremez/lolremez.vcxproj.filters b/tools/lolremez/lolremez.vcxproj.filters deleted file mode 100644 index a3c70a93..00000000 --- a/tools/lolremez/lolremez.vcxproj.filters +++ /dev/null @@ -1,12 +0,0 @@ - - - - - - - - - - - - \ No newline at end of file diff --git a/tools/lolremez/matrix.h b/tools/lolremez/matrix.h deleted file mode 100644 index e326cdd8..00000000 --- a/tools/lolremez/matrix.h +++ /dev/null @@ -1,97 +0,0 @@ -// -// LolRemez - Remez algorithm implementation -// -// Copyright © 2005—2015 Sam Hocevar -// -// This program is free software. It comes without any warranty, to -// the extent permitted by applicable law. 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 the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#pragma once - -using namespace lol; - -/* - * Arbitrarily-sized square matrices; for now this only supports - * naive inversion and is used for the Remez inversion method. - */ - -template -struct linear_system : public array2d -{ - inline linear_system(int cols) - { - ASSERT(cols > 0); - - this->resize(ivec2(cols)); - } - - void init(T const &x) - { - int const n = this->size().x; - - for (int j = 0; j < n; j++) - for (int i = 0; i < n; i++) - (*this)[i][j] = (i == j) ? x : (T)0; - } - - /* Naive matrix inversion */ - linear_system inverse() const - { - int const n = this->size().x; - linear_system a(*this), b(n); - - b.init((T)1); - - /* Inversion method: iterate through all columns and make sure - * all the terms are 1 on the diagonal and 0 everywhere else */ - for (int i = 0; i < n; i++) - { - /* If the expected coefficient is zero, add one of - * the other lines. The first we meet will do. */ - if (!a[i][i]) - { - for (int j = i + 1; j < n; j++) - { - if (!a[i][j]) - continue; - /* Add row j to row i */ - for (int k = 0; k < n; k++) - { - a[k][i] += a[k][j]; - b[k][i] += b[k][j]; - } - break; - } - } - - /* Now we know the diagonal term is non-zero. Get its inverse - * and use that to nullify all other terms in the column */ - T x = (T)1 / a[i][i]; - for (int j = 0; j < n; j++) - { - if (j == i) - continue; - T mul = x * a[i][j]; - for (int k = 0; k < n; k++) - { - a[k][j] -= mul * a[k][i]; - b[k][j] -= mul * b[k][i]; - } - } - - /* Finally, ensure the diagonal term is 1 */ - for (int k = 0; k < n; k++) - { - a[k][i] *= x; - b[k][i] *= x; - } - } - - return b; - } -}; - diff --git a/tools/lolremez/solver.cpp b/tools/lolremez/solver.cpp deleted file mode 100644 index 15b2cf8f..00000000 --- a/tools/lolremez/solver.cpp +++ /dev/null @@ -1,419 +0,0 @@ -// -// LolRemez - Remez algorithm implementation -// -// Copyright © 2005—2015 Sam Hocevar -// -// This program is free software. It comes without any warranty, to -// the extent permitted by applicable law. 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 the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#if HAVE_CONFIG_H -# include "config.h" -#endif - -#include - -#include - -#include -#include - -#include "matrix.h" -#include "solver.h" - -using lol::real; - -remez_solver::remez_solver(int order, int decimals) - : m_order(order), - m_decimals(decimals), - m_has_weight(false) -{ - /* Spawn 4 worker threads */ - for (int i = 0; i < 4; ++i) - { - auto th = new thread(std::bind(&remez_solver::worker_thread, this)); - m_workers.push(th); - } -} - -remez_solver::~remez_solver() -{ - /* Signal worker threads to quit, wait for worker threads to answer, - * and kill worker threads. */ - for (auto worker : m_workers) - UNUSED(worker), m_questions.push(-1); - - for (auto worker : m_workers) - UNUSED(worker), m_answers.pop(); - - for (auto worker : m_workers) - delete worker; -} - -void remez_solver::run(real a, real b, char const *func, char const *weight) -{ - m_func.parse(func); - - if (weight) - { - m_weight.parse(weight); - m_has_weight = true; - } - - m_k1 = (b + a) / 2; - m_k2 = (b - a) / 2; - m_epsilon = pow((real)10, (real)-(m_decimals + 2)); - - remez_init(); - print_poly(); - - for (int n = 0; ; n++) - { - real old_error = m_error; - - find_extrema(); - remez_step(); - - if (m_error >= (real)0 - && fabs(m_error - old_error) < m_error * m_epsilon) - break; - - print_poly(); - find_zeroes(); - } - - print_poly(); -} - -/* - * This is basically the first Remez step: we solve a system of - * order N+1 and get a good initial polynomial estimate. - */ -void remez_solver::remez_init() -{ - /* m_order + 1 zeroes of the error function */ - m_zeroes.resize(m_order + 1); - - /* m_order + 1 zeroes to find */ - m_zeroes_state.resize(m_order + 1); - - /* m_order + 2 control points */ - m_control.resize(m_order + 2); - - /* m_order extrema to find */ - m_extrema_state.resize(m_order); - - /* Initial estimates for the x_i where the error will be zero and - * precompute f(x_i). */ - array fxn; - for (int i = 0; i < m_order + 1; i++) - { - m_zeroes[i] = (real)(2 * i - m_order) / (real)(m_order + 1); - fxn.push(eval_func(m_zeroes[i])); - } - - /* We build a matrix of Chebyshev evaluations: row i contains the - * evaluations of x_i for polynomial order n = 0, 1, ... */ - linear_system system(m_order + 1); - for (int n = 0; n < m_order + 1; n++) - { - auto p = polynomial::chebyshev(n); - for (int i = 0; i < m_order + 1; i++) - system[i][n] = p.eval(m_zeroes[i]); - } - - /* Solve the system */ - system = system.inverse(); - - /* Compute new Chebyshev estimate */ - m_estimate = polynomial(); - for (int n = 0; n < m_order + 1; n++) - { - real weight = 0; - for (int i = 0; i < m_order + 1; i++) - weight += system[n][i] * fxn[i]; - - m_estimate += weight * polynomial::chebyshev(n); - } -} - -/* - * Every subsequent iteration of the Remez algorithm: we solve a system - * of order N+2 to both refine the estimate and compute the error. - */ -void remez_solver::remez_step() -{ - Timer t; - - /* Pick up x_i where error will be 0 and compute f(x_i) */ - array fxn; - for (int i = 0; i < m_order + 2; i++) - fxn.push(eval_func(m_control[i])); - - /* We build a matrix of Chebyshev evaluations: row i contains the - * evaluations of x_i for polynomial order n = 0, 1, ... */ - linear_system system(m_order + 2); - for (int n = 0; n < m_order + 1; n++) - { - auto p = polynomial::chebyshev(n); - for (int i = 0; i < m_order + 2; i++) - system[i][n] = p.eval(m_control[i]); - } - - /* The last line of the system is the oscillating error */ - for (int i = 0; i < m_order + 2; i++) - { - real error = fabs(eval_weight(m_control[i])); - system[i][m_order + 1] = (i & 1) ? error : -error; - } - - /* Solve the system */ - system = system.inverse(); - - /* Compute new polynomial estimate */ - m_estimate = polynomial(); - for (int n = 0; n < m_order + 1; n++) - { - real weight = 0; - for (int i = 0; i < m_order + 2; i++) - weight += system[n][i] * fxn[i]; - - m_estimate += weight * polynomial::chebyshev(n); - } - - /* Compute the error (FIXME: unused?) */ - real error = 0; - for (int i = 0; i < m_order + 2; i++) - error += system[m_order + 1][i] * fxn[i]; - - using std::printf; - printf(" -:- timing for inversion: %f ms\n", t.Get() * 1000.f); -} - -/* - * Find m_order + 1 zeroes of the error function. No need to compute the - * relative error: its zeroes are at the same place as the absolute error! - * - * The algorithm used here is naïve regula falsi. It still performs a lot - * better than the midpoint algorithm. - */ -void remez_solver::find_zeroes() -{ - Timer t; - - static real const limit = ldexp((real)1, -500); - static real const zero = (real)0; - - /* Initialise an [a,b] bracket for each zero we try to find */ - for (int i = 0; i < m_order + 1; i++) - { - point &a = m_zeroes_state[i].m1; - point &b = m_zeroes_state[i].m2; - - a.x = m_control[i]; - a.err = eval_estimate(a.x) - eval_func(a.x); - b.x = m_control[i + 1]; - b.err = eval_estimate(b.x) - eval_func(b.x); - - m_questions.push(i); - } - - /* Watch all brackets for updates from worker threads */ - for (int finished = 0; finished < m_order + 1; ) - { - int i = m_answers.pop(); - - point &a = m_zeroes_state[i].m1; - point &b = m_zeroes_state[i].m2; - point &c = m_zeroes_state[i].m3; - - if (c.err == zero || fabs(a.x - b.x) <= limit) - { - m_zeroes[i] = c.x; - ++finished; - continue; - } - - m_questions.push(i); - } - - using std::printf; - printf(" -:- timing for zeroes: %f ms\n", t.Get() * 1000.f); -} - -/* - * Find m_order extrema of the error function. We maximise the relative - * error, since its extrema are at slightly different locations than the - * absolute error’s. - * - * The algorithm used here is successive parabolic interpolation. FIXME: we - * could use Brent’s method instead, which combines parabolic interpolation - * and golden ratio search and has superlinear convergence. - */ -void remez_solver::find_extrema() -{ - Timer t; - - m_control[0] = -1; - m_control[m_order + 1] = 1; - m_error = 0; - - /* Initialise an [a,b,c] bracket for each extremum we try to find */ - for (int i = 0; i < m_order; i++) - { - point &a = m_extrema_state[i].m1; - point &b = m_extrema_state[i].m2; - point &c = m_extrema_state[i].m3; - - a.x = m_zeroes[i]; - b.x = m_zeroes[i + 1]; - c.x = a.x + (b.x - a.x) * real(rand(0.4f, 0.6f)); - - a.err = eval_error(a.x); - b.err = eval_error(b.x); - c.err = eval_error(c.x); - - m_questions.push(i + 1000); - } - - /* Watch all brackets for updates from worker threads */ - for (int finished = 0; finished < m_order; ) - { - int i = m_answers.pop(); - - point &a = m_extrema_state[i - 1000].m1; - point &b = m_extrema_state[i - 1000].m2; - point &c = m_extrema_state[i - 1000].m3; - - if (b.x - a.x <= m_epsilon) - { - m_control[i - 1000 + 1] = c.x; - if (c.err > m_error) - m_error = c.err; - ++finished; - continue; - } - - m_questions.push(i); - } - - using std::printf; - printf(" -:- timing for extrema: %f ms\n", t.Get() * 1000.f); - printf(" -:- error: "); - m_error.print(m_decimals); - printf("\n"); -} - -void remez_solver::print_poly() -{ - /* Transform our polynomial in the [-1..1] range into a polynomial - * in the [a..b] range by composing it with q: - * q(x) = 2x / (b-a) - (b+a) / (b-a) */ - polynomial q ({ -m_k1 / m_k2, real(1) / m_k2 }); - polynomial r = m_estimate.eval(q); - - using std::printf; - printf("\n"); - for (int j = 0; j < m_order + 1; j++) - { - if (j) - printf(" + x**%i * ", j); - r[j].print(m_decimals); - } - printf("\n\n"); -} - -real remez_solver::eval_estimate(real const &x) -{ - return m_estimate.eval(x); -} - -real remez_solver::eval_func(real const &x) -{ - return m_func.eval(x * m_k2 + m_k1); -} - -real remez_solver::eval_weight(real const &x) -{ - return m_has_weight ? m_weight.eval(x * m_k2 + m_k1) : real(1); -} - -real remez_solver::eval_error(real const &x) -{ - return fabs((eval_estimate(x) - eval_func(x)) / eval_weight(x)); -} - -void remez_solver::worker_thread() -{ - static real const zero = (real)0; - - for (;;) - { - int i = m_questions.pop(); - - if (i < 0) - { - m_answers.push(i); - break; - } - else if (i < 1000) - { - point &a = m_zeroes_state[i].m1; - point &b = m_zeroes_state[i].m2; - point &c = m_zeroes_state[i].m3; - - real s = abs(b.err) / (abs(a.err) + abs(b.err)); - real newc = b.x + s * (a.x - b.x); - - /* If the third point didn't change since last iteration, - * we may be at an inflection point. Use the midpoint to get - * out of this situation. */ - c.x = newc != c.x ? newc : (a.x + b.x) / 2; - c.err = eval_estimate(c.x) - eval_func(c.x); - - if ((a.err < zero && c.err < zero) - || (a.err > zero && c.err > zero)) - a = c; - else - b = c; - - m_answers.push(i); - } - else if (i < 2000) - { - point &a = m_extrema_state[i - 1000].m1; - point &b = m_extrema_state[i - 1000].m2; - point &c = m_extrema_state[i - 1000].m3; - point d; - - real d1 = c.x - a.x, d2 = c.x - b.x; - real k1 = d1 * (c.err - b.err); - real k2 = d2 * (c.err - a.err); - d.x = c.x - (d1 * k1 - d2 * k2) / (k1 - k2) / 2; - - /* If parabolic interpolation failed, pick a number - * inbetween. */ - if (d.x <= a.x || d.x >= b.x) - d.x = (a.x + b.x) / 2; - - d.err = eval_error(d.x); - - /* Update bracketing depending on the new point. */ - if (d.err < c.err) - { - (d.x > c.x ? b : a) = d; - } - else - { - (d.x > c.x ? a : b) = c; - c = d; - } - - m_answers.push(i); - } - } -} - diff --git a/tools/lolremez/solver.h b/tools/lolremez/solver.h deleted file mode 100644 index f5a84c12..00000000 --- a/tools/lolremez/solver.h +++ /dev/null @@ -1,75 +0,0 @@ -// -// LolRemez - Remez algorithm implementation -// -// Copyright © 2005—2015 Sam Hocevar -// -// This program is free software. It comes without any warranty, to -// the extent permitted by applicable law. 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 the WTFPL Task Force. -// See http://www.wtfpl.net/ for more details. -// - -#pragma once - -// -// The remez_solver class -// ---------------------- -// - -#include - -#include "expression.h" - -class remez_solver -{ -public: - remez_solver(int order, int decimals); - ~remez_solver(); - - void run(lol::real a, lol::real b, - char const *func, char const *weight = nullptr); - -private: - void remez_init(); - void remez_step(); - - void find_zeroes(); - void find_extrema(); - - void worker_thread(); - - void print_poly(); - - lol::real eval_estimate(lol::real const &x); - lol::real eval_func(lol::real const &x); - lol::real eval_weight(lol::real const &x); - lol::real eval_error(lol::real const &x); - -private: - /* User-defined parameters */ - expression m_func, m_weight; - int m_order, m_decimals; - bool m_has_weight; - - /* Solver state */ - lol::polynomial m_estimate; - - lol::array m_zeroes; - lol::array m_control; - - lol::real m_k1, m_k2, m_epsilon, m_error; - - struct point - { - lol::real x, err; - }; - - lol::array m_zeroes_state; - lol::array m_extrema_state; - - /* Threading information */ - lol::array m_workers; - lol::queue m_questions, m_answers; -}; -