adding hex-double header

libmandel/include/HexDouble.h

@@ -0,0 +1,121 @@
+#ifndef MANDEL_HEXDOUBLE_H
+#define MANDEL_HEXDOUBLE_H
+
+#include "PolyfloatUtil.h"
+
+namespace mnd
+{
+    struct HexDouble;
+}
+
+struct mnd::HexDouble
+{
+    double x[6];
+
+    inline HexDouble(double val) :
+        x{ val, 0.0, 0.0, 0.0, 0.0, 0.0 }
+    {}
+
+    inline HexDouble(double a, double b, double c, double d, double e, double f) :
+        x{ a, b, c, d, e, f }
+    {}
+
+    double operator[]           (int i) const   { return x[i]; }
+    const double& operator[]    (int i)         { return x[i]; }
+};
+
+
+inline mnd::HexDouble operator+(const mnd::HexDouble& a,
+    const mnd::HexDouble& b)
+{
+    auto[a0, a1] = mnd::pfu::twoSum(a[0], b[0]);
+    auto[b0, b1] = mnd::pfu::twoSum(a[1], b[1]);
+    auto[c0, c1] = mnd::pfu::twoSum(a[2], b[2]);
+    auto[d0, d1] = mnd::pfu::twoSum(a[3], b[3]);
+    auto[e0, e1] = mnd::pfu::twoSum(a[4], b[4]);
+
+    auto t0 = a0;
+    auto [t1, p1] = mnd::pfu::twoSum(a1, b0);
+    auto [t2, p2, p3] = mnd::pfu::threeSum(b1, c0, p1);
+    auto [t3, p4, p5, p6] = mnd::pfu::fourSum(c1, d0, p2, p3);
+    auto [t4, p7] = mnd::pfu::fiveTwoSum(d1, e0, p4, p5, p6);
+    auto t5 = a[5] + b[5] + e1 + p7;
+
+    auto[re0, er1] = mnd::pfu::quickTwoSum(t0, t1);
+    auto[re1, e2] = mnd::pfu::quickTwoSum(er1, t2);
+    auto[re2, e3] = mnd::pfu::quickTwoSum(e2, t3);
+    auto[re3, e4] = mnd::pfu::quickTwoSum(e3, t4);
+    auto[re4, re5] = mnd::pfu::quickTwoSum(e4, t5);
+
+    //return { t1, t2, t3, t4, t5, t6 };
+    return { re0, re1, re2, re3, re4, re5 };
+}
+
+
+inline bool operator>(const mnd::HexDouble& a, const mnd::HexDouble& 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;
+    return a[2] > b[2];
+}
+
+
+inline mnd::HexDouble operator-(const mnd::HexDouble& a,
+    const mnd::HexDouble& b)
+{
+    return a + mnd::HexDouble{ -b[0], -b[1], -b[2], -b[3], -b[4], -b[5] };
+}
+
+
+inline mnd::HexDouble operator*(const mnd::HexDouble& a,
+    const mnd::HexDouble& b)
+{
+    auto[p1_0, p2_0] = mnd::pfu::twoProd(a[0], b[0]);
+
+    auto[p2_1, p3_0] = mnd::pfu::twoProd(a[0], b[1]);
+    auto[p2_2, p3_1] = mnd::pfu::twoProd(a[1], b[0]);
+
+    auto[p3_2, p4_0] = mnd::pfu::twoProd(a[2], b[0]);
+    auto[p3_3, p4_1] = mnd::pfu::twoProd(a[1], b[1]);
+    auto[p3_4, p4_2] = mnd::pfu::twoProd(a[0], b[2]);
+
+    auto[p4_3, p5_0] = mnd::pfu::twoProd(a[3], b[0]);
+    auto[p4_4, p5_1] = mnd::pfu::twoProd(a[2], b[1]);
+    auto[p4_5, p5_2] = mnd::pfu::twoProd(a[1], b[2]);
+    auto[p4_6, p5_3] = mnd::pfu::twoProd(a[0], b[3]);
+
+    auto[p5_4, p6_0] = mnd::pfu::twoProd(a[4], b[0]);
+    auto[p5_5, p6_1] = mnd::pfu::twoProd(a[3], b[1]);
+    auto[p5_6, p6_2] = mnd::pfu::twoProd(a[2], b[2]);
+    auto[p5_7, p6_3] = mnd::pfu::twoProd(a[1], b[3]);
+    auto[p5_8, p6_4] = mnd::pfu::twoProd(a[0], b[4]);
+
+    auto t1 = p1_0;
+    auto[t2, tl3, tl4] = mnd::pfu::threeSum(p2_0, p2_1, p2_2);
+    auto[t3, tl4_2, tl5] = mnd::pfu::sixThreeSum(p3_0, p3_1, p3_2, p3_3, p3_4, tl3);
+    auto[t4, tl5_2, tl6] = mnd::pfu::nineThreeSum(p4_0, p4_1, p4_2, p4_3, p4_4, p4_5, p4_6, tl4, tl4_2);
+    auto[x1, x2, x3] = mnd::pfu::nineThreeSum(p5_0, p5_1, p5_2, p5_3, p5_4, p5_5, p5_6, p5_7, p5_8);
+    auto[t5, tl6_1, tl7] = mnd::pfu::sixThreeSum(x1, x2, x3, tl5, tl5_2, 0.0);
+
+    auto t6 = tl6 + tl6_1 + tl7 + p6_0 + p6_1 + p6_2 + p6_3 + p6_4 +
+        a[5] * b[0] + a[4] * b[1] + a[3] * b[2] + a[2] * b[3] + a[1] * b[4] + a[0] * b[5];
+
+    auto[re0, e1] = mnd::pfu::quickTwoSum(t1, t2);
+    auto[re1, e2] = mnd::pfu::quickTwoSum(e1, t3);
+    auto[re2, e3] = mnd::pfu::quickTwoSum(e2, t4);
+    auto[re3, e4] = mnd::pfu::quickTwoSum(e3, t5);
+    auto[re4, re5] = mnd::pfu::quickTwoSum(e4, t6);
+
+    //return { t1, t2, t3, t4, t5, t6 };
+    return { re0, re1, re2, re3, re4, re5 };
+}
+
+
+#endif // MANDEL_HEXDOUBLE_H
+