same

choosegenerators.cpp

@@ -296,11 +296,12 @@ void ChooseGenerators::on_compile_clicked()
 {
     QString formula = this->ui->formula->text();
     mnd::IterationFormula itf{ mnd::parse(formula.toStdString()) };
+    itf.optimize();
 
     std::string expr = mnd::toString(*itf.expr);
     printf("%s\n", expr.c_str()); fflush(stdout);
-    //chosenGenerator = std::make_unique<mnd::NaiveGenerator>(std::move(itf), mnd::getPrecision<double>());
-    //return;
+    chosenGenerator = std::make_unique<mnd::NaiveGenerator>(std::move(itf), mnd::getPrecision<double>());
+    return;
     mnd::ir::Formula irform = mnd::expand(itf);
     printf("%s\n", irform.toString().c_str()); fflush(stdout);
     auto cg = std::make_unique<mnd::CompiledGenerator>(mnd::compile(irform));

libmandel/include/IterationFormula.h

@@ -6,13 +6,15 @@
 #include <string>
 #include <stdexcept>
 
+#include "Types.h"
+
 namespace mnd
 {
     struct IterationFormula;
 
     struct Constant;
     struct Variable;
-    struct UnaryOperation;
+    struct Negation;
     struct BinaryOperation;
     struct Addition;
     struct Multiplication;
@@ -22,7 +24,7 @@ namespace mnd
     using Expression = std::variant<
             Constant,
             Variable,
-            UnaryOperation,
+            Negation,
             Addition,
             Multiplication,
             Division,
@@ -47,15 +49,22 @@ struct mnd::IterationFormula
     std::unique_ptr<Expression> expr;
     IterationFormula(Expression expr);
 
-    void constantPropagation
+    void optimize(void);
 };
 
 
 struct mnd::Constant
 {
-    double value;
-    inline Constant(double value) :
-        value{ value }
+    mnd::Real re;
+    mnd::Real im;
+    inline Constant(const mnd::Real& value) :
+        re{ value },
+        im{ 0.0 }
+    {}
+
+    inline Constant(const mnd::Real& re, const mnd::Real& im) :
+        re{ re },
+        im{ im }
     {}
 };
 
@@ -66,7 +75,7 @@ struct mnd::Variable
 };
 
 
-struct mnd::UnaryOperation
+struct mnd::Negation
 {
     std::unique_ptr<Expression> operand;
     /*inline UnaryOperation(const UnaryOperation& other) :
@@ -104,6 +113,8 @@ struct mnd::Division : mnd::BinaryOperation
 
 struct mnd::Pow : mnd::BinaryOperation 
 {
+    bool realExponent;
+    bool integerExponent;
 };
 
 

libmandel/include/IterationIR.h

@@ -24,6 +24,9 @@ namespace mnd
         struct Addition;
         struct Subtraction;
         struct Multiplication;
+        struct Division;
+//        struct CPow;
+//        struct RPow;
         struct Atan2;
         struct Pow;
         struct Cos;
@@ -38,6 +41,9 @@ namespace mnd
             Addition,
             Subtraction,
             Multiplication,
+            Division,
+//            CPow,
+//            RPow,
             Atan2,
             Pow,
             Cos,
@@ -71,7 +77,7 @@ struct mnd::ir::NodeBase
 struct mnd::ir::Constant : NodeBase
 {
     mnd::Real value;
-    inline Constant(double val) : value{ val } {}
+    inline Constant(const mnd::Real& val) : value{ val } {}
 };
 
 
@@ -125,6 +131,36 @@ struct mnd::ir::Multiplication : mnd::ir::BinaryOperation
 };
 
 
+struct mnd::ir::Division : mnd::ir::BinaryOperation
+{
+    inline Division(Node* left, Node* right) :
+        BinaryOperation{ left, right } {}
+};
+
+/*
+struct mnd::ir::CPow : mnd::ir::NodeBase
+{
+    Node* re;
+    Node* im;
+    Node* ere;
+    Node* eim;
+    inline CPow(Node* re, Node* im, Node* ere, Node* eim) :
+        re{ re }, im{ im }, ere{ ere }, eim{ eim }
+    {}
+};
+
+
+struct mnd::ir::RPow : mnd::ir::NodeBase
+{
+    Node* re;
+    Node* im;
+    Node* exponent;
+    inline RPow(Node* re, Node* im, Node* exponent) :
+        re{ re }, im{ im }, exponent{ exponent }
+    {}
+};*/
+
+
 struct mnd::ir::Atan2 : mnd::ir::BinaryOperation
 {
     inline Atan2(Node* left, Node* right) :

libmandel/src/IterationCompiler.cpp

@@ -116,11 +116,19 @@ namespace mnd
 
         Reg operator()(const ir::Multiplication& add) {
             auto res = cc.newXmmSd();
+            cc.comment("multiply");
             cc.movapd(res, visitNode(*add.left));
             cc.mulsd(res, visitNode(*add.right));
             return res;
         }
 
+        Reg operator()(const ir::Division& add) {
+            auto res = cc.newXmmSd();
+            cc.movapd(res, visitNode(*add.left));
+            cc.divsd(res, visitNode(*add.right));
+            return res;
+        }
+
         Reg operator()(const ir::Atan2& at2) {
             using namespace asmjit;
             auto y = visitNode(*at2.left);
@@ -286,10 +294,13 @@ namespace mnd
             }
             else {
                 std::string value = std::visit(*this, node);
-                std::string varname = createVarname();
-                code << "float " << varname << " = " << value << ";" << std::endl;
-                nodeData = varname;
-                return varname;
+                if (!std::get_if<ir::Variable>(&node)) {
+                    std::string varname = createVarname();
+                    code << "float " << varname << " = " << value << ";" << std::endl;
+                    nodeData = varname;
+                    return varname;
+                }
+                return value;
             }
         }
 
@@ -317,6 +328,10 @@ namespace mnd
             return "("s + visitNode(*a.left) + ") * (" + visitNode(*a.right) + ")";
         }
 
+        std::string operator()(const ir::Division& a) {
+            return "("s + visitNode(*a.left) + ") / (" + visitNode(*a.right) + ")";
+        }
+
         std::string operator()(const ir::Atan2& a) {
             return "atan2("s + visitNode(*a.left) + ", " + visitNode(*a.right) + ")";
         }
@@ -347,6 +362,7 @@ namespace mnd
         OpenClVisitor ocv;
         std::string newA = ocv.visitNode(*formula.newA);
         std::string newB = ocv.visitNode(*formula.newB);
+
         std::string prelude = 
 "__kernel void iterate(__global float* A, const int width, float xl, float yt, float pixelScaleX, float pixelScaleY, int max, int smooth, int julia, float juliaX, float juliaY) {\n"
 "   int index = get_global_id(0);\n"
@@ -386,7 +402,7 @@ namespace mnd
         code += "b = " + newB + ";\n";
         code += after;
         //code = mnd::getFloat_cl();
-        printf("cl: %s\n", code.c_str());
+        printf("cl: %s\n", code.c_str()); fflush(stdout);
         return code;
     }
 

libmandel/src/IterationFormula.cpp

@@ -4,6 +4,7 @@
 #include <vector>
 #include <stack>
 #include <regex>
+#include <optional>
 
 using mnd::ParseError;
 
@@ -14,6 +15,122 @@ mnd::IterationFormula::IterationFormula(mnd::Expression expr) :
 }
 
 
+struct SimpleOptimizer
+{
+    using Ret = std::optional<mnd::Expression>;
+    void visitExpr(std::unique_ptr<mnd::Expression>& expr)
+    {
+        Ret replacement = std::visit(*this, *expr);
+        if (replacement.has_value()) {
+            expr = std::make_unique<mnd::Expression>(std::move(replacement.value()));
+        }
+    }
+
+    Ret operator() (mnd::Constant& c)
+    {
+        return std::nullopt;
+    }
+
+    Ret operator() (mnd::Variable& v)
+    {
+        if (v.name == "i") {
+            return mnd::Constant{ 0.0, 1.0 };
+        }
+        else {
+            return std::nullopt;
+        }
+    }
+
+    Ret operator() (mnd::Negation& n)
+    {
+        visitExpr(n.operand);
+        auto* valConst = std::get_if<mnd::Constant>(n.operand.get());
+
+        if (valConst) {
+            return mnd::Constant {
+                -valConst->re,
+                -valConst->im
+            };
+        }
+        return std::nullopt;
+    }
+
+    Ret operator() (mnd::Addition& a)
+    {
+        visitExpr(a.left);
+        visitExpr(a.right);
+        auto* leftConst = std::get_if<mnd::Constant>(a.left.get());
+        auto* rightConst = std::get_if<mnd::Constant>(a.right.get());
+
+        if (leftConst && rightConst) {
+            if (a.subtraction) {
+                return mnd::Constant {
+                    leftConst->re - rightConst->re,
+                    leftConst->im - rightConst->im
+                };
+            }
+            else {
+                return mnd::Constant{
+                    leftConst->re + rightConst->re,
+                    leftConst->im + rightConst->im
+                };
+            }
+        }
+        return std::nullopt;
+    }
+
+    Ret operator() (mnd::Multiplication& a)
+    {
+        visitExpr(a.left);
+        visitExpr(a.right);
+        auto* leftConst = std::get_if<mnd::Constant>(a.left.get());
+        auto* rightConst = std::get_if<mnd::Constant>(a.right.get());
+
+        if (leftConst && rightConst) {
+            return mnd::Constant {
+                leftConst->re * rightConst->re - leftConst->im * rightConst->im,
+                (leftConst->re * rightConst->im + leftConst->im * rightConst->re) * 2
+            };
+        }
+        return std::nullopt;
+    }
+
+    Ret operator() (mnd::Division& a)
+    {
+        visitExpr(a.left);
+        visitExpr(a.right);
+        // TODO implement
+        return std::nullopt;
+    }
+
+    Ret operator() (mnd::Pow& a)
+    {
+        visitExpr(a.left);
+        visitExpr(a.right);
+        auto* leftConst = std::get_if<mnd::Constant>(a.left.get());
+        auto* rightConst = std::get_if<mnd::Constant>(a.right.get());
+
+        if (rightConst) {
+            if (rightConst->im == 0) {
+                a.realExponent = true;
+                if (int(rightConst->re) == rightConst->re) {
+                    a.integerExponent = true;
+                }
+            }
+        }
+
+        return std::nullopt;
+    }
+};
+
+
+void mnd::IterationFormula::optimize(void)
+{
+    SimpleOptimizer so;
+    so.visitExpr(this->expr);
+}
+
+
 static const std::string regexIdent = "[A-Za-z][A-Za-z0-9]*";
 static const std::string regexNum = "[1-9][0-9]*";
 static const std::string regexFloat = "(\\d*\\.?\\d+|\\d+\\.?\\d*)([eE][-+]\\d+)?";
@@ -183,13 +300,13 @@ namespace mnd
             std::stringstream ss;
             using T = std::decay_t<decltype(ex)>;
             if constexpr (std::is_same<T, mnd::Constant>::value) {
-                ss << "const[" << ex.value << "]";
+                ss << "const[" << ex.re << "+" << ex.im << "i" << "]";
                 return ss.str();
             }
             else if constexpr (std::is_same<T, mnd::Variable>::value) {
                 return std::string("var[") + ex.name + "]";
             }
-            else if constexpr (std::is_same<T, mnd::UnaryOperation>::value) {
+            else if constexpr (std::is_same<T, mnd::Negation>::value) {
                 return std::string("-(") + toString(*ex.operand) + ")";
             }
             else if constexpr (std::is_same<T, mnd::Addition>::value) {

libmandel/src/IterationGenerator.cpp

@@ -20,6 +20,7 @@ NaiveGenerator::NaiveGenerator(IterationFormula itf,
                                    const mnd::Real& prec) :
     IterationGenerator{ std::move(itf), prec }
 {
+    this->itf.optimize();
 }
 
 
@@ -60,7 +61,7 @@ void NaiveGenerator::generate(const mnd::MandelInfo& info, float* data)
                 if (std::abs(z) >= 4)
                     break;
             }
-            data[i + j * info.bWidth] = k;
+            data[i + j * info.bWidth] = float(k);
             /*if (info.smooth) {
                 if (k >= info.maxIter)
                     data[i + j * info.bWidth] = float(info.maxIter);
@@ -90,7 +91,7 @@ std::complex<double> NaiveGenerator::calc(mnd::Expression& expr, std::complex<do
     std::visit([this, &result, z, c] (auto&& ex) {
         using T = std::decay_t<decltype(ex)>;
         if constexpr (std::is_same<T, mnd::Constant>::value) {
-            result = ex.value;
+            result = std::complex{ ex.re, ex.im };
         }
         else if constexpr (std::is_same<T, mnd::Variable>::value) {
             if (ex.name == "z")
@@ -100,7 +101,7 @@ std::complex<double> NaiveGenerator::calc(mnd::Expression& expr, std::complex<do
             else if (ex.name == "i")
                 result = std::complex{ 0.0, 1.0 };
         }
-        else if constexpr (std::is_same<T, mnd::UnaryOperation>::value) {
+        else if constexpr (std::is_same<T, mnd::Negation>::value) {
             result = -calc(*ex.operand, z, c);
         }
         else if constexpr (std::is_same<T, mnd::Addition>::value) {

libmandel/src/IterationIR.cpp

@@ -14,15 +14,22 @@ namespace mnd
         using NodePair = std::pair<Node*, Node*>;
         util::Arena<Node>& arena;
 
+        Node* zero;
+        Node* half;
+        Node* one;
+
         ConvertVisitor(util::Arena<Node>& arena) :
             arena{ arena }
         {
+            zero = arena.allocate(ir::Constant{ 0.0 });
+            half = arena.allocate(ir::Constant{ 0.5 });
+            one = arena.allocate(ir::Constant{ 1.0 });
         }
 
         NodePair operator() (const Constant& c)
         {
-            Node* cnst = arena.allocate(ir::Constant{ c.value });
-            Node* zero = arena.allocate(ir::Constant{ 0.0 });
+            Node* cnst = arena.allocate(ir::Constant{ c.re });
+            Node* zero = arena.allocate(ir::Constant{ c.im });
 
             return { cnst, zero };
         }
@@ -43,16 +50,13 @@ namespace mnd
                 return { x, y };
             }
             else if (v.name == "i") {
-                Node* x = arena.allocate(ir::Constant{ 0.0 });
-                Node* y = arena.allocate(ir::Constant{ 1.0 });
-
-                return { x, y };
+                return { zero, one };
             }
             else
                 throw "unknown variable";
         }
 
-        NodePair operator() (const UnaryOperation& v)
+        NodePair operator() (const Negation& v)
         {
             auto [opa, opb] = std::visit(*this, *v.operand);
 
@@ -86,6 +90,11 @@ namespace mnd
             auto [a, b] = std::visit(*this, *mul.left);
             auto [c, d] = std::visit(*this, *mul.right);
 
+            return multiplication(a, b, c, d);
+        }
+
+        NodePair multiplication(Node* a, Node* b, Node* c, Node* d)
+        {
             Node* ac = arena.allocate(ir::Multiplication{ a, c });
             Node* bd = arena.allocate(ir::Multiplication{ b, d });
             Node* ad = arena.allocate(ir::Multiplication{ a, d });
@@ -97,6 +106,18 @@ namespace mnd
             return { newa, newb };
         }
 
+        NodePair sq(Node* a, Node* b)
+        {
+            Node* aa = arena.allocate(ir::Multiplication{ a, a });
+            Node* bb = arena.allocate(ir::Multiplication{ b, b });
+            Node* ab = arena.allocate(ir::Multiplication{ a, b });
+
+            Node* newa = arena.allocate(ir::Subtraction{ aa, bb });
+            Node* newb = arena.allocate(ir::Addition{ ab, ab });
+
+            return { newa, newb };
+        }
+
         NodePair operator() (const Division& mul)
         {
             // TODO implement
@@ -109,7 +130,14 @@ namespace mnd
             auto [a, b] = std::visit(*this, *p.left);
             auto [c, d] = std::visit(*this, *p.right);
 
-            auto half = arena.allocate(ir::Constant{ 0.5 });
+            if (p.integerExponent) {
+                if (auto* ex = std::get_if<ir::Constant>(c)) {
+                    return intPow({ a, b }, int(ex->value));
+                }
+            }
+            if (p.realExponent) {
+                return realPow({ a, b }, c);
+            }
 
             auto arg = arena.allocate(ir::Atan2{ b, a });
             auto aa = arena.allocate(ir::Multiplication{ a, a });
@@ -140,6 +168,60 @@ namespace mnd
 
             return { newA, newB };
         }
+
+        NodePair intPow(NodePair val, int exponent) {
+            auto [a, b] = val;
+
+            if (exponent < 0) {
+                // TODO implement
+                exponent = 0;
+                //return arena.allocate(ir::Division{ one });
+            }
+
+            if (exponent == 0)
+                return { one, zero };
+            else if (exponent == 1)
+                return val;
+            else if (exponent == 2)
+                return sq(a, b);
+            else {
+                bool isEven = (exponent % 2) == 0;
+                if (isEven) {
+                    NodePair square = sq(a, b);
+                    return intPow(square, exponent / 2);
+                }
+                else {
+                    int expm1 = exponent - 1;
+                    NodePair square = sq(a, b);
+                    auto[pa, pb] = intPow(square, expm1 / 2);
+                    return multiplication(pa, pb, a, b);
+                }
+            }
+        }
+
+        NodePair realPow(NodePair val, Node* exponent) {
+            auto [a, b] = val;
+
+            auto arg = arena.allocate(ir::Atan2{ b, a });
+            auto aa = arena.allocate(ir::Multiplication{ a, a });
+            auto bb = arena.allocate(ir::Multiplication{ b, b });
+            auto absSq = arena.allocate(ir::Addition{ aa, bb });
+
+            auto halfc = arena.allocate(ir::Multiplication{ exponent, half });
+
+            auto abspowc = arena.allocate(ir::Pow{ absSq, halfc });
+
+            auto newAbs = arena.allocate(ir::Multiplication{ abspowc, exponent });
+            auto newArg = arena.allocate(ir::Addition{ arg, exponent });
+
+            auto cosArg = arena.allocate(ir::Cos{ newArg });
+            auto sinArg = arena.allocate(ir::Sin{ newArg });
+
+            auto newA = arena.allocate(ir::Multiplication{ cosArg, newAbs });
+            auto newB = arena.allocate(ir::Multiplication{ sinArg, newAbs });
+
+            return { newA, newB };
+        }
     };
 
     ir::Formula expand(const mnd::IterationFormula& fmla)
@@ -182,6 +264,10 @@ std::string mnd::ir::Formula::toString(void) const
             return "(" + std::visit(*this, *n.left) + ") * (" + std::visit(*this, *n.right) + ")";
         }
 
+        std::string operator()(const ir::Division& n) {
+            return "(" + std::visit(*this, *n.left) + ") / (" + std::visit(*this, *n.right) + ")";
+        }
+
         std::string operator()(const ir::Atan2& n) {
             return "atan2(" + std::visit(*this, *n.left) + ", " + std::visit(*this, *n.right) + ")";
         }

mandelbench/CMakeLists.txt

@@ -1,5 +1,8 @@
 cmake_minimum_required(VERSION 3.9)
 
+set(Boost_USE_STATIC_LIBS ON)
+find_package(Boost 1.53)
+
 add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/../libmandel ${CMAKE_CURRENT_BINARY_DIR}/libmandel)
 set(CMAKE_CXX_STANDARD 17)
 
@@ -9,3 +12,8 @@ target_include_directories(mandelbench PUBLIC ../libmandel/include)
 
 target_link_libraries(mandelbench mandel)
 
+if(Boost_FOUND)
+    target_compile_definitions(mandelbench PUBLIC WITH_BOOST)
+    target_include_directories(mandelbench PRIVATE ${Boost_INCLUDE_DIRS})
+endif(Boost_FOUND)
+

mandelbench/mandelbench.cpp

@@ -7,7 +7,7 @@
 
 #include "Fixed.h"
 
-constexpr mnd::MandelViewport benchViewport(void)
+mnd::MandelViewport benchViewport(void)
 {
     return mnd::MandelViewport{ -1.250000598933854152929, 0.0001879894057291665530, 0.0000003839916666666565, 0.0000003839916666666565 };
 }
@@ -72,7 +72,7 @@ std::pair<long long, std::chrono::nanoseconds> measureMips(const std::function<s
 }
 
 
-double benchmark(mnd::Generator& generator)
+double benchmark(mnd::MandelGenerator& generator)
 {
     /*mnd::MandelInfo mi;
     mi.bWidth = 250;
@@ -97,6 +97,7 @@ double benchmark(mnd::Generator& generator)
 
     const mnd::MandelInfo& mi = benches[(testIndex >= benches.size()) ? (benches.size() - 1) : testIndex];
     auto data = std::make_unique<float[]>(mi.bWidth * mi.bHeight);
+
     auto [iters, time] = measureMips([&generator, &mi, &data]() { generator.generate(mi, data.get()); return std::make_pair(data.get(), mi.bWidth * mi.bHeight);  });
     //printf("bench time %d ms\n", time.count() / 1000 / 1000);
     //fflush(stdout);
@@ -120,20 +121,20 @@ int main()
 
     std::cout << "Benchmarking CPU [" << mc.getCpuInfo().getBrand() << "]" << std::endl;
 
-    REPORT_PERFORMANCE("float [MIps]: ", benchmark(mc.getCpuGeneratorFloat()));
-    REPORT_PERFORMANCE("double [MIps]: ", benchmark(mc.getCpuGeneratorDouble()));
-    REPORT_PERFORMANCE("fixed-point 128 bit [MIps]: ", benchmark(mc.getCpuGenerator128()));
+    REPORT_PERFORMANCE("float [MIps]: ", benchmark(*mc.getCpuGenerator(mnd::GeneratorType::FLOAT)));
+    REPORT_PERFORMANCE("double [MIps]: ", benchmark(*mc.getCpuGenerator(mnd::GeneratorType::DOUBLE)));
+    REPORT_PERFORMANCE("fixed-point 128 bit [MIps]: ", benchmark(*mc.getCpuGenerator(mnd::GeneratorType::FIXED128)));
     
 
     for (auto& device : mc.getDevices()) {
         std::cout << "Benchmarking Device [" << device.getName() << "]" << std::endl;
-        if (mnd::Generator* gpuf; gpuf = device.getGeneratorFloat()) {
+        if (mnd::MandelGenerator* gpuf; gpuf = device.getGenerator(mnd::GeneratorType::FLOAT)) {
             REPORT_PERFORMANCE("float [MIps]: ", benchmark(*gpuf));
         }
-        if (mnd::Generator* gpud; gpud = device.getGeneratorDouble()) {
+        if (mnd::MandelGenerator* gpud; gpud = device.getGenerator(mnd::GeneratorType::DOUBLE)) {
             REPORT_PERFORMANCE("double [MIps]: ", benchmark(*gpud));
         }
-        if (mnd::Generator* gpu128; gpu128 = device.getGenerator128()) {
+        if (mnd::MandelGenerator* gpu128; gpu128 = device.getGenerator(mnd::GeneratorType::FIXED128)) {
             REPORT_PERFORMANCE("fixed-point 128 bit [MIps]: ", benchmark(*gpu128));
         }
     }