diff --git a/tools/lolremez/expression.h b/tools/lolremez/expression.h new file mode 100644 index 00000000..592d8d27 --- /dev/null +++ b/tools/lolremez/expression.h @@ -0,0 +1,315 @@ +// +// 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 + +// +// Parser tools for a simple calculator grammar with + - * / +// + +#include "pegtl.hh" + +namespace grammar +{ + +using namespace pegtl; + +enum class op : uint8_t +{ + /* Variables and constants */ + x, + constant, + /* Unary functions/operators */ + plus, minus, + 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, +}; + +// Map operation enums to pegl::string<> rules +template struct r_call_string {}; + +template<> struct r_call_string : string<'s','q','r','t'> {}; +template<> struct r_call_string : string<'c','b','r','t'> {}; +template<> struct r_call_string : string<'e','x','p'> {}; +template<> struct r_call_string : string<'e','x','p','2'> {}; +template<> struct r_call_string : string<'l','o','g'> {}; +template<> struct r_call_string : string<'l','o','g','2'> {}; +template<> struct r_call_string : string<'l','o','g','1','0'> {}; +template<> struct r_call_string : string<'s','i','n'> {}; +template<> struct r_call_string : string<'c','o','s'> {}; +template<> struct r_call_string : string<'t','a','n'> {}; +template<> struct r_call_string : string<'a','s','i','n'> {}; +template<> struct r_call_string : string<'a','c','o','s'> {}; +template<> struct r_call_string : string<'a','t','a','n'> {}; +template<> struct r_call_string : string<'s','i','n','h'> {}; +template<> struct r_call_string : string<'c','o','s','h'> {}; +template<> struct r_call_string : string<'t','a','n','h'> {}; + +template<> struct r_call_string : string<'a','t','a','n','2'> {}; +template<> struct r_call_string : string<'p','o','w'> {}; +template<> struct r_call_string : string<'m','i','n'> {}; +template<> struct r_call_string : string<'m','a','x'> {}; + +struct expression +{ + /* + * Parse arithmetic expression in x, e.g. 2*x+3 + */ + void parse(std::string const &str) + { + m_ops.Empty(); + m_constants.Empty(); + + basic_parse_string(str, this); + } + + /* + * Evaluate expression at x + */ + lol::real eval(lol::real const &x) + { + lol::array stack; + lol::real tmp; + + for (int i = 0; i < m_ops.Count(); ++i) + { + switch (m_ops[i].m1) + { + case op::plus: + break; + case op::minus: + stack.Push(-stack.Pop()); + break; + case op::sqrt: + stack.Push(sqrt(stack.Pop())); + break; + case op::cbrt: + stack.Push(cbrt(stack.Pop())); + break; + case op::exp: + stack.Push(exp(stack.Pop())); + break; + case op::exp2: + stack.Push(exp2(stack.Pop())); + break; + case op::log: + stack.Push(log(stack.Pop())); + break; + case op::log2: + stack.Push(log2(stack.Pop())); + break; + case op::log10: + stack.Push(log10(stack.Pop())); + break; + case op::sin: + stack.Push(sin(stack.Pop())); + break; + case op::cos: + stack.Push(cos(stack.Pop())); + break; + case op::tan: + stack.Push(tan(stack.Pop())); + break; + case op::asin: + stack.Push(asin(stack.Pop())); + break; + case op::acos: + stack.Push(acos(stack.Pop())); + break; + case op::atan: + stack.Push(atan(stack.Pop())); + break; + case op::sinh: + stack.Push(sinh(stack.Pop())); + break; + case op::cosh: + stack.Push(cosh(stack.Pop())); + break; + case op::tanh: + stack.Push(tanh(stack.Pop())); + break; + + case op::add: + tmp = stack.Pop(); + stack.Push(stack.Pop() + tmp); + break; + case op::sub: + tmp = stack.Pop(); + stack.Push(stack.Pop() - tmp); + break; + case op::mul: + tmp = stack.Pop(); + stack.Push(stack.Pop() * tmp); + break; + case op::div: + tmp = stack.Pop(); + stack.Push(stack.Pop() / tmp); + break; + case op::atan2: + tmp = stack.Pop(); + stack.Push(atan2(stack.Pop(), tmp)); + break; + case op::pow: + tmp = stack.Pop(); + stack.Push(pow(stack.Pop(), tmp)); + break; + case op::min: + tmp = stack.Pop(); + stack.Push(min(stack.Pop(), tmp)); + break; + case op::max: + tmp = stack.Pop(); + stack.Push(max(stack.Pop(), tmp)); + break; + + case op::x: + stack.Push(x); + break; + + case op::constant: + stack.Push(m_constants[m_ops[i].m2]); + break; + } + } + + return stack.Pop(); + } + +private: + lol::array m_ops; + lol::array m_constants; + +private: + struct do_constant : action_base + { + static void apply(std::string const &ctx, expression *that) + { + /* FIXME: check if the constant is already in the list */ + that->m_ops.Push(op::constant, that->m_constants.Count()); + that->m_constants.Push(lol::real(ctx.c_str())); + } + }; + + template + struct do_op : action_base> + { + static void apply(std::string const &ctx, expression *that) + { + UNUSED(ctx); + that->m_ops.Push(OP, -1); + } + }; + + struct do_success : action_base + { + static void apply(std::string const &ctx, expression *that) + { + UNUSED(ctx, that); + } + }; + + struct r_expr; + + // Rule for unary function calls + template + struct r_call_unary : ifapply, + one<'('>, + r_expr, + one<')'>>, + do_op> {}; + + // Rule for binary function calls + template + struct r_call_binary : ifapply, + one<'('>, + r_expr, + one<','>, + r_expr, + one<')'>>, + do_op> {}; + + // r_constant <- + + struct r_constant : ifapply, do_constant> {}; + + // r_var <- "x" + struct r_var : ifapply, do_op> {}; + + // r_call <- "(" r_expr ")" + // / "(" r_expr "," r_expr ")" + struct r_call : sor, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_unary, + r_call_binary, + r_call_binary, + r_call_binary, + r_call_binary> {}; + + // r_parentheses <- "(" r_expr ")" + struct r_parentheses : seq, r_expr, one<')'>> {}; + + // r_terminal <- r_call / r_var / r_constant / r_parentheses + struct r_terminal : sor {}; + + // r_exponent <- ( "^" | "**" ) r_terminal + // r_factor <- r_terminal ( r_exponent ) * + struct r_exponent : ifapply, string<'*', '*'>>, + r_terminal>, do_op> {}; + struct r_factor : seq> {}; + + // r_mul <- "*" r_factor + // r_div <- "/" r_factor + // term <- r_factor ( r_mul / r_div ) * + struct r_mul : ifapply, r_factor>, do_op> {}; + struct r_div : ifapply, r_factor>, do_op> {}; + struct r_term : seq>> {}; + + // add <- "+" r_term + // sub <- "-" r_term + // r_expr <- r_term ( add / sub ) * + struct r_add : ifapply, r_term>, do_op> {}; + struct r_sub : ifapply, r_term>, do_op> {}; + struct r_expr : seq>> {}; + + // stmt <- r_expr + struct r_stmt : ifapply, do_success> {}; +}; + +} /* namespace grammar */ + +using grammar::expression; + diff --git a/tools/lolremez/lolremez.cpp b/tools/lolremez/lolremez.cpp index 536342ac..6e606304 100644 --- a/tools/lolremez/lolremez.cpp +++ b/tools/lolremez/lolremez.cpp @@ -1,11 +1,13 @@ // -// LolRemez - Remez algorithm implementation +// LolRemez - Remez algorithm implementation // -// Copyright: (c) 2005-2013 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://www.wtfpl.net/ for more details. +// 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 @@ -17,26 +19,68 @@ #include #include "solver.h" +#include "expression.h" using lol::real; +using lol::String; -/* See the tutorial at http://lolengine.net/wiki/doc/maths/remez */ -real f(real const &x) +void FAIL() { - return exp(x); -} + printf("Usage:\n"); + printf(" lolremez [-d degree] [-i xmin xmax] x-expression [x-error]\n"); + printf("Example:\n"); + printf(" lolremez -d 4 -i -1 1 \"atan(exp(1+x))\"\n"); + printf(" lolremez -d 4 -i -1 1 \"atan(exp(1+x))\" \"exp(1+x)\"\n"); -real g(real const &x) -{ - return exp(x); + exit(EXIT_FAILURE); } +/* See the tutorial at http://lolengine.net/wiki/doc/maths/remez */ int main(int argc, char **argv) { - UNUSED(argc, argv); + char const *xmin = "-1", *xmax = "1"; + char const *f = nullptr, *g = nullptr; + int degree = 4; + + for (int i = 1; i < argc; ++i) + { + if (argv[i] == String("-d")) + { + if (i + 1 >= argc) + FAIL(); + + degree = atoi(argv[++i]); + } + else if (argv[i] == String("-i")) + { + if (i + 2 >= argc) + FAIL(); + + xmin = argv[++i]; + xmax = argv[++i]; + i += 2; + } + else if (g) + { + FAIL(); + } + else if (f) + { + g = argv[i]; + } + else + { + f = argv[i]; + } + } + + if (!f || real(xmin) >= real(xmax)) + { + FAIL(); + } - RemezSolver solver(4, 40); - solver.Run(-1, 1, f, g); + RemezSolver solver(degree + 1, 20); + solver.Run(xmin, xmax, f, g); return 0; } diff --git a/tools/lolremez/matrix.h b/tools/lolremez/matrix.h index 5f9292a5..9b2ea249 100644 --- a/tools/lolremez/matrix.h +++ b/tools/lolremez/matrix.h @@ -1,11 +1,13 @@ // -// LolRemez - Remez algorithm implementation +// LolRemez - Remez algorithm implementation // -// Copyright: (c) 2010-2013 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://www.wtfpl.net/ for more details. +// 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 diff --git a/tools/lolremez/solver.cpp b/tools/lolremez/solver.cpp index a60d3c59..1c81db01 100644 --- a/tools/lolremez/solver.cpp +++ b/tools/lolremez/solver.cpp @@ -1,11 +1,13 @@ // -// LolRemez - Remez algorithm implementation +// LolRemez - Remez algorithm implementation // -// Copyright: (c) 2005-2013 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://www.wtfpl.net/ for more details. +// 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 @@ -26,18 +28,23 @@ using lol::real; RemezSolver::RemezSolver(int order, int decimals) : m_order(order), - m_decimals(decimals) + m_decimals(decimals), + m_has_weight(false) { } -void RemezSolver::Run(real a, real b, - RemezSolver::RealFunc *func, - RemezSolver::RealFunc *weight) +void RemezSolver::Run(real a, real b, char const *func, char const *weight) { using std::printf; - m_func = func; - m_weight = 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)); @@ -303,7 +310,7 @@ real RemezSolver::EvalEstimate(real const &x) real RemezSolver::EvalFunc(real const &x) { Timer t; - real ret = m_func(x * m_k2 + m_k1); + real ret = m_func.eval(x * m_k2 + m_k1); m_stats_func += t.Get(); return ret; } @@ -311,7 +318,7 @@ real RemezSolver::EvalFunc(real const &x) real RemezSolver::EvalWeight(real const &x) { Timer t; - real ret = m_weight ? m_weight(x * m_k2 + m_k1) : real(1); + real ret = m_has_weight ? m_weight.eval(x * m_k2 + m_k1) : real(1); m_stats_weight += t.Get(); return ret; } diff --git a/tools/lolremez/solver.h b/tools/lolremez/solver.h index 084a5a0d..1fb0922b 100644 --- a/tools/lolremez/solver.h +++ b/tools/lolremez/solver.h @@ -1,11 +1,13 @@ // -// LolRemez - Remez algorithm implementation +// LolRemez - Remez algorithm implementation // -// Copyright: (c) 2010-2013 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://www.wtfpl.net/ for more details. +// 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 @@ -17,15 +19,15 @@ #include +#include "expression.h" + class RemezSolver { public: - typedef lol::real RealFunc(lol::real const &x); - RemezSolver(int order, int decimals); void Run(lol::real a, lol::real b, - RealFunc *func, RealFunc *weight = nullptr); + char const *func, char const *weight = nullptr); private: void Init(); @@ -40,8 +42,9 @@ private: private: /* User-defined parameters */ - RealFunc *m_func, *m_weight; + expression m_func, m_weight; int m_order, m_decimals; + bool m_has_weight; /* Solver state */ lol::polynomial m_estimate;