lolremez: implement a mathematical expression evaluator.

tools/lolremez/expression.h

@@ -0,0 +1,315 @@
+//
+//  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
+
+//
+// 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<op OP> struct r_call_string {};
+
+template<> struct r_call_string<op::sqrt>  : string<'s','q','r','t'> {};
+template<> struct r_call_string<op::cbrt>  : string<'c','b','r','t'> {};
+template<> struct r_call_string<op::exp>   : string<'e','x','p'> {};
+template<> struct r_call_string<op::exp2>  : string<'e','x','p','2'> {};
+template<> struct r_call_string<op::log>   : string<'l','o','g'> {};
+template<> struct r_call_string<op::log2>  : string<'l','o','g','2'> {};
+template<> struct r_call_string<op::log10> : string<'l','o','g','1','0'> {};
+template<> struct r_call_string<op::sin>   : string<'s','i','n'> {};
+template<> struct r_call_string<op::cos>   : string<'c','o','s'> {};
+template<> struct r_call_string<op::tan>   : string<'t','a','n'> {};
+template<> struct r_call_string<op::asin>  : string<'a','s','i','n'> {};
+template<> struct r_call_string<op::acos>  : string<'a','c','o','s'> {};
+template<> struct r_call_string<op::atan>  : string<'a','t','a','n'> {};
+template<> struct r_call_string<op::sinh>  : string<'s','i','n','h'> {};
+template<> struct r_call_string<op::cosh>  : string<'c','o','s','h'> {};
+template<> struct r_call_string<op::tanh>  : string<'t','a','n','h'> {};
+
+template<> struct r_call_string<op::atan2> : string<'a','t','a','n','2'> {};
+template<> struct r_call_string<op::pow>   : string<'p','o','w'> {};
+template<> struct r_call_string<op::min>   : string<'m','i','n'> {};
+template<> struct r_call_string<op::max>   : 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<r_stmt>(str, this);
+    }
+        
+    /*
+     * Evaluate expression at x
+     */
+    lol::real eval(lol::real const &x)
+    {
+        lol::array<lol::real> 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<op, int> m_ops;
+    lol::array<lol::real> m_constants;
+
+private:
+    struct do_constant : action_base<do_constant>
+    {
+        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<op OP>
+    struct do_op : action_base<do_op<OP>>
+    {
+        static void apply(std::string const &ctx, expression *that)
+        {
+            UNUSED(ctx);
+            that->m_ops.Push(OP, -1);
+        }
+    };
+
+    struct do_success : action_base<do_success>
+    {
+        static void apply(std::string const &ctx, expression *that)
+        {
+            UNUSED(ctx, that);
+        }
+    };
+
+    struct r_expr;
+
+    // Rule for unary function calls
+    template<op OP>
+    struct r_call_unary : ifapply<seq<r_call_string<OP>,
+                                      one<'('>,
+                                      r_expr,
+                                      one<')'>>,
+                                  do_op<OP>> {};
+
+    // Rule for binary function calls
+    template<op OP>
+    struct r_call_binary : ifapply<seq<r_call_string<OP>,
+                                       one<'('>,
+                                       r_expr,
+                                       one<','>,
+                                       r_expr,
+                                       one<')'>>,
+                                   do_op<OP>> {};
+
+    // r_constant <- <digit> +
+    struct r_constant : ifapply<plus<digit>, do_constant> {};
+
+    // r_var <- "x"
+    struct r_var : ifapply<one<'x'>, do_op<op::x>> {};
+
+    // r_call <- <unary_op> "(" r_expr ")"
+    //         / <binary_op> "(" r_expr "," r_expr ")"
+    struct r_call : sor<r_call_unary<op::sqrt>,
+                        r_call_unary<op::cbrt>,
+                        r_call_unary<op::exp>,
+                        r_call_unary<op::exp2>,
+                        r_call_unary<op::log>,
+                        r_call_unary<op::log2>,
+                        r_call_unary<op::log10>,
+                        r_call_unary<op::sin>,
+                        r_call_unary<op::cos>,
+                        r_call_unary<op::tan>,
+                        r_call_unary<op::asin>,
+                        r_call_unary<op::acos>,
+                        r_call_unary<op::atan>,
+                        r_call_unary<op::sinh>,
+                        r_call_unary<op::cosh>,
+                        r_call_unary<op::tanh>,
+                        r_call_binary<op::atan2>,
+                        r_call_binary<op::pow>,
+                        r_call_binary<op::min>,
+                        r_call_binary<op::max>> {};
+
+    // r_parentheses <- "(" r_expr ")"
+    struct r_parentheses : seq<one<'('>, r_expr, one<')'>> {};
+
+    // r_terminal <- r_call / r_var / r_constant / r_parentheses
+    struct r_terminal : sor<r_call,
+                            r_var,
+                            r_constant,
+                            r_parentheses> {};
+
+    // r_exponent <- ( "^" | "**" ) r_terminal
+    // r_factor <- r_terminal ( r_exponent ) *
+    struct r_exponent : ifapply<seq<sor<one<'^'>, string<'*', '*'>>,
+                                    r_terminal>, do_op<op::pow>> {};
+    struct r_factor : seq<r_terminal,
+                          star<r_exponent>> {};
+
+    // r_mul <- "*" r_factor
+    // r_div <- "/" r_factor
+    // term <- r_factor ( r_mul / r_div ) *
+    struct r_mul : ifapply<seq<one<'*'>, r_factor>, do_op<op::mul>> {};
+    struct r_div : ifapply<seq<one<'/'>, r_factor>, do_op<op::div>> {};
+    struct r_term : seq<r_factor,
+                        star<sor<r_mul, r_div>>> {};
+
+    // add <- "+" r_term
+    // sub <- "-" r_term
+    // r_expr <- r_term ( add / sub ) *
+    struct r_add : ifapply<seq<one<'+'>, r_term>, do_op<op::add>> {};
+    struct r_sub : ifapply<seq<one<'-'>, r_term>, do_op<op::sub>> {};
+    struct r_expr : seq<r_term,
+                        star<sor<r_add, r_sub>>> {};
+
+    // stmt <- r_expr <end>
+    struct r_stmt : ifapply<seq<r_expr, eof>, do_success> {};
+};
+
+} /* namespace grammar */
+
+using grammar::expression;
+

tools/lolremez/lolremez.cpp

@@ -1,11 +1,13 @@
 //
-// LolRemez - Remez algorithm implementation
+//  LolRemez - Remez algorithm implementation
 //
-// Copyright: (c) 2005-2013 Sam Hocevar <sam@hocevar.net>
-//   This program is free software; you can redistribute it and/or
-//   modify it under the terms of the Do What The Fuck You Want To
-//   Public License, Version 2, as published by Sam Hocevar. See
-//   http://www.wtfpl.net/ for more details.
+//  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
@@ -17,26 +19,68 @@
 #include <lol/math/real.h>
 
 #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;
 }

tools/lolremez/matrix.h

@@ -1,11 +1,13 @@
 //
-// LolRemez - Remez algorithm implementation
+//  LolRemez - Remez algorithm implementation
 //
-// Copyright: (c) 2010-2013 Sam Hocevar <sam@hocevar.net>
-//   This program is free software; you can redistribute it and/or
-//   modify it under the terms of the Do What The Fuck You Want To
-//   Public License, Version 2, as published by Sam Hocevar. See
-//   http://www.wtfpl.net/ for more details.
+//  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

tools/lolremez/solver.cpp

@@ -1,11 +1,13 @@
 //
-// LolRemez - Remez algorithm implementation
+//  LolRemez - Remez algorithm implementation
 //
-// Copyright: (c) 2005-2013 Sam Hocevar <sam@hocevar.net>
-//   This program is free software; you can redistribute it and/or
-//   modify it under the terms of the Do What The Fuck You Want To
-//   Public License, Version 2, as published by Sam Hocevar. See
-//   http://www.wtfpl.net/ for more details.
+//  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
@@ -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;
 }

tools/lolremez/solver.h

@@ -1,11 +1,13 @@
 //
-// LolRemez - Remez algorithm implementation
+//  LolRemez - Remez algorithm implementation
 //
-// Copyright: (c) 2010-2013 Sam Hocevar <sam@hocevar.net>
-//   This program is free software; you can redistribute it and/or
-//   modify it under the terms of the Do What The Fuck You Want To
-//   Public License, Version 2, as published by Sam Hocevar. See
-//   http://www.wtfpl.net/ for more details.
+//  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
@@ -17,15 +19,15 @@
 
 #include <cstdio>
 
+#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<lol::real> m_estimate;