Nicolas Winkler 5 лет назад
Родитель
Сommit
c506f398a5

+ 14 - 0
libmandel/include/IterationGenerator.h

@@ -4,6 +4,7 @@
 #include "Generators.h"
 #include "ClGenerators.h"
 #include "IterationFormula.h"
+#include "IterationIR.h"
 
 #include <utility>
 #include <complex>
@@ -13,6 +14,7 @@ namespace mnd
     class IterationGenerator;
 
     class NaiveGenerator;
+    class NaiveIRGenerator;
     class CompiledGenerator;
     class CompiledClGenerator;
 
@@ -45,6 +47,18 @@ private:
 };
 
 
+class mnd::NaiveIRGenerator : public mnd::MandelGenerator
+{
+    const ir::Formula& form;
+public:
+    NaiveIRGenerator(const ir::Formula& irf, const mnd::Real& prec);
+    NaiveIRGenerator(NaiveIRGenerator&&) = default;
+
+    virtual void generate(const MandelInfo& info, float* data);
+    double calc(ir::Node* expr, double a, double b, double x, double y);
+};
+
+
 #if defined(__x86_64__) || defined(_M_X64)
 class mnd::CompiledGenerator : public mnd::MandelGenerator
 {

+ 29 - 8
libmandel/src/IterationCompiler.cpp

@@ -106,6 +106,13 @@ namespace mnd
             return res;
         }
 
+        static double myAtan2(double y, double x)
+        {
+            double result = ::atan2(y, x);
+            printf("atan2(%f, %f) = %f\n", y, x, result);
+            return result;
+        }
+
         Reg operator()(const ir::Atan2& at2) {
             using namespace asmjit;
             auto y = visitNode(*at2.left);
@@ -114,7 +121,7 @@ namespace mnd
             auto arg = cc.newXmmSd();
             double(*atanFunc)(double, double) = ::atan2;
             cc.comment("call atan2");
-            auto call = cc.call(imm(atanFunc), FuncSignatureT<double, double, double>(CallConv::kIdHostCDecl));
+            auto call = cc.call(imm(atanFunc), FuncSignatureT<double, double, double>(CallConv::kIdHost));
             call->setArg(0, y);
             call->setArg(1, x);
             call->setRet(0, arg);
@@ -129,7 +136,7 @@ namespace mnd
             auto arg = cc.newXmmSd();
             double(*powFunc)(double, double) = ::pow;
             cc.comment("call pow");
-            auto call = cc.call(imm(powFunc), FuncSignatureT<double, double, double>(CallConv::kIdHostCDecl));
+            auto call = cc.call(imm(powFunc), FuncSignatureT<double, double, double>(CallConv::kIdHost));
             call->setArg(0, a);
             call->setArg(1, b);
             call->setRet(0, arg);
@@ -143,7 +150,7 @@ namespace mnd
             auto arg = cc.newXmmSd();
             double(*cosFunc)(double) = ::cos;
             cc.comment("call cos");
-            auto call = cc.call(imm(cosFunc), FuncSignatureT<double, double>(CallConv::kIdHostCDecl));
+            auto call = cc.call(imm(cosFunc), FuncSignatureT<double, double>(CallConv::kIdHost));
             call->setArg(0, a);
             call->setRet(0, arg);
             return arg;
@@ -156,7 +163,7 @@ namespace mnd
             auto arg = cc.newXmmSd();
             double(*sinFunc)(double) = ::sin;
             cc.comment("call sin");
-            auto call = cc.call(imm(sinFunc), FuncSignatureT<double, double>(CallConv::kIdHostCDecl));
+            auto call = cc.call(imm(sinFunc), FuncSignatureT<double, double>(CallConv::kIdHost));
             call->setArg(0, a);
             call->setRet(0, arg);
             return arg;
@@ -169,7 +176,7 @@ namespace mnd
             auto arg = cc.newXmmSd();
             double(*expFunc)(double) = ::exp;
             cc.comment("call exp");
-            auto call = cc.call(imm(expFunc), FuncSignatureT<double, double>(CallConv::kIdHostCDecl));
+            auto call = cc.call(imm(expFunc), FuncSignatureT<double, double>(CallConv::kIdHost));
             call->setArg(0, a);
             call->setRet(0, arg);
             return arg;
@@ -182,7 +189,7 @@ namespace mnd
             auto arg = cc.newXmmSd();
             double(*logFunc)(double) = ::log;
             cc.comment("call log");
-            auto call = cc.call(imm(logFunc), FuncSignatureT<double, double>(CallConv::kIdHostCDecl));
+            auto call = cc.call(imm(logFunc), FuncSignatureT<double, double>(CallConv::kIdHost));
             call->setArg(0, a);
             call->setRet(0, arg);
             return arg;
@@ -190,6 +197,10 @@ namespace mnd
     };
 
 
+    static void printD(double d) {
+        printf("val: %f\n", d); fflush(stdout);
+    }
+
     CompiledGenerator compile(const ir::Formula& formula)
     {
         using namespace asmjit;
@@ -243,9 +254,11 @@ namespace mnd
         comp.movapd(bb, b);
         comp.mulsd(bb, b);
         comp.addsd(bb, aa);
+        //auto call = comp.call(imm(printD), FuncSignatureT<void, double>(CallConv::kIdHost));
+        //call->setArg(0, bb);
 
         comp.comisd(bb, sixteen);
-        comp.jle(endLoop);
+        comp.jnb(endLoop);
 
         comp.inc(k);
         comp.cmp(k, maxIter);
@@ -417,13 +430,21 @@ namespace mnd
         const IterationFormula& z0,
         const IterationFormula& zi)
     {
+
         //std::unique_ptr<mnd::MandelGenerator> ng = std::make_unique<NaiveGenerator>(z0.clone(), zi.clone(), mnd::getPrecision<double>());
+        
+        IterationFormula z0o = z0.clone();
+        IterationFormula zio = zi.clone();
+        z0o.optimize();
+        zio.optimize();
 
-        ir::Formula irf = mnd::expand(z0, zi);
+        ir::Formula irf = mnd::expand(z0o, zio);
         irf.optimize();
         printf("ir: %s\n", irf.toString().c_str()); fflush(stdout);
         auto dg = std::make_unique<CompiledGenerator>(compile(irf));
         printf("asm: %s\n", dg->dump().c_str()); fflush(stdout);
+        
+        //auto dg = std::make_unique<NaiveIRGenerator>(*irf, mnd::getPrecision<double>());
 
         std::vector<std::unique_ptr<mnd::MandelGenerator>> vec;
         //vec.push_back(std::move(ng));

+ 123 - 1
libmandel/src/IterationGenerator.cpp

@@ -9,6 +9,7 @@
 
 using mnd::IterationGenerator;
 using mnd::NaiveGenerator;
+using mnd::NaiveIRGenerator;
 using mnd::IterationFormula;
 
 
@@ -124,6 +125,127 @@ std::complex<double> NaiveGenerator::calc(mnd::Expression& expr, std::complex<do
 }
 
 
+NaiveIRGenerator::NaiveIRGenerator(const ir::Formula& irf,
+                                   const mnd::Real& prec) :
+    mnd::MandelGenerator{ prec },
+    form{ irf }
+{
+}
+
+
+void NaiveIRGenerator::generate(const mnd::MandelInfo& info, float* data)
+{
+    const MandelViewport& view = info.view;
+
+    const bool parallel = true;
+    using T = double;
+    T viewx = mnd::convert<T>(view.x);
+    T viewy = mnd::convert<T>(view.y);
+    T wpp = mnd::convert<T>(view.width / info.bWidth);
+    T hpp = mnd::convert<T>(view.height / info.bHeight);
+
+    if constexpr (parallel)
+        omp_set_num_threads(omp_get_num_procs());
+//#pragma omp parallel for schedule(static, 1) if (parallel)
+    for (long j = 0; j < info.bHeight; j++) {
+        T y = viewy + T(double(j)) * hpp;
+        long i = 0;
+        for (i; i < info.bWidth; i++) {
+            T x = viewx + T(double(i)) * wpp;
+
+            T a = calc(form.startA, x, y, 0, 0);
+            T b = calc(form.startB, x, y, 0, 0);
+
+            int k = 0;
+            for (k = 0; k < info.maxIter; k++) {
+                double newA = calc(form.newA, a, b, x, y);
+                double newB = calc(form.newB, a, b, x, y);
+                a = newA;
+                b = newB;
+                if (a * a + b * b >= 16.0)
+                    break;
+            }
+            data[i + j * info.bWidth] = float(k);
+        }
+    }
+}
+
+
+double NaiveIRGenerator::calc(ir::Node* expr, double a, double b, double x, double y)
+{
+    struct DoubleVisitor
+    {
+        double a, b, x, y;
+        double visitNode(ir::Node* n) {
+            return std::visit(*this, *n);
+        }
+
+        double operator()(const ir::Constant& c) {
+            return mnd::convert<double>(c.value);
+        }
+
+        double operator()(const ir::Variable& v) {
+            if (v.name == "z_re")
+                return a;
+            else if (v.name == "z_im")
+                return b;
+            else if (v.name == "c_re")
+                return x;
+            else if (v.name == "c_im")
+                return y;
+            else
+                return 0.0;
+        }
+
+        double operator()(const ir::Negation& n) {
+            return -visitNode(n.value);
+        }
+
+        double operator()(const ir::Addition& n) {
+            return visitNode(n.left) + visitNode(n.right);
+        }
+
+        double operator()(const ir::Subtraction& n) {
+            return visitNode(n.left) - visitNode(n.right);
+        }
+
+        double operator()(const ir::Multiplication& n) {
+            return visitNode(n.left) * visitNode(n.right);
+        }
+
+        double operator()(const ir::Division& n) {
+            return visitNode(n.left) / visitNode(n.right);
+        }
+
+        double operator()(const ir::Atan2& n) {
+            return ::atan2(visitNode(n.left), visitNode(n.right));
+        }
+
+        double operator()(const ir::Pow& n) {
+            return ::pow(visitNode(n.left), visitNode(n.right));
+        }
+
+        double operator()(const ir::Cos& n) {
+            return ::cos(visitNode(n.value));
+        }
+
+        double operator()(const ir::Sin& n) {
+            return ::sin(visitNode(n.value));
+        }
+
+        double operator()(const ir::Exp& n) {
+            return ::exp(visitNode(n.value));
+        }
+
+        double operator()(const ir::Ln& n) {
+            return ::log(visitNode(n.value));
+        }
+    };
+
+    DoubleVisitor dv{ a, b, x, y };
+    return dv.visitNode(expr);
+}
+
 
 using mnd::CompiledGenerator;
 using mnd::CompiledClGenerator;
@@ -172,7 +294,7 @@ void CompiledGenerator::generate(const mnd::MandelInfo& info, float* data)
     using IterFunc = int (*)(double, double, int);
 
     omp_set_num_threads(omp_get_num_procs());
-//#pragma omp parallel for schedule(static, 1)
+#pragma omp parallel for schedule(static, 1)
     for (int i = 0; i < info.bHeight; i++) {
         double y = mnd::convert<double>(info.view.y + info.view.height * i / info.bHeight);
         for (int j = 0; j < info.bWidth; j++) {