It will now be available at https://github.com/samhocevar/lolremezlegacy
@@ -58,10 +58,6 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Test", "Test", "{E4DFEBF9-C | |||||
EndProject | EndProject | ||||
Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "13_shader_builder", "..\doc\tutorial\13_shader_builder.vcxproj", "{F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}" | Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "13_shader_builder", "..\doc\tutorial\13_shader_builder.vcxproj", "{F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}" | ||||
EndProject | 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}" | Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Samples", "Samples", "{B6297FF2-63D0-41EE-BE13-EFF720C9B0FA}" | ||||
EndProject | EndProject | ||||
Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "15_lolimgui", "..\doc\tutorial\15_lolimgui.vcxproj", "{81C83B42-D00A-4FA3-9A3D-80F9D46524BF}" | 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|Win32.Build.0 = Release|Win32 | ||||
{F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}.Release|x64.ActiveCfg = Release|x64 | {F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}.Release|x64.ActiveCfg = Release|x64 | ||||
{F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4}.Release|x64.Build.0 = 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.ActiveCfg = Debug|ORBIS | ||||
{81C83B42-D00A-4FA3-9A3D-80F9D46524BF}.Debug|ORBIS.Build.0 = Debug|ORBIS | {81C83B42-D00A-4FA3-9A3D-80F9D46524BF}.Debug|ORBIS.Build.0 = Debug|ORBIS | ||||
{81C83B42-D00A-4FA3-9A3D-80F9D46524BF}.Debug|Win32.ActiveCfg = Debug|Win32 | {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} | {E05E23A5-67DE-42B5-98A3-E63CCE0CC0AF} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} | ||||
{E4DFEBF9-C310-462F-9876-7EB59C1E4D4E} = {1AFD580B-98B8-4689-B661-38C41132C60E} | {E4DFEBF9-C310-462F-9876-7EB59C1E4D4E} = {1AFD580B-98B8-4689-B661-38C41132C60E} | ||||
{F59FA82C-DDB9-4EE2-80AE-CB0E4C6567A4} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} | {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} | {B6297FF2-63D0-41EE-BE13-EFF720C9B0FA} = {1AFD580B-98B8-4689-B661-38C41132C60E} | ||||
{81C83B42-D00A-4FA3-9A3D-80F9D46524BF} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} | {81C83B42-D00A-4FA3-9A3D-80F9D46524BF} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} | ||||
{31B96262-1C41-43B9-BA38-27AA385B05DB} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} | {31B96262-1C41-43B9-BA38-27AA385B05DB} = {E74CF679-CA2A-47E9-B1F4-3779D6AC6B04} | ||||
@@ -21,8 +21,6 @@ AM_DEFAULT_VERBOSITY=0 | |||||
dnl Versioning of the separate software we ship | dnl Versioning of the separate software we ship | ||||
LOLUNIT_VERSION=0.1 | LOLUNIT_VERSION=0.1 | ||||
AC_SUBST(LOLUNIT_VERSION) | AC_SUBST(LOLUNIT_VERSION) | ||||
LOLREMEZ_VERSION=0.2 | |||||
AC_SUBST(LOLREMEZ_VERSION) | |||||
AC_SUBST(lol_srcdir, '${top_srcdir}') | AC_SUBST(lol_srcdir, '${top_srcdir}') | ||||
AC_SUBST(lol_builddir, '${top_builddir}') | AC_SUBST(lol_builddir, '${top_builddir}') | ||||
@@ -262,7 +260,6 @@ AC_CONFIG_FILES( | |||||
doc/samples/sandbox/Makefile | doc/samples/sandbox/Makefile | ||||
doc/tutorial/Makefile | doc/tutorial/Makefile | ||||
tools/Makefile | tools/Makefile | ||||
tools/lolremez/Makefile | |||||
tools/lolunit/Makefile | tools/lolunit/Makefile | ||||
tools/vimlol/Makefile | tools/vimlol/Makefile | ||||
tools/vslol/Makefile | tools/vslol/Makefile | ||||
@@ -2,7 +2,6 @@ | |||||
include $(top_srcdir)/build/autotools/common.am | include $(top_srcdir)/build/autotools/common.am | ||||
SUBDIRS = | SUBDIRS = | ||||
SUBDIRS += lolremez | |||||
SUBDIRS += lolunit | SUBDIRS += lolunit | ||||
SUBDIRS += vimlol | SUBDIRS += vimlol | ||||
SUBDIRS += vslol | SUBDIRS += vslol | ||||
@@ -1,2 +0,0 @@ | |||||
# Our binaries | |||||
lolremez |
@@ -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@ | |||||
@@ -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 | |||||
@@ -1,349 +0,0 @@ | |||||
// | |||||
// LolRemez — Remez algorithm implementation | |||||
// | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | |||||
// 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<lol::real> 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<id, int> m_ops; | |||||
lol::array<lol::real> m_constants; | |||||
private: | |||||
struct r_expr; | |||||
// r_ <- <blank> * | |||||
struct _ : star<space> {}; | |||||
// r_constant <- <digit> + ( "." <digit> * ) ? ( [eE] [+-] ? <digit> + ) ? | |||||
struct r_constant : seq<plus<digit>, | |||||
opt<seq<one<'.'>, | |||||
star<digit>>>, | |||||
opt<seq<one<'e', 'E'>, | |||||
opt<one<'+', '-'>>, | |||||
plus<digit>>>> {}; | |||||
// r_var <- "x" | |||||
struct r_var : pegtl_string_t("x") {}; | |||||
// r_binary_call <- <r_binary_fun> "(" r_expr "," r_expr ")" | |||||
struct r_binary_fun : sor<pegtl_string_t("atan2"), | |||||
pegtl_string_t("pow"), | |||||
pegtl_string_t("min"), | |||||
pegtl_string_t("max")> {}; | |||||
struct r_binary_call : seq<r_binary_fun, | |||||
_, one<'('>, | |||||
_, r_expr, | |||||
_, one<','>, | |||||
_, r_expr, | |||||
_, one<')'>> {}; | |||||
// r_unary_call <- <r_unary_fun> "(" r_expr ")" | |||||
struct r_unary_fun : sor<pegtl_string_t("abs"), | |||||
pegtl_string_t("sqrt"), | |||||
pegtl_string_t("cbrt"), | |||||
pegtl_string_t("exp"), | |||||
pegtl_string_t("exp2"), | |||||
pegtl_string_t("log"), | |||||
pegtl_string_t("log2"), | |||||
pegtl_string_t("log10"), | |||||
pegtl_string_t("sin"), | |||||
pegtl_string_t("cos"), | |||||
pegtl_string_t("tan"), | |||||
pegtl_string_t("asin"), | |||||
pegtl_string_t("acos"), | |||||
pegtl_string_t("atan"), | |||||
pegtl_string_t("sinh"), | |||||
pegtl_string_t("cosh"), | |||||
pegtl_string_t("tanh")> {}; | |||||
struct r_unary_call : seq<r_unary_fun, | |||||
_, one<'('>, | |||||
_, r_expr, | |||||
_, one<')'>> {}; | |||||
// r_call <- r_binary_call / r_unary_call | |||||
struct r_call : sor<r_binary_call, | |||||
r_unary_call> {}; | |||||
// r_parentheses <- "(" r_expr ")" | |||||
struct r_parentheses : seq<one<'('>, | |||||
pad<r_expr, space>, | |||||
one<')'>> {}; | |||||
// r_terminal <- r_call / r_var / r_constant / r_parentheses | |||||
struct r_terminal : sor<r_call, | |||||
r_var, | |||||
r_constant, | |||||
r_parentheses> {}; | |||||
// r_signed <- "-" r_signed / "+" r_signed / r_terminal | |||||
struct r_signed; | |||||
struct r_minus : seq<one<'-'>, _, r_signed> {}; | |||||
struct r_signed : sor<r_minus, | |||||
seq<one<'+'>, _, r_signed>, | |||||
r_terminal> {}; | |||||
// r_exponent <- ( "^" / "**" ) r_signed | |||||
struct r_exponent : seq<pad<sor<one<'^'>, | |||||
string<'*', '*'>>, space>, | |||||
r_signed> {}; | |||||
// r_factor <- r_signed ( r_exponent ) * | |||||
struct r_factor : seq<r_signed, | |||||
star<r_exponent>> {}; | |||||
// 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_factor, | |||||
star<sor<r_mul, r_div>>> {}; | |||||
// 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_term, | |||||
star<sor<r_add, r_sub>>> {}; | |||||
// r_stmt <- r_expr <end> | |||||
struct r_stmt : must<pad<r_expr, space>, pegtl::eof> {}; | |||||
// | |||||
// Default actions | |||||
// | |||||
template<typename R> | |||||
struct action : nothing<R> {}; | |||||
template<id OP> | |||||
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<r_stmt, action>(str, "expression", this); | |||||
} | |||||
}; | |||||
// | |||||
// Rule specialisations for simple operators | |||||
// | |||||
template<> struct expression::action<expression::r_var> : generic_action<id::x> {}; | |||||
template<> struct expression::action<expression::r_exponent> : generic_action<id::pow> {}; | |||||
template<> struct expression::action<expression::r_mul> : generic_action<id::mul> {}; | |||||
template<> struct expression::action<expression::r_div> : generic_action<id::div> {}; | |||||
template<> struct expression::action<expression::r_add> : generic_action<id::add> {}; | |||||
template<> struct expression::action<expression::r_sub> : generic_action<id::sub> {}; | |||||
template<> struct expression::action<expression::r_minus> : generic_action<id::minus> {}; | |||||
// | |||||
// Rule specialisations for unary and binary function calls | |||||
// | |||||
template<> | |||||
struct expression::action<expression::r_binary_call> | |||||
{ | |||||
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<expression::r_unary_call> | |||||
{ | |||||
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<expression::r_constant> | |||||
{ | |||||
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; | |||||
@@ -1,129 +0,0 @@ | |||||
// | |||||
// LolRemez - Remez algorithm implementation | |||||
// | |||||
// Copyright © 2005—2016 Sam Hocevar <sam@hocevar.net> | |||||
// | |||||
// 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 <lol/engine.h> | |||||
#include <lol/math/real.h> | |||||
#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 <sam@hocevar.net>\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 <sam@hocevar.net>.\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> degree of final polynomial\n"); | |||||
printf(" -r, --range <xmin>:<xmax> 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 <sam@hocevar.net>.\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<String> 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; | |||||
} | |||||
@@ -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 |
@@ -1,73 +0,0 @@ | |||||
<?xml version="1.0" encoding="utf-8"?> | |||||
<Project DefaultTargets="Build" ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003"> | |||||
<PropertyGroup Label="LolMacros"> | |||||
<LolDir Condition="Exists('$(SolutionDir)\lol')">$(SolutionDir)\lol</LolDir> | |||||
<LolDir Condition="!Exists('$(SolutionDir)\lol')">$(SolutionDir)\..</LolDir> | |||||
</PropertyGroup> | |||||
<ItemGroup Label="ProjectConfigurations"> | |||||
<ProjectConfiguration Include="Debug|ORBIS"> | |||||
<Configuration>Debug</Configuration> | |||||
<Platform>ORBIS</Platform> | |||||
</ProjectConfiguration> | |||||
<ProjectConfiguration Include="Debug|Win32"> | |||||
<Configuration>Debug</Configuration> | |||||
<Platform>Win32</Platform> | |||||
</ProjectConfiguration> | |||||
<ProjectConfiguration Include="Debug|x64"> | |||||
<Configuration>Debug</Configuration> | |||||
<Platform>x64</Platform> | |||||
</ProjectConfiguration> | |||||
<ProjectConfiguration Include="Release|ORBIS"> | |||||
<Configuration>Release</Configuration> | |||||
<Platform>ORBIS</Platform> | |||||
</ProjectConfiguration> | |||||
<ProjectConfiguration Include="Release|Win32"> | |||||
<Configuration>Release</Configuration> | |||||
<Platform>Win32</Platform> | |||||
</ProjectConfiguration> | |||||
<ProjectConfiguration Include="Release|x64"> | |||||
<Configuration>Release</Configuration> | |||||
<Platform>x64</Platform> | |||||
</ProjectConfiguration> | |||||
</ItemGroup> | |||||
<ItemGroup> | |||||
<ClInclude Include="expression.h" /> | |||||
<ClInclude Include="matrix.h" /> | |||||
<ClInclude Include="solver.h" /> | |||||
</ItemGroup> | |||||
<ItemGroup> | |||||
<ClCompile Include="lolremez.cpp" /> | |||||
<ClCompile Include="solver.cpp" /> | |||||
</ItemGroup> | |||||
<ItemGroup> | |||||
<ProjectReference Include="$(LolDir)\src\lol-core.vcxproj"> | |||||
<Project>{9e62f2fe-3408-4eae-8238-fd84238ceeda}</Project> | |||||
</ProjectReference> | |||||
<ProjectReference Include="$(LolDir)\src\3rdparty\lol-bullet.vcxproj"> | |||||
<Project>{83d3b207-c601-4025-8f41-01dedc354661}</Project> | |||||
</ProjectReference> | |||||
<ProjectReference Include="$(LolDir)\src\3rdparty\lol-lua.vcxproj"> | |||||
<Project>{d84021ca-b233-4e0f-8a52-071b83bbccc4}</Project> | |||||
</ProjectReference> | |||||
</ItemGroup> | |||||
<PropertyGroup Label="Globals"> | |||||
<ProjectGuid>{73F1A804-1116-46C3-922A-9C0ADEB33F52}</ProjectGuid> | |||||
<ConfigurationType>Application</ConfigurationType> | |||||
<Keyword>Win32Proj</Keyword> | |||||
</PropertyGroup> | |||||
<Import Project="$(LolDir)\build\msbuild\lol.config.props" /> | |||||
<ImportGroup Label="ExtensionSettings"> | |||||
<Import Project="$(LolDir)\build\msbuild\lolfx.props" /> | |||||
</ImportGroup> | |||||
<ImportGroup Label="PropertySheets"> | |||||
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" /> | |||||
<Import Project="$(LolDir)\build\msbuild\lol.vars.props" /> | |||||
</ImportGroup> | |||||
<PropertyGroup Label="UserMacros" /> | |||||
<Import Project="$(LolDir)\build\msbuild\lol.rules.props" /> | |||||
<ItemDefinitionGroup /> | |||||
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" /> | |||||
<ImportGroup Label="ExtensionTargets"> | |||||
<Import Project="$(LolDir)\build\msbuild\lolfx.targets" /> | |||||
</ImportGroup> | |||||
</Project> |
@@ -1,12 +0,0 @@ | |||||
<?xml version="1.0" encoding="utf-8"?> | |||||
<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003"> | |||||
<ItemGroup> | |||||
<ClCompile Include="lolremez.cpp" /> | |||||
<ClCompile Include="solver.cpp" /> | |||||
</ItemGroup> | |||||
<ItemGroup> | |||||
<ClInclude Include="expression.h" /> | |||||
<ClInclude Include="matrix.h" /> | |||||
<ClInclude Include="solver.h" /> | |||||
</ItemGroup> | |||||
</Project> |
@@ -1,97 +0,0 @@ | |||||
// | |||||
// LolRemez - Remez algorithm implementation | |||||
// | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | |||||
// 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<typename T> | |||||
struct linear_system : public array2d<T> | |||||
{ | |||||
inline linear_system<T>(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<T> 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; | |||||
} | |||||
}; | |||||
@@ -1,419 +0,0 @@ | |||||
// | |||||
// LolRemez - Remez algorithm implementation | |||||
// | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | |||||
// 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 <functional> | |||||
#include <lol/engine.h> | |||||
#include <lol/math/real.h> | |||||
#include <lol/math/polynomial.h> | |||||
#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<real> 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<real> system(m_order + 1); | |||||
for (int n = 0; n < m_order + 1; n++) | |||||
{ | |||||
auto p = polynomial<real>::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<real>(); | |||||
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<real>::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<real> 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<real> system(m_order + 2); | |||||
for (int n = 0; n < m_order + 1; n++) | |||||
{ | |||||
auto p = polynomial<real>::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<real>(); | |||||
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<real>::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<real> q ({ -m_k1 / m_k2, real(1) / m_k2 }); | |||||
polynomial<real> 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); | |||||
} | |||||
} | |||||
} | |||||
@@ -1,75 +0,0 @@ | |||||
// | |||||
// LolRemez - Remez algorithm implementation | |||||
// | |||||
// Copyright © 2005—2015 Sam Hocevar <sam@hocevar.net> | |||||
// | |||||
// 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 <cstdio> | |||||
#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<lol::real> m_estimate; | |||||
lol::array<lol::real> m_zeroes; | |||||
lol::array<lol::real> m_control; | |||||
lol::real m_k1, m_k2, m_epsilon, m_error; | |||||
struct point | |||||
{ | |||||
lol::real x, err; | |||||
}; | |||||
lol::array<point, point, point> m_zeroes_state; | |||||
lol::array<point, point, point> m_extrema_state; | |||||
/* Threading information */ | |||||
lol::array<lol::thread *> m_workers; | |||||
lol::queue<int> m_questions, m_answers; | |||||
}; | |||||