OctaDouble.h 6.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194
  1. #ifndef MANDEL_OCTADOUBLE_H
  2. #define MANDEL_OCTADOUBLE_H
  3. #include "PolyfloatUtil.h"
  4. #ifndef _MSC_VER
  5. # define likely(x) __builtin_expect(!!(x), 1)
  6. #else
  7. # define likely(x) (x)
  8. #endif
  9. namespace mnd
  10. {
  11. struct OctaDouble;
  12. }
  13. struct mnd::OctaDouble
  14. {
  15. double x[8];
  16. inline OctaDouble(double val) :
  17. x{ val, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
  18. {}
  19. inline OctaDouble(double a, double b, double c, double d, double e, double f, double g, double h) :
  20. x{ a, b, c, d, e, f, g, h }
  21. {}
  22. double operator[] (int i) const { return x[i]; }
  23. const double& operator[] (int i) { return x[i]; }
  24. };
  25. inline mnd::OctaDouble operator+(const mnd::OctaDouble& a,
  26. const mnd::OctaDouble& b)
  27. {
  28. auto[a0, a1] = mnd::pfu::twoSum(a[0], b[0]);
  29. auto[b0, b1] = mnd::pfu::twoSum(a[1], b[1]);
  30. auto[c0, c1] = mnd::pfu::twoSum(a[2], b[2]);
  31. auto[d0, d1] = mnd::pfu::twoSum(a[3], b[3]);
  32. auto[e0, e1] = mnd::pfu::twoSum(a[4], b[4]);
  33. auto[f0, f1] = mnd::pfu::twoSum(a[5], b[5]);
  34. auto[g0, g1] = mnd::pfu::twoSum(a[6], b[6]);
  35. if (likely(a1 + b0 != b0)) {
  36. auto t0 = a0;
  37. auto [t1, p1] = mnd::pfu::twoSum(a1, b0);
  38. auto [t2, p2] = mnd::pfu::threeTwoSum(b1, c0, p1);
  39. auto [t3, p3] = mnd::pfu::threeTwoSum(c1, d0, p2);
  40. auto [t4, p4] = mnd::pfu::threeTwoSum(d1, e0, p3);
  41. auto [t5, p5] = mnd::pfu::threeTwoSum(e1, f0, p4);
  42. auto [t6, p6] = mnd::pfu::threeTwoSum(f1, g0, p5);
  43. auto t7 = a[7] + b[7] + g1 + p6;
  44. auto[re0, er1] = mnd::pfu::quickTwoSum(t0, t1);
  45. auto[re1, e2] = mnd::pfu::quickTwoSum(er1, t2);
  46. auto[re2, e3] = mnd::pfu::quickTwoSum(e2, t3);
  47. auto[re3, e4] = mnd::pfu::quickTwoSum(e3, t4);
  48. auto[re4, e5] = mnd::pfu::quickTwoSum(e4, t5);
  49. auto[re5, e6] = mnd::pfu::quickTwoSum(e5, t6);
  50. auto[re6, re7] = mnd::pfu::quickTwoSum(e6, t7);
  51. //return { t1, t2, t3, t4, t5, t6 };
  52. return { re0, re1, re2, re3, re4, re5, re6, re7 };
  53. }
  54. auto t0 = a0;
  55. auto [t1, p1] = mnd::pfu::twoSum(a1, b0);
  56. auto [t2, p2, p3] = mnd::pfu::threeSum(b1, c0, p1);
  57. auto [t3, p4, p5, p6] = mnd::pfu::fourSum(c1, d0, p2, p3);
  58. auto [t4, p7, p8] = mnd::pfu::sixThreeSum(d1, e0, p4, p5, p6, 0.0);
  59. auto [t5, p9, p10] = mnd::pfu::sixThreeSum(e1, f0, p7, p8, 0.0, 0.0);
  60. auto [t6, p11] = mnd::pfu::fiveTwoSum(f1, g0, p9, p10, 0.0);
  61. auto t7 = a[7] + b[7] + g1 + p11;
  62. auto[re0, er1] = mnd::pfu::quickTwoSum(t0, t1);
  63. auto[re1, e2] = mnd::pfu::quickTwoSum(er1, t2);
  64. auto[re2, e3] = mnd::pfu::quickTwoSum(e2, t3);
  65. auto[re3, e4] = mnd::pfu::quickTwoSum(e3, t4);
  66. auto[re4, e5] = mnd::pfu::quickTwoSum(e4, t5);
  67. auto[re5, e6] = mnd::pfu::quickTwoSum(e5, t6);
  68. auto[re6, re7] = mnd::pfu::quickTwoSum(e6, t7);
  69. //return { t1, t2, t3, t4, t5, t6 };
  70. return { re0, re1, re2, re3, re4, re5, re6, re7 };
  71. }
  72. inline bool operator>(const mnd::OctaDouble& a, const mnd::OctaDouble& b)
  73. {
  74. if (a[0] > b[0])
  75. return true;
  76. if (a[0] < b[0])
  77. return false;
  78. if (a[1] > b[1])
  79. return true;
  80. if (a[1] < b[1])
  81. return false;
  82. return a[2] > b[2];
  83. }
  84. inline mnd::OctaDouble operator-(const mnd::OctaDouble& a,
  85. const mnd::OctaDouble& b)
  86. {
  87. return a + mnd::OctaDouble{ -b[0], -b[1], -b[2], -b[3], -b[4], -b[5], -b[6], -b[7] };
  88. }
  89. inline mnd::OctaDouble operator*(const mnd::OctaDouble& a,
  90. const mnd::OctaDouble& b)
  91. {
  92. auto[p1_0, p2_0] = mnd::pfu::twoProd(a[0], b[0]);
  93. auto[p2_1, p3_0] = mnd::pfu::twoProd(a[0], b[1]);
  94. auto[p2_2, p3_1] = mnd::pfu::twoProd(a[1], b[0]);
  95. auto[p3_2, p4_0] = mnd::pfu::twoProd(a[2], b[0]);
  96. auto[p3_3, p4_1] = mnd::pfu::twoProd(a[1], b[1]);
  97. auto[p3_4, p4_2] = mnd::pfu::twoProd(a[0], b[2]);
  98. auto[p4_3, p5_0] = mnd::pfu::twoProd(a[3], b[0]);
  99. auto[p4_4, p5_1] = mnd::pfu::twoProd(a[2], b[1]);
  100. auto[p4_5, p5_2] = mnd::pfu::twoProd(a[1], b[2]);
  101. auto[p4_6, p5_3] = mnd::pfu::twoProd(a[0], b[3]);
  102. auto[p5_4, p6_0] = mnd::pfu::twoProd(a[4], b[0]);
  103. auto[p5_5, p6_1] = mnd::pfu::twoProd(a[3], b[1]);
  104. auto[p5_6, p6_2] = mnd::pfu::twoProd(a[2], b[2]);
  105. auto[p5_7, p6_3] = mnd::pfu::twoProd(a[1], b[3]);
  106. auto[p5_8, p6_4] = mnd::pfu::twoProd(a[0], b[4]);
  107. auto[p6_5, p7_0] = mnd::pfu::twoProd(a[5], b[0]);
  108. auto[p6_6, p7_1] = mnd::pfu::twoProd(a[4], b[1]);
  109. auto[p6_7, p7_2] = mnd::pfu::twoProd(a[3], b[2]);
  110. auto[p6_8, p7_3] = mnd::pfu::twoProd(a[2], b[3]);
  111. auto[p6_9, p7_4] = mnd::pfu::twoProd(a[1], b[4]);
  112. auto[p6_10, p7_5] = mnd::pfu::twoProd(a[0], b[5]);
  113. auto[p7_6, p8_0] = mnd::pfu::twoProd(a[6], b[0]);
  114. auto[p7_7, p8_1] = mnd::pfu::twoProd(a[5], b[1]);
  115. auto[p7_8, p8_2] = mnd::pfu::twoProd(a[4], b[2]);
  116. auto[p7_9, p8_3] = mnd::pfu::twoProd(a[3], b[3]);
  117. auto[p7_10, p8_4] = mnd::pfu::twoProd(a[2], b[4]);
  118. auto[p7_11, p8_5] = mnd::pfu::twoProd(a[1], b[5]);
  119. auto[p7_12, p8_6] = mnd::pfu::twoProd(a[0], b[6]);
  120. auto t1 = p1_0;
  121. auto[t2, tl3, tl4] = mnd::pfu::threeSum(p2_0, p2_1, p2_2);
  122. auto[t3, tl4_2, tl5] = mnd::pfu::sixThreeSum(p3_0, p3_1, p3_2, p3_3, p3_4, tl3);
  123. 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);
  124. 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);
  125. auto[t5, tl6_1, tl7] = mnd::pfu::sixThreeSum(x1, x2, x3, tl5, tl5_2, 0.0);
  126. auto[x4, x5, x6] = mnd::pfu::nineThreeSum(p6_0, p6_1, p6_2, p6_3, p6_4, p6_5, p6_6, p6_7, p6_8);
  127. auto[t6, tl7_1, tl8] = mnd::pfu::nineThreeSum(x4, x5, x6, tl6, tl6_1, p6_9, p6_10, 0.0, 0.0);
  128. auto[x7, x8, x9] = mnd::pfu::nineThreeSum(p7_0, p7_1, p7_2, p7_3, p7_4, p7_5, p7_6, p7_7, p7_8);
  129. auto[t7, tl8_1, tl9] = mnd::pfu::nineThreeSum(x7, x8, x9, tl7, tl7_1, p7_9, p7_10, p7_11, p7_12);
  130. auto t8 = tl8 + tl8_1 + tl9 + p8_0 + p8_1 + p8_2 + p8_3 + p8_4 + p8_5 + p8_6 +
  131. a[7] * b[0] + a[6] * b[1] + a[5] * b[2] + a[5] * b[3] +
  132. a[3] * b[4] + a[2] * b[5] + a[1] * b[6] + a[0] * b[7];
  133. auto[re0, e1] = mnd::pfu::quickTwoSum(t1, t2);
  134. auto[re1, e2] = mnd::pfu::quickTwoSum(e1, t3);
  135. auto[re2, e3] = mnd::pfu::quickTwoSum(e2, t4);
  136. auto[re3, e4] = mnd::pfu::quickTwoSum(e3, t5);
  137. auto[re4, e5] = mnd::pfu::quickTwoSum(e4, t6);
  138. auto[re5, e6] = mnd::pfu::quickTwoSum(e5, t7);
  139. auto[re6, re7] = mnd::pfu::quickTwoSum(e6, t8);
  140. return { re0, re1, re2, re3, re4, re5, re6, re7 };
  141. /*double r[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
  142. for (int i = 0; i < 8; i++) {
  143. #pragma GCC unroll(8)
  144. for (int j = 0; i + j < 8; j++) {
  145. auto[s, e] = mnd::pfu::twoProd(a[i], b[j]);
  146. #pragma GCC unroll(8)
  147. for (int k = 0; k < 8; k++) {
  148. auto[x1, x2, x3] = mnd::pfu::threeSum(r[k], s, e);
  149. r[k] = x1;
  150. s = x2;
  151. e = x3;
  152. }
  153. }
  154. }
  155. return mnd::OctaDouble{ r[0], r[1], r[2], r[3], r[4], r[5], r[6], r[7] };
  156. */
  157. }
  158. #endif // MANDEL_OCTADOUBLE_H