CpuGeneratorsAVX.cpp 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511
  1. #include "CpuGenerators.h"
  2. #include <immintrin.h>
  3. #include <omp.h>
  4. #include <cmath>
  5. #include <utility>
  6. #include <memory>
  7. using mnd::CpuGenerator;
  8. namespace mnd
  9. {
  10. template class CpuGenerator<float, mnd::X86_AVX, false>;
  11. template class CpuGenerator<float, mnd::X86_AVX, true>;
  12. template class CpuGenerator<double, mnd::X86_AVX, false>;
  13. template class CpuGenerator<double, mnd::X86_AVX, true>;
  14. template class CpuGenerator<DoubleDouble, mnd::X86_AVX, false>;
  15. template class CpuGenerator<DoubleDouble, mnd::X86_AVX, true>;
  16. }
  17. template<bool parallel>
  18. void CpuGenerator<float, mnd::X86_AVX, parallel>::generate(const mnd::MandelInfo& info, float* data)
  19. {
  20. using T = float;
  21. const MandelViewport& view = info.view;
  22. const float dppf = float(view.width / info.bWidth);
  23. const float viewxf = float(view.x);
  24. __m256 viewx = _mm256_set1_ps(viewxf);
  25. __m256 dpp = _mm256_set1_ps(dppf);
  26. T jX = mnd::convert<T>(info.juliaX);
  27. T jY = mnd::convert<T>(info.juliaY);
  28. __m256 juliaX = { jX, jX, jX, jX, jX, jX, jX, jX };
  29. __m256 juliaY = { jY, jY, jY, jY, jY, jY, jY, jY };
  30. if constexpr(parallel)
  31. omp_set_num_threads(omp_get_num_procs());
  32. #pragma omp parallel for schedule(static, 1) if (parallel)
  33. for (long j = 0; j < info.bHeight; j++) {
  34. T y = T(view.y) + T(j) * T(view.height / info.bHeight);
  35. __m256 ys = _mm256_set1_ps(y);
  36. for (long i = 0; i < info.bWidth; i += 16) {
  37. __m256 pixc = { float(i), float(i + 1), float(i + 2), float(i + 3), float(i + 4), float(i + 5), float(i + 6), float(i + 7) };
  38. __m256 pixc2 = { float(i + 8), float(i + 9), float(i + 10), float(i + 11), float(i + 12), float(i + 13), float(i + 14), float(i + 15) };
  39. __m256 xs = _mm256_add_ps(_mm256_mul_ps(dpp, pixc), viewx);
  40. __m256 xs2 = _mm256_add_ps(_mm256_mul_ps(dpp, pixc2), viewx);
  41. __m256 counter = _mm256_setzero_ps();
  42. __m256 adder = _mm256_set1_ps(1);
  43. __m256 resultsa = _mm256_setzero_ps();
  44. __m256 resultsb = _mm256_setzero_ps();
  45. __m256 counter2 = _mm256_setzero_ps();
  46. __m256 adder2 = _mm256_set1_ps(1);
  47. __m256 resultsa2 = _mm256_setzero_ps();
  48. __m256 resultsb2 = _mm256_setzero_ps();
  49. __m256 threshold = _mm256_set1_ps(16);
  50. __m256 a = xs;
  51. __m256 a2 = xs2;
  52. __m256 b = ys;
  53. __m256 b2 = ys;
  54. __m256 cx = info.julia ? juliaX : xs;
  55. __m256 cx2 = info.julia ? juliaX : xs2;
  56. __m256 cy = info.julia ? juliaY : ys;
  57. if (info.smooth) {
  58. for (int k = 0; k < info.maxIter; k++) {
  59. __m256 aa = _mm256_mul_ps(a, a);
  60. __m256 aa2 = _mm256_mul_ps(a2, a2);
  61. __m256 bb = _mm256_mul_ps(b, b);
  62. __m256 bb2 = _mm256_mul_ps(b2, b2);
  63. __m256 abab = _mm256_mul_ps(a, b); abab = _mm256_add_ps(abab, abab);
  64. __m256 abab2 = _mm256_mul_ps(a2, b2); abab2 = _mm256_add_ps(abab2, abab2);
  65. a = _mm256_add_ps(_mm256_sub_ps(aa, bb), cx);
  66. a2 = _mm256_add_ps(_mm256_sub_ps(aa2, bb2), cx2);
  67. b = _mm256_add_ps(abab, cy);
  68. b2 = _mm256_add_ps(abab2, cy);
  69. __m256 cmp = _mm256_cmp_ps(_mm256_add_ps(aa, bb), threshold, _CMP_LE_OQ);
  70. __m256 cmp2 = _mm256_cmp_ps(_mm256_add_ps(aa2, bb2), threshold, _CMP_LE_OQ);
  71. resultsa = _mm256_or_ps(_mm256_andnot_ps(cmp, resultsa), _mm256_and_ps(cmp, a));
  72. resultsb = _mm256_or_ps(_mm256_andnot_ps(cmp, resultsb), _mm256_and_ps(cmp, b));
  73. resultsa2 = _mm256_or_ps(_mm256_andnot_ps(cmp2, resultsa2), _mm256_and_ps(cmp2, a2));
  74. resultsb2 = _mm256_or_ps(_mm256_andnot_ps(cmp2, resultsb2), _mm256_and_ps(cmp2, b2));
  75. adder = _mm256_and_ps(adder, cmp);
  76. counter = _mm256_add_ps(counter, adder);
  77. adder2 = _mm256_and_ps(adder2, cmp2);
  78. counter2 = _mm256_add_ps(counter2, adder2);
  79. if ((k & 0x7) == 0 && _mm256_testz_ps(cmp, cmp) != 0 && _mm256_testz_ps(cmp2, cmp2) != 0) {
  80. break;
  81. }
  82. }
  83. }
  84. else {
  85. for (int k = 0; k < info.maxIter; k++) {
  86. __m256 aa = _mm256_mul_ps(a, a);
  87. __m256 aa2 = _mm256_mul_ps(a2, a2);
  88. __m256 bb = _mm256_mul_ps(b, b);
  89. __m256 bb2 = _mm256_mul_ps(b2, b2);
  90. __m256 abab = _mm256_mul_ps(a, b); abab = _mm256_add_ps(abab, abab);
  91. __m256 abab2 = _mm256_mul_ps(a2, b2); abab2 = _mm256_add_ps(abab2, abab2);
  92. a = _mm256_add_ps(_mm256_sub_ps(aa, bb), cx);
  93. a2 = _mm256_add_ps(_mm256_sub_ps(aa2, bb2), cx2);
  94. b = _mm256_add_ps(abab, cy);
  95. b2 = _mm256_add_ps(abab2, cy);
  96. __m256 cmp = _mm256_cmp_ps(_mm256_add_ps(aa, bb), threshold, _CMP_LE_OQ);
  97. __m256 cmp2 = _mm256_cmp_ps(_mm256_add_ps(aa2, bb2), threshold, _CMP_LE_OQ);
  98. adder = _mm256_and_ps(adder, cmp);
  99. counter = _mm256_add_ps(counter, adder);
  100. adder2 = _mm256_and_ps(adder2, cmp2);
  101. counter2 = _mm256_add_ps(counter2, adder2);
  102. if ((k & 0x7) == 0 && _mm256_testz_ps(cmp, cmp) != 0 && _mm256_testz_ps(cmp2, cmp2) != 0) {
  103. break;
  104. }
  105. }
  106. }
  107. auto alignVec = [](float* data) -> float* {
  108. void* aligned = data;
  109. ::size_t length = 64;
  110. std::align(32, 8 * sizeof(float), aligned, length);
  111. return static_cast<float*>(aligned);
  112. };
  113. float resData[64];
  114. float* ftRes = alignVec(resData);
  115. float* resa = ftRes + 16;
  116. float* resb = resa + 16;
  117. _mm256_store_ps(ftRes, counter);
  118. _mm256_store_ps(ftRes + 8, counter2);
  119. _mm256_store_ps(resa, resultsa);
  120. _mm256_store_ps(resa + 8, resultsa2);
  121. _mm256_store_ps(resb, resultsb);
  122. _mm256_store_ps(resb + 8, resultsb2);
  123. for (int k = 0; k < 16 && i + k < info.bWidth; k++) {
  124. if (info.smooth) {
  125. data[i + k + j * info.bWidth] = ftRes[k] <= 0 ? info.maxIter :
  126. ftRes[k] >= info.maxIter ? info.maxIter :
  127. ((float)ftRes[k]) + 1 - ::log(::log(resa[k] * resa[k] + resb[k] * resb[k]) / 2) / ::log(2.0f);
  128. }
  129. else {
  130. data[i + k + j * info.bWidth] = ftRes[k] <= 0 ? info.maxIter : ftRes[k];
  131. }
  132. }
  133. }
  134. }
  135. }
  136. template<bool parallel>
  137. void CpuGenerator<double, mnd::X86_AVX, parallel>::generate(const mnd::MandelInfo& info, float* data)
  138. {
  139. using T = double;
  140. const MandelViewport& view = info.view;
  141. const double dppf = double(view.width / info.bWidth);
  142. const double viewxf = double(view.x);
  143. __m256d viewx = { viewxf, viewxf, viewxf, viewxf };
  144. __m256d dpp = { dppf, dppf, dppf, dppf };
  145. T jX = mnd::convert<T>(info.juliaX);
  146. T jY = mnd::convert<T>(info.juliaY);
  147. __m256d juliaX = { jX, jX, jX, jX };
  148. __m256d juliaY = { jY, jY, jY, jY };
  149. if constexpr(parallel)
  150. omp_set_num_threads(omp_get_num_procs());
  151. #pragma omp parallel for schedule(static, 1) if (parallel)
  152. for (long j = 0; j < info.bHeight; j++) {
  153. T y = T(view.y + T(j) * view.height / info.bHeight);
  154. __m256d ys = { y, y, y, y };
  155. for (long i = 0; i < info.bWidth; i += 8) {
  156. __m256d pixc = { double(i), double(i + 1), double(i + 2), double(i + 3) };
  157. __m256d pixc2 = { double(i + 4), double(i + 5), double(i + 6), double(i + 7) };
  158. __m256d xs = _mm256_add_pd(_mm256_mul_pd(dpp, pixc), viewx);
  159. __m256d xs2 = _mm256_add_pd(_mm256_mul_pd(dpp, pixc2), viewx);
  160. int itRes[4] = { 0, 0, 0, 0 };
  161. __m256d threshold = { 16.0, 16.0, 16.0, 16.0 };
  162. __m256d counter = { 0, 0, 0, 0 };
  163. __m256d adder = { 1, 1, 1, 1 };
  164. __m256d counter2 = { 0, 0, 0, 0 };
  165. __m256d adder2 = { 1, 1, 1, 1 };
  166. __m256d resultsa = { 0, 0, 0, 0 };
  167. __m256d resultsb = { 0, 0, 0, 0 };
  168. __m256d resultsa2 = { 0, 0, 0, 0 };
  169. __m256d resultsb2 = { 0, 0, 0, 0 };
  170. __m256d a = xs;
  171. __m256d b = ys;
  172. __m256d a2 = xs2;
  173. __m256d b2 = ys;
  174. __m256d cx = info.julia ? juliaX : xs;
  175. __m256d cx2 = info.julia ? juliaX : xs2;
  176. __m256d cy = info.julia ? juliaY : ys;
  177. if (info.smooth) {
  178. for (int k = 0; k < info.maxIter; k++) {
  179. __m256d aa = _mm256_mul_pd(a, a);
  180. __m256d aa2 = _mm256_mul_pd(a2, a2);
  181. __m256d bb = _mm256_mul_pd(b, b);
  182. __m256d bb2 = _mm256_mul_pd(b2, b2);
  183. __m256d abab = _mm256_mul_pd(a, b); abab = _mm256_add_pd(abab, abab);
  184. __m256d abab2 = _mm256_mul_pd(a2, b2); abab2 = _mm256_add_pd(abab2, abab2);
  185. a = _mm256_add_pd(_mm256_sub_pd(aa, bb), cx);
  186. a2 = _mm256_add_pd(_mm256_sub_pd(aa2, bb2), cx2);
  187. b = _mm256_add_pd(abab, cy);
  188. b2 = _mm256_add_pd(abab2, cy);
  189. __m256d cmp = _mm256_cmp_pd(_mm256_add_pd(aa, bb), threshold, _CMP_LE_OQ);
  190. __m256d cmp2 = _mm256_cmp_pd(_mm256_add_pd(aa2, bb2), threshold, _CMP_LE_OQ);
  191. resultsa = _mm256_or_pd(_mm256_andnot_pd(cmp, resultsa), _mm256_and_pd(cmp, a));
  192. resultsb = _mm256_or_pd(_mm256_andnot_pd(cmp, resultsb), _mm256_and_pd(cmp, b));
  193. resultsa2 = _mm256_or_pd(_mm256_andnot_pd(cmp2, resultsa2), _mm256_and_pd(cmp2, a2));
  194. resultsb2 = _mm256_or_pd(_mm256_andnot_pd(cmp2, resultsb2), _mm256_and_pd(cmp2, b2));
  195. adder = _mm256_and_pd(adder, cmp);
  196. counter = _mm256_add_pd(counter, adder);
  197. adder2 = _mm256_and_pd(adder2, cmp2);
  198. counter2 = _mm256_add_pd(counter2, adder2);
  199. if ((k & 0x7) == 0 && _mm256_testz_si256(_mm256_castpd_si256(cmp), _mm256_castpd_si256(cmp)) != 0 &&
  200. _mm256_testz_si256(_mm256_castpd_si256(cmp2), _mm256_castpd_si256(cmp2)) != 0) {
  201. break;
  202. }
  203. }
  204. }
  205. else {
  206. for (int k = 0; k < info.maxIter; k++) {
  207. __m256d aa = _mm256_mul_pd(a, a);
  208. __m256d aa2 = _mm256_mul_pd(a2, a2);
  209. __m256d bb = _mm256_mul_pd(b, b);
  210. __m256d bb2 = _mm256_mul_pd(b2, b2);
  211. __m256d abab = _mm256_mul_pd(a, b); abab = _mm256_add_pd(abab, abab);
  212. __m256d abab2 = _mm256_mul_pd(a2, b2); abab2 = _mm256_add_pd(abab2, abab2);
  213. a = _mm256_add_pd(_mm256_sub_pd(aa, bb), cx);
  214. a2 = _mm256_add_pd(_mm256_sub_pd(aa2, bb2), cx2);
  215. b = _mm256_add_pd(abab, cy);
  216. b2 = _mm256_add_pd(abab2, cy);
  217. __m256d cmp = _mm256_cmp_pd(_mm256_add_pd(aa, bb), threshold, _CMP_LE_OQ);
  218. __m256d cmp2 = _mm256_cmp_pd(_mm256_add_pd(aa2, bb2), threshold, _CMP_LE_OQ);
  219. adder = _mm256_and_pd(adder, cmp);
  220. counter = _mm256_add_pd(counter, adder);
  221. adder2 = _mm256_and_pd(adder2, cmp2);
  222. counter2 = _mm256_add_pd(counter2, adder2);
  223. if ((k & 0x7) == 0 && _mm256_testz_si256(_mm256_castpd_si256(cmp), _mm256_castpd_si256(cmp)) != 0 &&
  224. _mm256_testz_si256(_mm256_castpd_si256(cmp2), _mm256_castpd_si256(cmp2)) != 0) {
  225. break;
  226. }
  227. }
  228. }
  229. auto alignVec = [](double* data) -> double* {
  230. void* aligned = data;
  231. ::size_t length = 64;
  232. std::align(32, 4 * sizeof(double), aligned, length);
  233. return static_cast<double*>(aligned);
  234. };
  235. double resData[8];
  236. double* ftRes = alignVec(resData);
  237. double* resa = (double*) &resultsa;
  238. double* resb = (double*) &resultsb;
  239. _mm256_store_pd(ftRes, counter);
  240. for (int k = 0; k < 4 && i + k < info.bWidth; k++) {
  241. if (info.smooth)
  242. data[i + k + j * info.bWidth] = ftRes[k] <= 0 ? float(info.maxIter) :
  243. ftRes[k] >= info.maxIter ? float(info.maxIter) :
  244. float(((float)ftRes[k]) + 1 - ::log(::log(resa[k] * resa[k] + resb[k] * resb[k]) / 2) / ::logf(2.0f));
  245. else
  246. data[i + k + j * info.bWidth] = ftRes[k] > 0 ? float(ftRes[k]) : info.maxIter;
  247. }
  248. resa = (double*) &resultsa2;
  249. resb = (double*) &resultsb2;
  250. _mm256_store_pd(ftRes, counter2);
  251. i += 4;
  252. for (int k = 0; k < 4 && i + k < info.bWidth; k++) {
  253. if (info.smooth)
  254. data[i + k + j * info.bWidth] = ftRes[k] <= 0 ? float(info.maxIter) :
  255. ftRes[k] >= info.maxIter ? float(info.maxIter) :
  256. float(((float)ftRes[k]) + 1 - ::log(::log(resa[k] * resa[k] + resb[k] * resb[k]) / 2) / ::logf(2.0f));
  257. else
  258. data[i + k + j * info.bWidth] = ftRes[k] > 0 ? float(ftRes[k]) : info.maxIter;
  259. }
  260. i -= 4;
  261. }
  262. }
  263. }
  264. struct VecPair
  265. {
  266. __m256d a;
  267. __m256d b;
  268. };
  269. static inline VecPair quickTwoSum(__m256d a, __m256d b)
  270. {
  271. __m256d s = _mm256_add_pd(a, b);
  272. __m256d e = _mm256_sub_pd(b, _mm256_sub_pd(s, a));
  273. return { s, e };
  274. }
  275. static inline VecPair quickTwoDiff(__m256d a, __m256d b)
  276. {
  277. __m256d s = _mm256_sub_pd(a, b);
  278. __m256d e = _mm256_sub_pd(_mm256_sub_pd(a, s), b);
  279. return { s, e };
  280. }
  281. static inline VecPair twoSum(__m256d a, __m256d b)
  282. {
  283. __m256d s = _mm256_add_pd(a, b);
  284. __m256d bb = _mm256_sub_pd(s, a);
  285. __m256d e = _mm256_add_pd(_mm256_sub_pd(a, _mm256_sub_pd(s, bb)), _mm256_sub_pd(b, bb));
  286. return { s, e };
  287. }
  288. static inline VecPair twoDiff(__m256d a, __m256d b)
  289. {
  290. __m256d s = _mm256_sub_pd(a, b);
  291. __m256d bb = _mm256_sub_pd(s, a);
  292. __m256d e = _mm256_sub_pd(_mm256_sub_pd(a, _mm256_sub_pd(s, bb)), _mm256_add_pd(b, bb));
  293. return { s, e };
  294. }
  295. static inline VecPair split(__m256d a)
  296. {
  297. /*
  298. // -- this should never happen when doing mandelbrot calculations,
  299. // so we omit this check.
  300. if (a > _QD_SPLIT_THRESH || a < -_QD_SPLIT_THRESH) {
  301. a *= 3.7252902984619140625e-09; // 2^-28
  302. temp = _QD_SPLITTER * a;
  303. hi = temp - (temp - a);
  304. lo = a - hi;
  305. hi *= 268435456.0; // 2^28
  306. lo *= 268435456.0; // 2^28
  307. } else {
  308. temp = _QD_SPLITTER * a;
  309. hi = temp - (temp - a);
  310. lo = a - hi;
  311. }
  312. */
  313. static const __m256d SPLITTER = { 134217729.0, 134217729.0, 134217729.0, 134217729.0 };
  314. __m256d temp = _mm256_mul_pd(SPLITTER, a);
  315. __m256d hi = _mm256_sub_pd(temp, _mm256_sub_pd(temp, a));
  316. __m256d lo = _mm256_sub_pd(a, hi);
  317. return { hi, lo };
  318. }
  319. static inline VecPair twoProd(__m256d a, __m256d b)
  320. {
  321. __m256d p = _mm256_mul_pd(a, b);
  322. auto[a_hi, a_lo] = split(a);
  323. auto[b_hi, b_lo] = split(b);
  324. __m256d err = _mm256_add_pd(_mm256_add_pd(_mm256_sub_pd(_mm256_mul_pd(a_hi, b_hi), p), _mm256_add_pd(_mm256_mul_pd(a_hi, b_lo), _mm256_mul_pd(a_lo, b_hi))), _mm256_mul_pd(a_lo, b_lo));
  325. return { p, err };
  326. }
  327. struct AvxDoubleDouble
  328. {
  329. __m256d x[2];
  330. inline AvxDoubleDouble(__m256d a, __m256d b) :
  331. x{ a, b }
  332. {}
  333. inline AvxDoubleDouble operator + (const AvxDoubleDouble& sm) const
  334. {
  335. auto[s, e] = twoSum(x[0], sm.x[0]);
  336. e = _mm256_add_pd(e, _mm256_add_pd(x[1], sm.x[1]));
  337. auto[r1, r2] = quickTwoSum(s, e);
  338. return AvxDoubleDouble{ r1, r2 };
  339. }
  340. inline AvxDoubleDouble operator - (const AvxDoubleDouble& sm) const
  341. {
  342. auto[s, e] = twoDiff(x[0], sm.x[0]);
  343. e = _mm256_add_pd(e, x[1]);
  344. e = _mm256_sub_pd(e, sm.x[1]);
  345. auto[r1, r2] = quickTwoSum(s, e);
  346. return AvxDoubleDouble{ r1, r2 };
  347. }
  348. inline AvxDoubleDouble operator * (const AvxDoubleDouble& sm) const
  349. {
  350. auto[p1, p2] = twoProd(this->x[0], sm.x[0]);
  351. p2 = _mm256_add_pd(p2,
  352. _mm256_add_pd(_mm256_mul_pd(sm.x[1], x[0]), _mm256_mul_pd(sm.x[0], x[1])) );
  353. auto[r1, r2] = quickTwoSum(p1, p2);
  354. return AvxDoubleDouble{ r1, r2 };
  355. }
  356. };
  357. template<bool parallel>
  358. void CpuGenerator<mnd::DoubleDouble, mnd::X86_AVX, parallel>::generate(const mnd::MandelInfo& info, float* data)
  359. {
  360. const MandelViewport& view = info.view;
  361. using T = DoubleDouble;
  362. T viewx = mnd::convert<T>(view.x);
  363. T viewy = mnd::convert<T>(view.y);
  364. T wpp = mnd::convert<T>(view.width / info.bWidth);
  365. T hpp = mnd::convert<T>(view.height / info.bHeight);
  366. T jX = mnd::convert<T>(info.juliaX);
  367. T jY = mnd::convert<T>(info.juliaY);
  368. AvxDoubleDouble juliaX = { __m256d{ jX.x[0], jX.x[0], jX.x[0], jX.x[0] }, __m256d{ jX.x[1], jX.x[1], jX.x[1], jX.x[1] } };
  369. AvxDoubleDouble juliaY = { __m256d{ jY.x[0], jY.x[0], jY.x[0], jY.x[0] }, __m256d{ jY.x[1], jY.x[1], jY.x[1], jY.x[1] } };
  370. if constexpr(parallel)
  371. omp_set_num_threads(omp_get_num_procs());
  372. #pragma omp parallel for schedule(static, 1) if (parallel)
  373. for (long j = 0; j < info.bHeight; j++) {
  374. T y = viewy + T(double(j)) * hpp;
  375. __m256d y0s = { y.x[0], y.x[0], y.x[0], y.x[0] };
  376. __m256d y1s = { y.x[1], y.x[1], y.x[1], y.x[1] };
  377. AvxDoubleDouble ys{ y0s, y1s };
  378. for (long i = 0; i < info.bWidth; i += 4) {
  379. T x1 = viewx + T(double(i)) * wpp;
  380. T x2 = x1 + wpp;
  381. T x3 = x2 + wpp;
  382. T x4 = x3 + wpp;
  383. __m256d x0s = {
  384. x1.x[0], x2.x[0], x3.x[0], x4.x[0],
  385. };
  386. __m256d x1s = {
  387. x1.x[1], x2.x[1], x3.x[1], x4.x[1],
  388. };
  389. AvxDoubleDouble xs{ x0s, x1s };
  390. AvxDoubleDouble cx = info.julia ? juliaX : xs;
  391. AvxDoubleDouble cy = info.julia ? juliaY : ys;
  392. int itRes[4] = { 0, 0, 0, 0 };
  393. __m256d threshold = { 16.0, 16.0, 16.0, 16.0 };
  394. __m256d counter = { 0, 0, 0, 0 };
  395. __m256d adder = { 1, 1, 1, 1 };
  396. AvxDoubleDouble a = xs;
  397. AvxDoubleDouble b = ys;
  398. __m256d resultsa;
  399. __m256d resultsb;
  400. for (int k = 0; k < info.maxIter; k++) {
  401. AvxDoubleDouble aa = a * a;
  402. AvxDoubleDouble bb = b * b;
  403. AvxDoubleDouble abab = a * b; abab = abab + abab;
  404. a = aa - bb + cx;
  405. b = abab + cy;
  406. __m256d cmp = _mm256_cmp_pd(_mm256_add_pd(aa.x[0], bb.x[0]), threshold, _CMP_LE_OQ);
  407. if (info.smooth) {
  408. resultsa = _mm256_or_pd(_mm256_andnot_pd(cmp, resultsa), _mm256_and_pd(cmp, a.x[0]));
  409. resultsb = _mm256_or_pd(_mm256_andnot_pd(cmp, resultsb), _mm256_and_pd(cmp, b.x[0]));
  410. }
  411. adder = _mm256_and_pd(adder, cmp);
  412. counter = _mm256_add_pd(counter, adder);
  413. if (_mm256_testz_si256(_mm256_castpd_si256(cmp), _mm256_castpd_si256(cmp)) != 0) {
  414. break;
  415. }
  416. }
  417. auto alignVec = [](double* data) -> double* {
  418. void* aligned = data;
  419. ::size_t length = 64;
  420. std::align(32, 4 * sizeof(double), aligned, length);
  421. return static_cast<double*>(aligned);
  422. };
  423. double resData[8];
  424. double* ftRes = alignVec(resData);
  425. double* resa = (double*) &resultsa;
  426. double* resb = (double*) &resultsb;
  427. _mm256_store_pd(ftRes, counter);
  428. for (int k = 0; k < 4 && i + k < info.bWidth; k++) {
  429. if (info.smooth)
  430. data[i + k + j * info.bWidth] = float(ftRes[k] <= 0 ? info.maxIter :
  431. ftRes[k] >= info.maxIter ? info.maxIter :
  432. ((float)ftRes[k]) + 1 - ::logf(::logf(float(resa[k] * resa[k] + resb[k] * resb[k])) / 2) / ::logf(2.0f));
  433. else
  434. data[i + k + j * info.bWidth] = ftRes[k] > 0 ? float(ftRes[k]) : info.maxIter;
  435. }
  436. }
  437. }
  438. }