iteration formula

Almond.cpp

@@ -7,6 +7,18 @@
 
 #include <cmath>
 
+mnd::MandelGenerator* generateTest()
+{
+    mnd::Variable z{ "z" };
+    mnd::Variable c{ "c" };
+    mnd::Pow m{ std::make_unique<mnd::Expression>(z), std::make_unique<mnd::Expression>(mnd::Constant{ 2.01 }) };
+    //mnd::Multiplication m2{ std::make_unique<mnd::Expression>(std::move(m)), std::make_unique<mnd::Expression>(z) };
+    mnd::Addition a{ std::make_unique<mnd::Expression>(std::move(m)), std::make_unique<mnd::Expression>(c) };
+
+    mnd::IterationFormula itf{ std::make_unique<mnd::Expression>(std::move(a)) };
+    return new mnd::NaiveGenerator(std::move(itf), 1.0e-10);
+}
+
 Almond::Almond(QWidget* parent) :
     QMainWindow{ parent },
     mandelContext{ mnd::initializeContext() }
@@ -26,6 +38,8 @@ Almond::Almond(QWidget* parent) :
     mandelGeneratorSave = &mandelContext.getDefaultGenerator();
     mandelViewSave = mw->getViewport();
 
+    mw->setGenerator(generateTest());
+
     QObject::connect(mw.get(), &MandelWidget::pointSelected, this, &Almond::pointSelected);
     ui.mainContainer->addWidget(mw.get());
     ui.maxIterations->setValidator(new QIntValidator(1, 1000000000, this));

libmandel/include/IterationFormula.h

@@ -1,12 +1,13 @@
 #ifndef MANDEL_ITERATIONFORMULA_H
 #define MANDEL_ITERATIONFORMULA_H
 
+#include <variant>
+#include <memory>
 
 namespace mnd
 {
     struct IterationFormula;
 
-    struct Expression;
     struct Constant;
     struct Variable;
     struct UnaryOperation;
@@ -15,6 +16,18 @@ namespace mnd
     struct Multiplication;
     struct Division;
     struct Pow;
+
+    using Expression = std::variant<
+            Constant,
+            Variable,
+            UnaryOperation,
+            Addition,
+            Multiplication,
+            Division,
+            Pow
+    >;
+
+    class FormulaVisitor;
 }
 
 
@@ -24,30 +37,25 @@ struct mnd::IterationFormula
 };
 
 
-struct mnd::Expression
-{
-};
-
-
-struct mnd::Constant : mnd::Expression
+struct mnd::Constant
 {
     double value;
 };
 
 
-struct mnd::Variable : mnd::Expression
+struct mnd::Variable
 {
     std::string name;
 };
 
 
-struct mnd::UnaryOperation : mnd::Expression
+struct mnd::UnaryOperation
 {
     std::unique_ptr<Expression> operand;
 };
 
 
-struct mnd::BinaryOperation : mnd::Expression
+struct mnd::BinaryOperation
 {
     std::unique_ptr<Expression> left;
     std::unique_ptr<Expression> right;

libmandel/include/IterationGenerator.h

@@ -4,20 +4,36 @@
 #include "Generators.h"
 #include "IterationFormula.h"
 
+#include <utility>
+#include <complex>
 
 namespace mnd
 {
     class IterationGenerator;
+
+    class NaiveGenerator;
 }
 
 
 class mnd::IterationGenerator : public mnd::MandelGenerator
 {
+protected:
     IterationFormula itf;
 public:
     IterationGenerator(IterationFormula itf, const mnd::Real& prec);
 };
 
 
+class mnd::NaiveGenerator : public mnd::IterationGenerator
+{
+public:
+    NaiveGenerator(IterationFormula itf, const mnd::Real& prec);
+
+    virtual void generate(const MandelInfo& info, float* data);
+private:
+    std::complex<double> iterate(std::complex<double> z, std::complex<double> c);
+    std::complex<double> calc(mnd::Expression& expr, std::complex<double> z, std::complex<double> c);
+};
+
 #endif // MANDEL_ITERATIONGENERATOR_H
 

libmandel/include/Mandel.h

@@ -8,6 +8,7 @@
 
 #include "MandelUtil.h"
 #include "Generators.h"
+#include "IterationGenerator.h"
 #include "CpuGenerators.h"
 #include "Hardware.h"
 

libmandel/src/IterationGenerator.cpp

@@ -1,12 +1,14 @@
 #include "IterationGenerator.h"
 
+#include <omp.h>
 
 using mnd::IterationGenerator;
+using mnd::NaiveGenerator;
 using mnd::IterationFormula;
 
 
 IterationGenerator::IterationGenerator(IterationFormula itf,
-                                       const mnd::Real& prec) :
+                                   const mnd::Real& prec) :
     mnd::MandelGenerator{ prec },
     itf{ std::move(itf) }
 {
@@ -14,5 +16,105 @@ IterationGenerator::IterationGenerator(IterationFormula itf,
 
 
 
+NaiveGenerator::NaiveGenerator(IterationFormula itf,
+                                   const mnd::Real& prec) :
+    IterationGenerator{ std::move(itf), prec }
+{
+}
+
+
+void NaiveGenerator::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);
+
+    T juliaX = mnd::convert<T>(info.juliaX);
+    T juliaY = mnd::convert<T>(info.juliaY);
+
+    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 cx = info.julia ? juliaX : x;
+            T cy = info.julia ? juliaY : y;
+            std::complex<double> z{ x, y };
+            std::complex<double> c{ cx, cy };
+
+            int k = 0;
+            for (k = 0; k < info.maxIter; k++) {
+                z = this->iterate(z, c);
+                if (std::abs(z) >= 4)
+                    break;
+            }
+            data[i + j * info.bWidth] = k;
+            /*if (info.smooth) {
+                if (k >= info.maxIter)
+                    data[i + j * info.bWidth] = float(info.maxIter);
+                else {
+                    float aapp = mnd::convert<float>(a);
+                    float bapp = mnd::convert<float>(b);
+                    data[i + j * info.bWidth] = ((float) k) + 1 - ::logf(::logf(aapp * aapp + bapp * bapp) / 2) / ::logf(2.0f);
+                }
+            }
+            else
+                data[i + j * info.bWidth] = k;*/
+        }
+    }
+}
+
+
+std::complex<double> NaiveGenerator::iterate(std::complex<double> z, std::complex<double> c)
+{
+    auto& expr = *itf.expr;
+    return calc(expr, z, c);
+}
+
+
+std::complex<double> NaiveGenerator::calc(mnd::Expression& expr, std::complex<double> z, std::complex<double> c)
+{
+    std::complex<double> result = 0;
+    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;
+        }
+        else if constexpr (std::is_same<T, mnd::Variable>::value) {
+            if (ex.name == "z")
+                result = z;
+            else if (ex.name == "c")
+                result = c;
+        }
+        else if constexpr (std::is_same<T, mnd::UnaryOperation>::value) {
+            result = -calc(*ex.operand, z, c);
+        }
+        else if constexpr (std::is_same<T, mnd::Addition>::value) {
+            if (ex.subtraction)
+                result = calc(*ex.left, z, c) - calc(*ex.right, z, c);
+            else
+                result = calc(*ex.left, z, c) + calc(*ex.right, z, c);
+        }
+        else if constexpr (std::is_same<T, mnd::Multiplication>::value) {
+            result = calc(*ex.left, z, c) * calc(*ex.right, z, c);
+        }
+        else if constexpr (std::is_same<T, mnd::Division>::value) {
+            result = calc(*ex.left, z, c) / calc(*ex.right, z, c);
+        }
+        else if constexpr (std::is_same<T, mnd::Pow>::value) {
+            result = std::pow(calc(*ex.left, z, c), calc(*ex.right, z, c));
+        }
+    }, expr);
+    return result;
+}