removed dependency to qd_real #18

libmandel/include/QuadDouble.h

@@ -0,0 +1,119 @@
+#ifndef MANDEL_QUADDOUBLE_H
+#define MANDEL_QUADDOUBLE_H
+
+#include "PolyfloatUtil.h"
+
+namespace mnd
+{
+    struct QuadDouble;
+
+    inline mnd::pfu::Quadruple<double> renorm2(double x0, double x1, double x2, double x3, double x4)
+    {
+        auto [st1, t4] = mnd::pfu::quickTwoSum(x3, x4);
+        auto [st2, t3] = mnd::pfu::quickTwoSum(x2, st1);
+        auto [st3, t2] = mnd::pfu::quickTwoSum(x1, st2);
+        auto [t0, t1] = mnd::pfu::quickTwoSum(x0, st3);
+
+        auto e = t0;
+
+        auto [r0, e1] = mnd::pfu::quickTwoSum(e, t1);
+        auto [r1, e2] = mnd::pfu::quickTwoSum(e1, t2);
+        auto [r2, e3] = mnd::pfu::quickTwoSum(e2, t3);
+        auto r3 = e3 + t4;
+        return { r0, r1, r2, r3 };
+    }
+}
+
+struct mnd::QuadDouble
+{
+    double x[4];
+
+    inline QuadDouble(double val) :
+        x{ val, 0.0, 0.0, 0.0 }
+    {}
+
+    inline QuadDouble(double a, double b, double c, double d) :
+        x{ a, b, c, d }
+    {}
+
+    double operator[]           (int i) const   { return x[i]; }
+    const double& operator[]    (int i)         { return x[i]; }
+};
+
+
+inline mnd::QuadDouble operator+(const mnd::QuadDouble& a,
+    const mnd::QuadDouble& b)
+{
+    auto[s0, e0] = mnd::pfu::twoSum(a[0], b[0]);
+    auto[s1, e1] = mnd::pfu::twoSum(a[1], b[1]);
+    auto[s2, e2] = mnd::pfu::twoSum(a[2], b[2]);
+    auto[s3, e3] = mnd::pfu::twoSum(a[3], b[3]);
+    auto r0 = s0;
+
+    auto [r1, t0] = mnd::pfu::twoSum(s1, e0);
+    auto [r2, t1, t2] = mnd::pfu::threeSum(s2, e1, t0);
+    auto [r3, t3, _t4] = mnd::pfu::threeSum(s3, e2, t1);
+    auto [r4, _t5, _t6] = mnd::pfu::threeSum(e3, t3, t2);
+
+    auto[re0, er1] = mnd::pfu::quickTwoSum(r0, r1);
+    auto[re1, er2] = mnd::pfu::quickTwoSum(er1, r2);
+    auto[re2, er3] = mnd::pfu::quickTwoSum(er2, r3);
+    auto[re3, er4] = mnd::pfu::quickTwoSum(er3, r4);
+    return { re0, re1, re2, re3 };
+}
+
+
+inline bool operator>(const mnd::QuadDouble& a, const mnd::QuadDouble& b)
+{
+    if (a[0] > b[0])
+        return true;
+    if (a[0] < b[0])
+        return false;
+    if (a[1] > b[1])
+        return true;
+    if (a[1] < b[1])
+        return false;
+    if (a[2] > b[2])
+        return true;
+    if (a[2] < b[2])
+        return false;
+    return a[3] > b[3];
+}
+
+
+inline mnd::QuadDouble operator-(const mnd::QuadDouble& a,
+    const mnd::QuadDouble& b)
+{
+    return a + mnd::QuadDouble{ -b[0], -b[1], -b[2], -b[3] };
+}
+
+
+inline mnd::QuadDouble operator*(const mnd::QuadDouble& a,
+    const mnd::QuadDouble& b)
+{
+    auto[a0, b0] = mnd::pfu::twoProd(a[0], b[0]);
+    auto[b1, c0] = mnd::pfu::twoProd(a[0], b[1]);
+    auto[b2, c1] = mnd::pfu::twoProd(a[1], b[0]);
+    auto[c2, d0] = mnd::pfu::twoProd(a[0], b[2]);
+    auto[c3, d1] = mnd::pfu::twoProd(a[1], b[1]);
+    auto[c4, d2] = mnd::pfu::twoProd(a[2], b[0]);
+
+    auto d5 = a[3] * b[0];
+    auto d6 = a[2] * b[1];
+    auto d7 = a[1] * b[2];
+    auto d8 = a[0] * b[3];
+
+    auto r0 = a0;
+    auto[r1, c5, d3] = mnd::pfu::threeSum(b0, b1, b2);
+    auto[r2, d4, e0] = mnd::pfu::sixThreeSum(c0, c1, c2, c3, c4, c5);
+    auto[r3, e1, ex_] = mnd::pfu::nineThreeSum(d0, d1, d2, d3, d4, d5, d6, d7, d8);
+    auto r4 = e0 + e1;
+
+    auto [n0, n1, n2, n3] = mnd::renorm2(r0, r1, r2, r3, r4);
+
+    return { n0, n1, n2, n3 };
+}
+
+
+#endif // MANDEL_QUADDOUBLE_H
+

libmandel/include/Types.h

@@ -21,10 +21,11 @@
 #   include <boost/functional/hash.hpp>
 #endif
 
-#include <qd/dd_real.h>
-#include <qd/qd_real.h>
+//#include <qd/dd_real.h>
+//#include <qd/qd_real.h>
 
 #include "LightDoubleDouble.h"
+#include "QuadDouble.h"
 #include "TripleDouble.h"
 #include "TripleFloat.h"
 #include "HexDouble.h"
@@ -74,7 +75,7 @@ namespace mnd
 #endif
 
     using DoubleDouble = mnd::LightDoubleDouble;
-    using QuadDouble = qd_real;
+    //using QuadDouble = qd_real;
 
     /*inline DoubleDouble abs(const DoubleDouble& x) { return ::abs(x); }
     inline DoubleDouble sqrt(const DoubleDouble& x) { return ::sqrt(x); }
@@ -88,7 +89,7 @@ namespace mnd
     inline DoubleDouble exp(const DoubleDouble& x) { return ::exp(x); }*/
 
 
-    inline QuadDouble abs(const QuadDouble& x) { return ::abs(x); }
+    /*inline QuadDouble abs(const QuadDouble& x) { return ::abs(x); }
     inline QuadDouble sqrt(const QuadDouble& x) { return ::sqrt(x); }
     inline QuadDouble floor(const QuadDouble& x) { return ::floor(x); }
     inline QuadDouble log(const QuadDouble& x) { return ::log(x); }
@@ -97,7 +98,7 @@ namespace mnd
     inline QuadDouble atan2(const QuadDouble& y, const QuadDouble& x) { return ::atan2(y, x); }
     inline QuadDouble cos(const QuadDouble& x) { return ::cos(x); }
     inline QuadDouble sin(const QuadDouble& x) { return ::sin(x); }
-    inline QuadDouble exp(const QuadDouble& x) { return ::exp(x); }
+    inline QuadDouble exp(const QuadDouble& x) { return ::exp(x); }*/
 
     inline double abs(double x) { return ::abs(x); }
     inline float abs(float x) { return ::abs(x); }