implementing gen

libmandel/include/LightDoubleDouble.h

@@ -6,26 +6,33 @@
 namespace mnd
 {
     struct LightDoubleDouble;
+    struct LightDoubleFloat;
 
     namespace ldd
     {
-        using DoublePair = std::pair<double, double>;
-        inline static DoublePair quickTwoSum(double a, double b)
+        template<typename T>
+        using Pair = std::pair<T, T>;
+        using DoublePair = Pair<double>;
+        using FloatPair = Pair<float>;
+
+        template<typename T>
+        inline static Pair<T> quickTwoSum(T a, T b)
         {
-            double s = a + b;
-            double e = b - (s - a);
+            T s = a + b;
+            T e = b - (s - a);
             return { s, e };
         }
 
-        inline DoublePair twoSum(double a, double b)
+        template<typename T>
+        inline Pair<T> twoSum(T a, T b)
         {
-            double s = a + b;
-            double v = s - a;
-            double e = (a - (s - v)) + (b - v);
+            T s = a + b;
+            T v = s - a;
+            T e = (a - (s - v)) + (b - v);
             return { s, e };
         }
 
-        inline DoublePair split(double a)
+        inline Pair<double> split(double a)
         {
             static const double splitter = double((1ULL << 27) + 1);
             double t = splitter * a;
@@ -34,19 +41,21 @@ namespace mnd
             return { ahi, alo };
         }
 
-        inline DoublePair twoProd(double a, double b)
+        template<typename T>
+        inline Pair<T> twoProd(T a, T b)
         {
-            double p = a * b;
+            T p = a * b;
             auto [ahi, alo] = split(a);
             auto [bhi, blo] = split(b);
-            double e = ((ahi * bhi - p) + ahi * blo + alo * bhi) + alo * blo;
+            T e = ((ahi * bhi - p) + ahi * blo + alo * bhi) + alo * blo;
             return { p, e };
         }
 
-        inline DoublePair twoProdFma(double a, double b)
+        template<typename T>
+        inline Pair<T> twoProdFma(T a, T b)
         {
-            double p = a * b;
-            double e = std::fma(a, b, -p);
+            T p = a * b;
+            T e = std::fma(a, b, -p);
             return { p, e };
         }
     }
@@ -102,4 +111,53 @@ inline mnd::LightDoubleDouble operator*(const mnd::LightDoubleDouble& a,
 }
 
 
+struct mnd::LightDoubleFloat
+{
+    float x[2];
+
+    inline LightDoubleFloat(double val)
+    {
+        x[0] = float(val);
+        x[1] = float(val - double(x[0]));
+    }
+
+    inline LightDoubleFloat(float u, float l) :
+        x{ u, l }
+    {}
+
+    inline LightDoubleFloat(mnd::ldd::FloatPair dp) :
+        x{ dp.first, dp.second }
+    {}
+
+    float operator[]           (int i) const   { return x[i]; }
+    const float& operator[]    (int i)         { return x[i]; }
+
+private:
+};
+
+
+inline mnd::LightDoubleFloat operator+(const mnd::LightDoubleFloat& a,
+    const mnd::LightDoubleFloat& b)
+{
+    auto[s, e] = mnd::ldd::twoSum(a[0], b[0]);
+    e += a[1] + b[1];
+    return mnd::ldd::quickTwoSum(s, e);
+}
+
+inline mnd::LightDoubleFloat operator-(const mnd::LightDoubleFloat& a,
+    const mnd::LightDoubleFloat& b)
+{
+    auto[s, e] = mnd::ldd::twoSum(a[0], -b[0]);
+    e += a[1] - b[1];
+    return mnd::ldd::quickTwoSum(s, e);
+}
+
+inline mnd::LightDoubleFloat operator*(const mnd::LightDoubleFloat& a,
+    const mnd::LightDoubleFloat& b)
+{
+    auto[p1, p2] = mnd::ldd::twoProd(a[0], b[0]);
+    p2 += a[0] * b[1] + a[1] * b[0];
+    return mnd::ldd::quickTwoSum(p1, p2);
+}
+
 #endif // MANDEL_LIGHTDOUBLEDOUBLE_H

libmandel/src/ClGenerators.cpp

@@ -238,63 +238,22 @@ void ClGeneratorDoubleFloat::generate(const mnd::MandelInfo& info, float* data)
     };
 
     Buffer buffer_A(context, CL_MEM_WRITE_ONLY, bufferSize);
-    double pixelScX = double(info.view.width / info.bWidth);
-    double pixelScY = double(info.view.height / info.bHeight);
-
-    auto[x1, x2] = splitDouble(double(info.view.x));
-    auto[y1, y2] = splitDouble(double(info.view.y));
-    auto[w1, w2] = splitDouble(pixelScX);
-    auto[h1, h2] = splitDouble(pixelScY);
-
-    /*
-    for (int px = 0; px < info.bWidth; px++) {
-        for (int py = 0; py < info.bHeight; py++) {
-            std::pair<float, float> xl = { x1, x2 };
-            std::pair<float, float> yt = { y1, y2 };
-            std::pair<float, float> pixelScaleX = { w1, w2 };
-            std::pair<float, float> pixelScaleY = { h1, h2 };
-
-            std::pair<float, float> a = add(mulFloat(pixelScaleX, (float) px), xl); // pixelScaleX * px + xl
-            std::pair<float, float> b = add(mulFloat(pixelScaleY, (float) py), yt); // pixelScaleY * py + yt
-            std::pair<float, float> ca = a;
-            std::pair<float, float> cb = b;
-
-            int n = 0;
-            while (n < info.maxIter - 1) {
-                std::pair<float, float> aa = mul(a, a);
-                std::pair<float, float> bb = mul(b, b);
-                std::pair<float, float> ab = mul(a, b);
-                if (aa.first + bb.first > 16) break;
-                std::pair<float, float> minusbb = { -bb.first, -bb.second };
-                a = add(add(aa, minusbb), ca);
-                b = add(add(ab, ab), cb);
-                n++;
-            }
-
-            // N + 1 - log (log  |Z(N)|) / log 2
-            if (n >= info.maxIter - 1)
-                data[px + py * info.bWidth] = info.maxIter;
-            else {
-                if (info.smooth)
-                    data[px + py * info.bWidth] = ((float) n) + 1 - log(log(a.first * a.first + b.first * b.first ) / 2) / log(2.0f);
-                else
-                    data[px + py * info.bWidth] = ((float)n);
-            }
-        }
-    }
-    return;
-    */
-    
+    mnd::LightDoubleFloat pixelScX = double(info.view.width / info.bWidth);
+    mnd::LightDoubleFloat pixelScY = double(info.view.height / info.bHeight);
+
+    mnd::LightDoubleFloat x = double(info.view.x);
+    mnd::LightDoubleFloat y = double(info.view.y);
+
     kernel.setArg(0, buffer_A);
     kernel.setArg(1, int(info.bWidth));
-    kernel.setArg(2, x1);
-    kernel.setArg(3, x2);
-    kernel.setArg(4, y1);
-    kernel.setArg(5, y2);
-    kernel.setArg(6, w1);
-    kernel.setArg(7, w2);
-    kernel.setArg(8, h1);
-    kernel.setArg(9, h2);
+    kernel.setArg(2, x[0]);
+    kernel.setArg(3, x[1]);
+    kernel.setArg(4, y[0]);
+    kernel.setArg(5, y[1]);
+    kernel.setArg(6, pixelScX[0]);
+    kernel.setArg(7, pixelScX[1]);
+    kernel.setArg(8, pixelScY[0]);
+    kernel.setArg(9, pixelScY[1]);
     kernel.setArg(10, int(info.maxIter));
     kernel.setArg(11, int(info.smooth ? 1 : 0));
 

libmandel/src/IterationGenerator.cpp

@@ -183,15 +183,30 @@ void NaiveIRGenerator<U>::generate(const mnd::MandelInfo& info, float* data)
 
 
 template<typename T>
+struct EvalNode
+{
+
+};
+
+
+template<typename T>
 double NaiveIRGenerator<T>::calc(mnd::ir::Node* expr, double a, double b, double x, double y)
 {
     struct DoubleVisitor
     {
         double a, b, x, y;
         double visitNode(ir::Node* n) {
+            auto& nodeData = getNodeData(n);
+            if (std::get_a)
             return std::visit(*this, *n);
         }
 
+        std::any& getNodeData(ir::Node* n) {
+            return std::visit([](auto& n) {
+                return n.nodeData;
+            }, *n);
+        }
+
         double operator()(const ir::Constant& c) {
             return mnd::convert<double>(c.value);
         }