Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: SSE2 helper
5 : * Author: Even Rouault <even dot rouault at spatialys dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2014, Even Rouault <even dot rouault at spatialys dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #ifndef GDALSSE_PRIV_H_INCLUDED
14 : #define GDALSSE_PRIV_H_INCLUDED
15 :
16 : #ifndef DOXYGEN_SKIP
17 :
18 : #include "cpl_port.h"
19 :
20 : /* We restrict to 64bit processors because they are guaranteed to have SSE2 */
21 : /* Could possibly be used too on 32bit, but we would need to check at runtime */
22 : #if (defined(__x86_64) || defined(_M_X64) || defined(USE_SSE2)) && \
23 : !defined(USE_SSE2_EMULATION)
24 :
25 : #include <string.h>
26 :
27 : #ifdef USE_NEON_OPTIMIZATIONS
28 : #include "include_sse2neon.h"
29 : #else
30 : /* Requires SSE2 */
31 : #include <emmintrin.h>
32 :
33 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
34 : #include <smmintrin.h>
35 : #endif
36 : #endif
37 :
38 : #include "gdal_priv_templates.hpp"
39 :
40 534518 : static inline __m128i GDALCopyInt16ToXMM(const void *ptr)
41 : {
42 : unsigned short s;
43 534518 : memcpy(&s, ptr, 2);
44 1069040 : return _mm_cvtsi32_si128(s);
45 : }
46 :
47 418974635 : static inline __m128i GDALCopyInt32ToXMM(const void *ptr)
48 : {
49 : GInt32 i;
50 418974635 : memcpy(&i, ptr, 4);
51 837950280 : return _mm_cvtsi32_si128(i);
52 : }
53 :
54 272220 : static inline __m128i GDALCopyInt64ToXMM(const void *ptr)
55 : {
56 : #if defined(__i386__) || defined(_M_IX86)
57 : return _mm_loadl_epi64(static_cast<const __m128i *>(ptr));
58 : #else
59 : GInt64 i;
60 272220 : memcpy(&i, ptr, 8);
61 544440 : return _mm_cvtsi64_si128(i);
62 : #endif
63 : }
64 :
65 : #ifndef GDALCopyXMMToInt16_defined
66 : #define GDALCopyXMMToInt16_defined
67 :
68 : static inline void GDALCopyXMMToInt16(const __m128i xmm, void *pDest)
69 : {
70 : GInt16 i = static_cast<GInt16>(_mm_extract_epi16(xmm, 0));
71 : memcpy(pDest, &i, 2);
72 : }
73 : #endif
74 :
75 : class XMMReg4Int;
76 :
77 : class XMMReg4Double;
78 :
79 : class XMMReg4Float
80 : {
81 : public:
82 : __m128 xmm;
83 :
84 26427500 : XMMReg4Float()
85 : #if !defined(_MSC_VER)
86 26427500 : : xmm(_mm_undefined_ps())
87 : #endif
88 : {
89 26427500 : }
90 :
91 35391 : XMMReg4Float(const XMMReg4Float &other) : xmm(other.xmm)
92 : {
93 35391 : }
94 :
95 695 : static inline XMMReg4Float Zero()
96 : {
97 695 : XMMReg4Float reg;
98 695 : reg.Zeroize();
99 695 : return reg;
100 : }
101 :
102 281530 : static inline XMMReg4Float Set1(float f)
103 : {
104 281530 : XMMReg4Float reg;
105 281530 : reg.xmm = _mm_set1_ps(f);
106 281530 : return reg;
107 : }
108 :
109 86824 : static inline XMMReg4Float LoadAllVal(const float *ptr)
110 : {
111 86824 : return Load4Val(ptr);
112 : }
113 :
114 86824 : static inline XMMReg4Float Load4Val(const float *ptr)
115 : {
116 86824 : XMMReg4Float reg;
117 86824 : reg.nsLoad4Val(ptr);
118 86824 : return reg;
119 : }
120 :
121 1616420 : static inline XMMReg4Float Load4Val(const unsigned char *ptr)
122 : {
123 1616420 : XMMReg4Float reg;
124 1616420 : reg.nsLoad4Val(ptr);
125 1616420 : return reg;
126 : }
127 :
128 : static inline XMMReg4Float Load4Val(const short *ptr)
129 : {
130 : XMMReg4Float reg;
131 : reg.nsLoad4Val(ptr);
132 : return reg;
133 : }
134 :
135 : static inline XMMReg4Float Load4Val(const unsigned short *ptr)
136 : {
137 : XMMReg4Float reg;
138 : reg.nsLoad4Val(ptr);
139 : return reg;
140 : }
141 :
142 : static inline XMMReg4Float Load4Val(const int *ptr)
143 : {
144 : XMMReg4Float reg;
145 : reg.nsLoad4Val(ptr);
146 : return reg;
147 : }
148 :
149 1616420 : static inline XMMReg4Float Equals(const XMMReg4Float &expr1,
150 : const XMMReg4Float &expr2)
151 : {
152 1616420 : XMMReg4Float reg;
153 1616420 : reg.xmm = _mm_cmpeq_ps(expr1.xmm, expr2.xmm);
154 1616420 : return reg;
155 : }
156 :
157 : static inline XMMReg4Float NotEquals(const XMMReg4Float &expr1,
158 : const XMMReg4Float &expr2)
159 : {
160 : XMMReg4Float reg;
161 : reg.xmm = _mm_cmpneq_ps(expr1.xmm, expr2.xmm);
162 : return reg;
163 : }
164 :
165 538808 : static inline XMMReg4Float Lesser(const XMMReg4Float &expr1,
166 : const XMMReg4Float &expr2)
167 : {
168 538808 : XMMReg4Float reg;
169 538808 : reg.xmm = _mm_cmplt_ps(expr1.xmm, expr2.xmm);
170 538808 : return reg;
171 : }
172 :
173 : static inline XMMReg4Float Greater(const XMMReg4Float &expr1,
174 : const XMMReg4Float &expr2)
175 : {
176 : XMMReg4Float reg;
177 : reg.xmm = _mm_cmpgt_ps(expr1.xmm, expr2.xmm);
178 : return reg;
179 : }
180 :
181 : static inline XMMReg4Float And(const XMMReg4Float &expr1,
182 : const XMMReg4Float &expr2)
183 : {
184 : XMMReg4Float reg;
185 : reg.xmm = _mm_and_ps(expr1.xmm, expr2.xmm);
186 : return reg;
187 : }
188 :
189 2155230 : static inline XMMReg4Float Ternary(const XMMReg4Float &cond,
190 : const XMMReg4Float &true_expr,
191 : const XMMReg4Float &false_expr)
192 : {
193 2155230 : XMMReg4Float reg;
194 4310460 : reg.xmm = _mm_or_ps(_mm_and_ps(cond.xmm, true_expr.xmm),
195 2155230 : _mm_andnot_ps(cond.xmm, false_expr.xmm));
196 2155230 : return reg;
197 : }
198 :
199 1077620 : static inline XMMReg4Float Min(const XMMReg4Float &expr1,
200 : const XMMReg4Float &expr2)
201 : {
202 1077620 : XMMReg4Float reg;
203 1077620 : reg.xmm = _mm_min_ps(expr1.xmm, expr2.xmm);
204 1077620 : return reg;
205 : }
206 :
207 1663420 : static inline XMMReg4Float Max(const XMMReg4Float &expr1,
208 : const XMMReg4Float &expr2)
209 : {
210 1663420 : XMMReg4Float reg;
211 1663420 : reg.xmm = _mm_max_ps(expr1.xmm, expr2.xmm);
212 1663420 : return reg;
213 : }
214 :
215 88312 : inline void nsLoad4Val(const float *ptr)
216 : {
217 88312 : xmm = _mm_loadu_ps(ptr);
218 88312 : }
219 :
220 372 : static inline void Load16Val(const float *ptr, XMMReg4Float &r0,
221 : XMMReg4Float &r1, XMMReg4Float &r2,
222 : XMMReg4Float &r3)
223 : {
224 372 : r0.nsLoad4Val(ptr);
225 372 : r1.nsLoad4Val(ptr + 4);
226 372 : r2.nsLoad4Val(ptr + 8);
227 372 : r3.nsLoad4Val(ptr + 12);
228 372 : }
229 :
230 1024 : inline void nsLoad4Val(const int *ptr)
231 : {
232 : const __m128i xmm_i =
233 1024 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr));
234 1024 : xmm = _mm_cvtepi32_ps(xmm_i);
235 1024 : }
236 :
237 256 : static inline void Load16Val(const int *ptr, XMMReg4Float &r0,
238 : XMMReg4Float &r1, XMMReg4Float &r2,
239 : XMMReg4Float &r3)
240 : {
241 256 : r0.nsLoad4Val(ptr);
242 256 : r1.nsLoad4Val(ptr + 4);
243 256 : r2.nsLoad4Val(ptr + 8);
244 256 : r3.nsLoad4Val(ptr + 12);
245 256 : }
246 :
247 2156260 : static inline __m128i cvtepu8_epi32(__m128i x)
248 : {
249 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
250 : return _mm_cvtepu8_epi32(x);
251 : #else
252 6468770 : return _mm_unpacklo_epi16(_mm_unpacklo_epi8(x, _mm_setzero_si128()),
253 2156260 : _mm_setzero_si128());
254 : #endif
255 : }
256 :
257 1616420 : inline void nsLoad4Val(const unsigned char *ptr)
258 : {
259 1616420 : const __m128i xmm_i = GDALCopyInt32ToXMM(ptr);
260 1616420 : xmm = _mm_cvtepi32_ps(cvtepu8_epi32(xmm_i));
261 1616420 : }
262 :
263 269404 : static inline void Load8Val(const unsigned char *ptr, XMMReg4Float &r0,
264 : XMMReg4Float &r1)
265 : {
266 269404 : const __m128i xmm_i = GDALCopyInt64ToXMM(ptr);
267 269404 : r0.xmm = _mm_cvtepi32_ps(cvtepu8_epi32(xmm_i));
268 269404 : r1.xmm = _mm_cvtepi32_ps(cvtepu8_epi32(_mm_srli_si128(xmm_i, 4)));
269 269404 : }
270 :
271 256 : static inline void Load16Val(const unsigned char *ptr, XMMReg4Float &r0,
272 : XMMReg4Float &r1, XMMReg4Float &r2,
273 : XMMReg4Float &r3)
274 : {
275 : const __m128i xmm_i =
276 256 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr));
277 256 : r0.xmm = _mm_cvtepi32_ps(cvtepu8_epi32(xmm_i));
278 256 : r1.xmm = _mm_cvtepi32_ps(cvtepu8_epi32(_mm_srli_si128(xmm_i, 4)));
279 256 : r2.xmm = _mm_cvtepi32_ps(cvtepu8_epi32(_mm_srli_si128(xmm_i, 8)));
280 256 : r3.xmm = _mm_cvtepi32_ps(cvtepu8_epi32(_mm_srli_si128(xmm_i, 12)));
281 256 : }
282 :
283 1024 : static inline __m128i cvtepi16_epi32(__m128i x)
284 : {
285 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
286 : return _mm_cvtepi16_epi32(x);
287 : #else
288 : /* 0|0|0|0|b|b|a|a --> 0|0|0|0|sign(b)|b|sign(a)|a */
289 1024 : return _mm_srai_epi32(
290 : /* 0|0|0|0|0|0|b|a --> 0|0|0|0|b|b|a|a */
291 1024 : _mm_unpacklo_epi16(x, x), 16);
292 : #endif
293 : }
294 :
295 : inline void nsLoad4Val(const short *ptr)
296 : {
297 : const __m128i xmm_i = GDALCopyInt64ToXMM(ptr);
298 : xmm = _mm_cvtepi32_ps(cvtepi16_epi32(xmm_i));
299 : }
300 :
301 512 : static inline void Load8Val(const short *ptr, XMMReg4Float &r0,
302 : XMMReg4Float &r1)
303 : {
304 : const __m128i xmm_i =
305 512 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr));
306 512 : r0.xmm = _mm_cvtepi32_ps(cvtepi16_epi32(xmm_i));
307 512 : r1.xmm = _mm_cvtepi32_ps(cvtepi16_epi32(_mm_srli_si128(xmm_i, 8)));
308 512 : }
309 :
310 256 : static inline void Load16Val(const short *ptr, XMMReg4Float &r0,
311 : XMMReg4Float &r1, XMMReg4Float &r2,
312 : XMMReg4Float &r3)
313 : {
314 256 : Load8Val(ptr, r0, r1);
315 256 : Load8Val(ptr + 8, r2, r3);
316 256 : }
317 :
318 1024 : static inline __m128i cvtepu16_epi32(__m128i x)
319 : {
320 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
321 : return _mm_cvtepu16_epi32(x);
322 : #else
323 1024 : return _mm_unpacklo_epi16(
324 1024 : x, _mm_setzero_si128()); /* 0|0|0|0|0|0|b|a --> 0|0|0|0|0|b|0|a */
325 : #endif
326 : }
327 :
328 : inline void nsLoad4Val(const unsigned short *ptr)
329 : {
330 : const __m128i xmm_i = GDALCopyInt64ToXMM(ptr);
331 : xmm = _mm_cvtepi32_ps(cvtepu16_epi32(xmm_i));
332 : }
333 :
334 512 : static inline void Load8Val(const unsigned short *ptr, XMMReg4Float &r0,
335 : XMMReg4Float &r1)
336 : {
337 : const __m128i xmm_i =
338 512 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr));
339 512 : r0.xmm = _mm_cvtepi32_ps(cvtepu16_epi32(xmm_i));
340 512 : r1.xmm = _mm_cvtepi32_ps(cvtepu16_epi32(_mm_srli_si128(xmm_i, 8)));
341 512 : }
342 :
343 256 : static inline void Load16Val(const unsigned short *ptr, XMMReg4Float &r0,
344 : XMMReg4Float &r1, XMMReg4Float &r2,
345 : XMMReg4Float &r3)
346 : {
347 256 : Load8Val(ptr, r0, r1);
348 256 : Load8Val(ptr + 8, r2, r3);
349 256 : }
350 :
351 695 : inline void Zeroize()
352 : {
353 695 : xmm = _mm_setzero_ps();
354 695 : }
355 :
356 1077620 : inline XMMReg4Float &operator=(const XMMReg4Float &other)
357 : {
358 1077620 : xmm = other.xmm;
359 1077620 : return *this;
360 : }
361 :
362 59849 : inline XMMReg4Float &operator+=(const XMMReg4Float &other)
363 : {
364 59849 : xmm = _mm_add_ps(xmm, other.xmm);
365 59849 : return *this;
366 : }
367 :
368 10853 : inline XMMReg4Float &operator-=(const XMMReg4Float &other)
369 : {
370 10853 : xmm = _mm_sub_ps(xmm, other.xmm);
371 10853 : return *this;
372 : }
373 :
374 : inline XMMReg4Float &operator*=(const XMMReg4Float &other)
375 : {
376 : xmm = _mm_mul_ps(xmm, other.xmm);
377 : return *this;
378 : }
379 :
380 3470160 : inline XMMReg4Float operator+(const XMMReg4Float &other) const
381 : {
382 3470160 : XMMReg4Float ret;
383 3470160 : ret.xmm = _mm_add_ps(xmm, other.xmm);
384 3470160 : return ret;
385 : }
386 :
387 4892680 : inline XMMReg4Float operator-(const XMMReg4Float &other) const
388 : {
389 4892680 : XMMReg4Float ret;
390 4892680 : ret.xmm = _mm_sub_ps(xmm, other.xmm);
391 4892680 : return ret;
392 : }
393 :
394 5670030 : inline XMMReg4Float operator*(const XMMReg4Float &other) const
395 : {
396 5670030 : XMMReg4Float ret;
397 5670030 : ret.xmm = _mm_mul_ps(xmm, other.xmm);
398 5670030 : return ret;
399 : }
400 :
401 538808 : inline XMMReg4Float operator/(const XMMReg4Float &other) const
402 : {
403 538808 : XMMReg4Float ret;
404 538808 : ret.xmm = _mm_div_ps(xmm, other.xmm);
405 538808 : return ret;
406 : }
407 :
408 538808 : inline XMMReg4Float inverse() const
409 : {
410 538808 : XMMReg4Float ret;
411 1077620 : ret.xmm = _mm_div_ps(_mm_set1_ps(1.0f), xmm);
412 538808 : return ret;
413 : }
414 :
415 : inline XMMReg4Int truncate_to_int() const;
416 :
417 21706 : inline XMMReg4Float cast_to_float() const
418 : {
419 21706 : return *this;
420 : }
421 :
422 : inline XMMReg4Double cast_to_double() const;
423 :
424 46991 : inline XMMReg4Float approx_inv_sqrt(const XMMReg4Float &one,
425 : const XMMReg4Float &half) const
426 : {
427 46991 : __m128 reg = xmm;
428 93982 : __m128 reg_half = _mm_mul_ps(reg, half.xmm);
429 : // Compute rough approximation of 1 / sqrt(b) with _mm_rsqrt_ps
430 46991 : reg = _mm_rsqrt_ps(reg);
431 : // And perform one step of Newton-Raphson approximation to improve it
432 : // approx_inv_sqrt_x = approx_inv_sqrt_x*(1.5 -
433 : // 0.5*x*approx_inv_sqrt_x*approx_inv_sqrt_x);
434 93982 : const __m128 one_and_a_half = _mm_add_ps(one.xmm, half.xmm);
435 140973 : reg = _mm_mul_ps(
436 : reg, _mm_sub_ps(one_and_a_half,
437 : _mm_mul_ps(reg_half, _mm_mul_ps(reg, reg))));
438 46991 : XMMReg4Float ret;
439 46991 : ret.xmm = reg;
440 46991 : return ret;
441 : }
442 :
443 46991 : inline void StoreAllVal(float *ptr) const
444 : {
445 46991 : Store4Val(ptr);
446 46991 : }
447 :
448 49823 : inline void Store4Val(float *ptr) const
449 : {
450 49823 : _mm_storeu_ps(ptr, xmm);
451 49823 : }
452 :
453 : inline void Store4ValAligned(float *ptr) const
454 : {
455 : _mm_store_ps(ptr, xmm);
456 : }
457 :
458 : inline operator float() const
459 : {
460 : return _mm_cvtss_f32(xmm);
461 : }
462 : };
463 :
464 : class XMMReg4Int
465 : {
466 : public:
467 : __m128i xmm;
468 :
469 3127700 : XMMReg4Int()
470 : #if !defined(_MSC_VER)
471 3127700 : : xmm(_mm_undefined_si128())
472 : #endif
473 : {
474 3127700 : }
475 :
476 36138 : XMMReg4Int(const XMMReg4Int &other) : xmm(other.xmm)
477 : {
478 36138 : }
479 :
480 : inline XMMReg4Int &operator=(const XMMReg4Int &other)
481 : {
482 : xmm = other.xmm;
483 : return *this;
484 : }
485 :
486 : static inline XMMReg4Int Zero()
487 : {
488 : XMMReg4Int reg;
489 : reg.xmm = _mm_setzero_si128();
490 : return reg;
491 : }
492 :
493 : static inline XMMReg4Int Set1(int i)
494 : {
495 : XMMReg4Int reg;
496 : reg.xmm = _mm_set1_epi32(i);
497 : return reg;
498 : }
499 :
500 289104 : static inline XMMReg4Int LoadAllVal(const int *ptr)
501 : {
502 289104 : return Load4Val(ptr);
503 : }
504 :
505 289104 : static inline XMMReg4Int Load4Val(const int *ptr)
506 : {
507 289104 : XMMReg4Int reg;
508 289104 : reg.nsLoad4Val(ptr);
509 289104 : return reg;
510 : }
511 :
512 289104 : inline void nsLoad4Val(const int *ptr)
513 : {
514 289104 : xmm = _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr));
515 289104 : }
516 :
517 : static inline XMMReg4Int Equals(const XMMReg4Int &expr1,
518 : const XMMReg4Int &expr2)
519 : {
520 : XMMReg4Int reg;
521 : reg.xmm = _mm_cmpeq_epi32(expr1.xmm, expr2.xmm);
522 : return reg;
523 : }
524 :
525 : static inline XMMReg4Int Ternary(const XMMReg4Int &cond,
526 : const XMMReg4Int &true_expr,
527 : const XMMReg4Int &false_expr)
528 : {
529 : XMMReg4Int reg;
530 : reg.xmm = _mm_or_si128(_mm_and_si128(cond.xmm, true_expr.xmm),
531 : _mm_andnot_si128(cond.xmm, false_expr.xmm));
532 : return reg;
533 : }
534 :
535 180690 : inline XMMReg4Int &operator+=(const XMMReg4Int &other)
536 : {
537 180690 : xmm = _mm_add_epi32(xmm, other.xmm);
538 180690 : return *this;
539 : }
540 :
541 36138 : inline XMMReg4Int &operator-=(const XMMReg4Int &other)
542 : {
543 36138 : xmm = _mm_sub_epi32(xmm, other.xmm);
544 36138 : return *this;
545 : }
546 :
547 : inline XMMReg4Int operator+(const XMMReg4Int &other) const
548 : {
549 : XMMReg4Int ret;
550 : ret.xmm = _mm_add_epi32(xmm, other.xmm);
551 : return ret;
552 : }
553 :
554 144552 : inline XMMReg4Int operator-(const XMMReg4Int &other) const
555 : {
556 144552 : XMMReg4Int ret;
557 144552 : ret.xmm = _mm_sub_epi32(xmm, other.xmm);
558 144552 : return ret;
559 : }
560 :
561 : XMMReg4Double cast_to_double() const;
562 :
563 611084 : XMMReg4Float cast_to_float() const
564 : {
565 611084 : XMMReg4Float ret;
566 611084 : ret.xmm = _mm_cvtepi32_ps(xmm);
567 611084 : return ret;
568 : }
569 : };
570 :
571 2694040 : inline XMMReg4Int XMMReg4Float::truncate_to_int() const
572 : {
573 2694040 : XMMReg4Int ret;
574 2694040 : ret.xmm = _mm_cvttps_epi32(xmm);
575 2694040 : return ret;
576 : }
577 :
578 : class XMMReg8Byte
579 : {
580 : public:
581 : __m128i xmm;
582 :
583 5779430 : XMMReg8Byte()
584 : #if !defined(_MSC_VER)
585 5779430 : : xmm(_mm_undefined_si128())
586 : #endif
587 : {
588 5779430 : }
589 :
590 : XMMReg8Byte(const XMMReg8Byte &other) : xmm(other.xmm)
591 : {
592 : }
593 :
594 269404 : static inline XMMReg8Byte Zero()
595 : {
596 269404 : XMMReg8Byte reg;
597 269404 : reg.xmm = _mm_setzero_si128();
598 269404 : return reg;
599 : }
600 :
601 269404 : static inline XMMReg8Byte Set1(char i)
602 : {
603 269404 : XMMReg8Byte reg;
604 269404 : reg.xmm = _mm_set1_epi8(i);
605 269404 : return reg;
606 : }
607 :
608 1347020 : static inline XMMReg8Byte Equals(const XMMReg8Byte &expr1,
609 : const XMMReg8Byte &expr2)
610 : {
611 1347020 : XMMReg8Byte reg;
612 1347020 : reg.xmm = _mm_cmpeq_epi8(expr1.xmm, expr2.xmm);
613 1347020 : return reg;
614 : }
615 :
616 489656 : static inline XMMReg8Byte Or(const XMMReg8Byte &expr1,
617 : const XMMReg8Byte &expr2)
618 : {
619 489656 : XMMReg8Byte reg;
620 489656 : reg.xmm = _mm_or_si128(expr1.xmm, expr2.xmm);
621 489656 : return reg;
622 : }
623 :
624 1248720 : static inline XMMReg8Byte Ternary(const XMMReg8Byte &cond,
625 : const XMMReg8Byte &true_expr,
626 : const XMMReg8Byte &false_expr)
627 : {
628 1248720 : XMMReg8Byte reg;
629 2497430 : reg.xmm = _mm_or_si128(_mm_and_si128(cond.xmm, true_expr.xmm),
630 1248720 : _mm_andnot_si128(cond.xmm, false_expr.xmm));
631 1248720 : return reg;
632 : }
633 :
634 538808 : inline XMMReg8Byte operator+(const XMMReg8Byte &other) const
635 : {
636 538808 : XMMReg8Byte ret;
637 538808 : ret.xmm = _mm_add_epi8(xmm, other.xmm);
638 538808 : return ret;
639 : }
640 :
641 269404 : inline XMMReg8Byte operator-(const XMMReg8Byte &other) const
642 : {
643 269404 : XMMReg8Byte ret;
644 269404 : ret.xmm = _mm_sub_epi8(xmm, other.xmm);
645 269404 : return ret;
646 : }
647 :
648 1347020 : static inline XMMReg8Byte Pack(const XMMReg4Int &r0, const XMMReg4Int &r1)
649 : {
650 1347020 : XMMReg8Byte reg;
651 1347020 : reg.xmm = _mm_packs_epi32(r0.xmm, r1.xmm);
652 1347020 : reg.xmm = _mm_packus_epi16(reg.xmm, reg.xmm);
653 1347020 : return reg;
654 : }
655 :
656 371338 : inline void Store8Val(unsigned char *ptr) const
657 : {
658 371338 : GDALCopyXMMToInt64(xmm, reinterpret_cast<GInt64 *>(ptr));
659 371338 : }
660 : };
661 :
662 : class XMMReg2Double
663 : {
664 : public:
665 : __m128d xmm;
666 :
667 4716800084 : XMMReg2Double()
668 : #if !defined(_MSC_VER)
669 4716800084 : : xmm(_mm_undefined_pd())
670 : #endif
671 : {
672 4716800084 : }
673 :
674 : XMMReg2Double(double val) : xmm(_mm_load_sd(&val))
675 : {
676 : }
677 :
678 2677980 : XMMReg2Double(const XMMReg2Double &other) : xmm(other.xmm)
679 : {
680 2677980 : }
681 :
682 222 : static inline XMMReg2Double Set1(double d)
683 : {
684 222 : XMMReg2Double reg;
685 222 : reg.xmm = _mm_set1_pd(d);
686 222 : return reg;
687 : }
688 :
689 1845030 : static inline XMMReg2Double Zero()
690 : {
691 1845030 : XMMReg2Double reg;
692 1845030 : reg.Zeroize();
693 1845030 : return reg;
694 : }
695 :
696 : static inline XMMReg2Double Load1ValHighAndLow(const double *ptr)
697 : {
698 : XMMReg2Double reg;
699 : reg.nsLoad1ValHighAndLow(ptr);
700 : return reg;
701 : }
702 :
703 145743 : static inline XMMReg2Double Load2Val(const double *ptr)
704 : {
705 145743 : XMMReg2Double reg;
706 145743 : reg.nsLoad2Val(ptr);
707 145743 : return reg;
708 : }
709 :
710 768 : static inline XMMReg2Double Load2Val(const float *ptr)
711 : {
712 768 : XMMReg2Double reg;
713 768 : reg.nsLoad2Val(ptr);
714 768 : return reg;
715 : }
716 :
717 23062900 : static inline XMMReg2Double Load2ValAligned(const double *ptr)
718 : {
719 23062900 : XMMReg2Double reg;
720 23062900 : reg.nsLoad2ValAligned(ptr);
721 23062900 : return reg;
722 : }
723 :
724 534521 : static inline XMMReg2Double Load2Val(const unsigned char *ptr)
725 : {
726 534521 : XMMReg2Double reg;
727 534519 : reg.nsLoad2Val(ptr);
728 534519 : return reg;
729 : }
730 :
731 18492 : static inline XMMReg2Double Load2Val(const short *ptr)
732 : {
733 18492 : XMMReg2Double reg;
734 18492 : reg.nsLoad2Val(ptr);
735 18492 : return reg;
736 : }
737 :
738 29184 : static inline XMMReg2Double Load2Val(const unsigned short *ptr)
739 : {
740 29184 : XMMReg2Double reg;
741 29184 : reg.nsLoad2Val(ptr);
742 29184 : return reg;
743 : }
744 :
745 : static inline XMMReg2Double Load2Val(const int *ptr)
746 : {
747 : XMMReg2Double reg;
748 : reg.nsLoad2Val(ptr);
749 : return reg;
750 : }
751 :
752 2 : static inline XMMReg2Double Equals(const XMMReg2Double &expr1,
753 : const XMMReg2Double &expr2)
754 : {
755 2 : XMMReg2Double reg;
756 2 : reg.xmm = _mm_cmpeq_pd(expr1.xmm, expr2.xmm);
757 2 : return reg;
758 : }
759 :
760 2653732 : static inline XMMReg2Double NotEquals(const XMMReg2Double &expr1,
761 : const XMMReg2Double &expr2)
762 : {
763 2653732 : XMMReg2Double reg;
764 2665002 : reg.xmm = _mm_cmpneq_pd(expr1.xmm, expr2.xmm);
765 2668712 : return reg;
766 : }
767 :
768 6 : static inline XMMReg2Double Greater(const XMMReg2Double &expr1,
769 : const XMMReg2Double &expr2)
770 : {
771 6 : XMMReg2Double reg;
772 6 : reg.xmm = _mm_cmpgt_pd(expr1.xmm, expr2.xmm);
773 6 : return reg;
774 : }
775 :
776 2663510 : static inline XMMReg2Double And(const XMMReg2Double &expr1,
777 : const XMMReg2Double &expr2)
778 : {
779 2663510 : XMMReg2Double reg;
780 2667100 : reg.xmm = _mm_and_pd(expr1.xmm, expr2.xmm);
781 2666600 : return reg;
782 : }
783 :
784 2 : static inline XMMReg2Double Ternary(const XMMReg2Double &cond,
785 : const XMMReg2Double &true_expr,
786 : const XMMReg2Double &false_expr)
787 : {
788 2 : XMMReg2Double reg;
789 4 : reg.xmm = _mm_or_pd(_mm_and_pd(cond.xmm, true_expr.xmm),
790 2 : _mm_andnot_pd(cond.xmm, false_expr.xmm));
791 2 : return reg;
792 : }
793 :
794 8270212 : static inline XMMReg2Double Min(const XMMReg2Double &expr1,
795 : const XMMReg2Double &expr2)
796 : {
797 8270212 : XMMReg2Double reg;
798 8288652 : reg.xmm = _mm_min_pd(expr1.xmm, expr2.xmm);
799 8363502 : return reg;
800 : }
801 :
802 121168001 : inline void nsLoad1ValHighAndLow(const double *ptr)
803 : {
804 121168001 : xmm = _mm_load1_pd(ptr);
805 121168001 : }
806 :
807 535477017 : inline void nsLoad2Val(const double *ptr)
808 : {
809 535477017 : xmm = _mm_loadu_pd(ptr);
810 535477017 : }
811 :
812 235933000 : inline void nsLoad2ValAligned(const double *ptr)
813 : {
814 235933000 : xmm = _mm_load_pd(ptr);
815 235933000 : }
816 :
817 768 : inline void nsLoad2Val(const float *ptr)
818 : {
819 1536 : xmm = _mm_cvtps_pd(_mm_castsi128_ps(GDALCopyInt64ToXMM(ptr)));
820 768 : }
821 :
822 2048 : inline void nsLoad2Val(const int *ptr)
823 : {
824 2048 : xmm = _mm_cvtepi32_pd(GDALCopyInt64ToXMM(ptr));
825 2048 : }
826 :
827 534518 : inline void nsLoad2Val(const unsigned char *ptr)
828 : {
829 534518 : __m128i xmm_i = GDALCopyInt16ToXMM(ptr);
830 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
831 : xmm_i = _mm_cvtepu8_epi32(xmm_i);
832 : #else
833 1069040 : xmm_i = _mm_unpacklo_epi8(xmm_i, _mm_setzero_si128());
834 1069040 : xmm_i = _mm_unpacklo_epi16(xmm_i, _mm_setzero_si128());
835 : #endif
836 534522 : xmm = _mm_cvtepi32_pd(xmm_i);
837 534522 : }
838 :
839 2048212 : inline void nsLoad2Val(const short *ptr)
840 : {
841 2048212 : __m128i xmm_i = GDALCopyInt32ToXMM(ptr);
842 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
843 : xmm_i = _mm_cvtepi16_epi32(xmm_i);
844 : #else
845 2048212 : xmm_i = _mm_unpacklo_epi16(
846 : xmm_i, xmm_i); /* 0|0|0|0|0|0|b|a --> 0|0|0|0|b|b|a|a */
847 2048212 : xmm_i = _mm_srai_epi32(
848 : xmm_i, 16); /* 0|0|0|0|b|b|a|a --> 0|0|0|0|sign(b)|b|sign(a)|a */
849 : #endif
850 2048212 : xmm = _mm_cvtepi32_pd(xmm_i);
851 2048212 : }
852 :
853 44347302 : inline void nsLoad2Val(const unsigned short *ptr)
854 : {
855 44347302 : __m128i xmm_i = GDALCopyInt32ToXMM(ptr);
856 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
857 : xmm_i = _mm_cvtepu16_epi32(xmm_i);
858 : #else
859 89418504 : xmm_i = _mm_unpacklo_epi16(
860 : xmm_i,
861 : _mm_setzero_si128()); /* 0|0|0|0|0|0|b|a --> 0|0|0|0|0|b|0|a */
862 : #endif
863 44837402 : xmm = _mm_cvtepi32_pd(xmm_i);
864 44837402 : }
865 :
866 370096001 : static inline void Load4Val(const unsigned char *ptr, XMMReg2Double &low,
867 : XMMReg2Double &high)
868 : {
869 370096001 : __m128i xmm_i = GDALCopyInt32ToXMM(ptr);
870 : #if defined(__SSE4_1__) || defined(__AVX__) || defined(USE_NEON_OPTIMIZATIONS)
871 : xmm_i = _mm_cvtepu8_epi32(xmm_i);
872 : #else
873 743631002 : xmm_i = _mm_unpacklo_epi8(xmm_i, _mm_setzero_si128());
874 744907002 : xmm_i = _mm_unpacklo_epi16(xmm_i, _mm_setzero_si128());
875 : #endif
876 373159001 : low.xmm = _mm_cvtepi32_pd(xmm_i);
877 373397001 : high.xmm =
878 373159001 : _mm_cvtepi32_pd(_mm_shuffle_epi32(xmm_i, _MM_SHUFFLE(3, 2, 3, 2)));
879 373397001 : }
880 :
881 : static inline void Load4Val(const short *ptr, XMMReg2Double &low,
882 : XMMReg2Double &high)
883 : {
884 : low.nsLoad2Val(ptr);
885 : high.nsLoad2Val(ptr + 2);
886 : }
887 :
888 : static inline void Load4Val(const unsigned short *ptr, XMMReg2Double &low,
889 : XMMReg2Double &high)
890 : {
891 : low.nsLoad2Val(ptr);
892 : high.nsLoad2Val(ptr + 2);
893 : }
894 :
895 : static inline void Load4Val(const double *ptr, XMMReg2Double &low,
896 : XMMReg2Double &high)
897 : {
898 : low.nsLoad2Val(ptr);
899 : high.nsLoad2Val(ptr + 2);
900 : }
901 :
902 38657 : static inline void Load4Val(const float *ptr, XMMReg2Double &low,
903 : XMMReg2Double &high)
904 : {
905 38657 : __m128 temp1 = _mm_loadu_ps(ptr);
906 38657 : __m128 temp2 = _mm_shuffle_ps(temp1, temp1, _MM_SHUFFLE(3, 2, 3, 2));
907 38657 : low.xmm = _mm_cvtps_pd(temp1);
908 38657 : high.xmm = _mm_cvtps_pd(temp2);
909 38657 : }
910 :
911 505071000 : inline void Zeroize()
912 : {
913 505071000 : xmm = _mm_setzero_pd();
914 505071000 : }
915 :
916 1609220037 : inline XMMReg2Double &operator=(const XMMReg2Double &other)
917 : {
918 1609220037 : xmm = other.xmm;
919 1609220037 : return *this;
920 : }
921 :
922 1176810003 : inline XMMReg2Double &operator+=(const XMMReg2Double &other)
923 : {
924 1176810003 : xmm = _mm_add_pd(xmm, other.xmm);
925 1176810003 : return *this;
926 : }
927 :
928 22726802 : inline XMMReg2Double &operator*=(const XMMReg2Double &other)
929 : {
930 22726802 : xmm = _mm_mul_pd(xmm, other.xmm);
931 22726802 : return *this;
932 : }
933 :
934 242721009 : inline XMMReg2Double operator+(const XMMReg2Double &other) const
935 : {
936 242721009 : XMMReg2Double ret;
937 242716009 : ret.xmm = _mm_add_pd(xmm, other.xmm);
938 242716009 : return ret;
939 : }
940 :
941 2 : inline XMMReg2Double operator-(const XMMReg2Double &other) const
942 : {
943 2 : XMMReg2Double ret;
944 2 : ret.xmm = _mm_sub_pd(xmm, other.xmm);
945 2 : return ret;
946 : }
947 :
948 1228170002 : inline XMMReg2Double operator*(const XMMReg2Double &other) const
949 : {
950 1228170002 : XMMReg2Double ret;
951 1232550002 : ret.xmm = _mm_mul_pd(xmm, other.xmm);
952 1232550002 : return ret;
953 : }
954 :
955 2653452 : inline XMMReg2Double operator/(const XMMReg2Double &other) const
956 : {
957 2653452 : XMMReg2Double ret;
958 2656472 : ret.xmm = _mm_div_pd(xmm, other.xmm);
959 2656472 : return ret;
960 : }
961 :
962 244057001 : inline double GetHorizSum() const
963 : {
964 : __m128d xmm2;
965 244057001 : xmm2 = _mm_shuffle_pd(
966 : xmm, xmm,
967 : _MM_SHUFFLE2(0, 1)); /* transfer high word into low word of xmm2 */
968 734181003 : return _mm_cvtsd_f64(_mm_add_sd(xmm, xmm2));
969 : }
970 :
971 7200 : inline void Store2Val(double *ptr) const
972 : {
973 7200 : _mm_storeu_pd(ptr, xmm);
974 7200 : }
975 :
976 : inline void Store2ValAligned(double *ptr) const
977 : {
978 : _mm_store_pd(ptr, xmm);
979 : }
980 :
981 101510002 : inline void Store2Val(float *ptr) const
982 : {
983 203190004 : __m128i xmm_i = _mm_castps_si128(_mm_cvtpd_ps(xmm));
984 101680002 : GDALCopyXMMToInt64(xmm_i, reinterpret_cast<GInt64 *>(ptr));
985 101540002 : }
986 :
987 : inline void Store2Val(unsigned char *ptr) const
988 : {
989 : __m128i tmp = _mm_cvttpd_epi32(_mm_add_pd(
990 : xmm,
991 : _mm_set1_pd(0.5))); /* Convert the 2 double values to 2 integers */
992 : tmp = _mm_packs_epi32(tmp, tmp);
993 : tmp = _mm_packus_epi16(tmp, tmp);
994 : GDALCopyXMMToInt16(tmp, reinterpret_cast<GInt16 *>(ptr));
995 : }
996 :
997 5089662 : inline void Store2Val(unsigned short *ptr) const
998 : {
999 5089662 : __m128i tmp = _mm_cvttpd_epi32(_mm_add_pd(
1000 5089662 : xmm,
1001 : _mm_set1_pd(0.5))); /* Convert the 2 double values to 2 integers */
1002 : // X X X X 0 B 0 A --> X X X X A A B A
1003 5075702 : tmp = _mm_shufflelo_epi16(tmp, 0 | (2 << 2));
1004 5144302 : GDALCopyXMMToInt32(tmp, reinterpret_cast<GInt32 *>(ptr));
1005 5074742 : }
1006 :
1007 8 : inline void StoreMask(unsigned char *ptr) const
1008 : {
1009 8 : _mm_storeu_si128(reinterpret_cast<__m128i *>(ptr),
1010 8 : _mm_castpd_si128(xmm));
1011 8 : }
1012 :
1013 : inline operator double() const
1014 : {
1015 : return _mm_cvtsd_f64(xmm);
1016 : }
1017 : };
1018 :
1019 : #else
1020 :
1021 : #ifndef NO_WARN_USE_SSE2_EMULATION
1022 : #warning "Software emulation of SSE2 !"
1023 : #endif
1024 :
1025 : class XMMReg2Double
1026 : {
1027 : public:
1028 : double low;
1029 : double high;
1030 :
1031 : // cppcheck-suppress uninitMemberVar
1032 : XMMReg2Double() = default;
1033 :
1034 : explicit XMMReg2Double(double val)
1035 : {
1036 : low = val;
1037 : high = 0.0;
1038 : }
1039 :
1040 : XMMReg2Double(const XMMReg2Double &other) : low(other.low), high(other.high)
1041 : {
1042 : }
1043 :
1044 : static inline XMMReg2Double Zero()
1045 : {
1046 : XMMReg2Double reg;
1047 : reg.Zeroize();
1048 : return reg;
1049 : }
1050 :
1051 : static inline XMMReg2Double Set1(double d)
1052 : {
1053 : XMMReg2Double reg;
1054 : reg.low = d;
1055 : reg.high = d;
1056 : return reg;
1057 : }
1058 :
1059 : static inline XMMReg2Double Load1ValHighAndLow(const double *ptr)
1060 : {
1061 : XMMReg2Double reg;
1062 : reg.nsLoad1ValHighAndLow(ptr);
1063 : return reg;
1064 : }
1065 :
1066 2 : static inline XMMReg2Double Equals(const XMMReg2Double &expr1,
1067 : const XMMReg2Double &expr2)
1068 : {
1069 : XMMReg2Double reg;
1070 :
1071 2 : if (expr1.low == expr2.low)
1072 2 : memset(&(reg.low), 0xFF, sizeof(double));
1073 : else
1074 0 : reg.low = 0;
1075 :
1076 2 : if (expr1.high == expr2.high)
1077 2 : memset(&(reg.high), 0xFF, sizeof(double));
1078 : else
1079 0 : reg.high = 0;
1080 :
1081 2 : return reg;
1082 : }
1083 :
1084 2 : static inline XMMReg2Double NotEquals(const XMMReg2Double &expr1,
1085 : const XMMReg2Double &expr2)
1086 : {
1087 : XMMReg2Double reg;
1088 :
1089 2 : if (expr1.low != expr2.low)
1090 0 : memset(&(reg.low), 0xFF, sizeof(double));
1091 : else
1092 2 : reg.low = 0;
1093 :
1094 2 : if (expr1.high != expr2.high)
1095 0 : memset(&(reg.high), 0xFF, sizeof(double));
1096 : else
1097 2 : reg.high = 0;
1098 :
1099 2 : return reg;
1100 : }
1101 :
1102 6 : static inline XMMReg2Double Greater(const XMMReg2Double &expr1,
1103 : const XMMReg2Double &expr2)
1104 : {
1105 : XMMReg2Double reg;
1106 :
1107 6 : if (expr1.low > expr2.low)
1108 2 : memset(&(reg.low), 0xFF, sizeof(double));
1109 : else
1110 4 : reg.low = 0;
1111 :
1112 6 : if (expr1.high > expr2.high)
1113 2 : memset(&(reg.high), 0xFF, sizeof(double));
1114 : else
1115 4 : reg.high = 0;
1116 :
1117 6 : return reg;
1118 : }
1119 :
1120 : static inline XMMReg2Double And(const XMMReg2Double &expr1,
1121 : const XMMReg2Double &expr2)
1122 : {
1123 : XMMReg2Double reg;
1124 : int low1[2], high1[2];
1125 : int low2[2], high2[2];
1126 : memcpy(low1, &expr1.low, sizeof(double));
1127 : memcpy(high1, &expr1.high, sizeof(double));
1128 : memcpy(low2, &expr2.low, sizeof(double));
1129 : memcpy(high2, &expr2.high, sizeof(double));
1130 : low1[0] &= low2[0];
1131 : low1[1] &= low2[1];
1132 : high1[0] &= high2[0];
1133 : high1[1] &= high2[1];
1134 : memcpy(®.low, low1, sizeof(double));
1135 : memcpy(®.high, high1, sizeof(double));
1136 : return reg;
1137 : }
1138 :
1139 2 : static inline XMMReg2Double Ternary(const XMMReg2Double &cond,
1140 : const XMMReg2Double &true_expr,
1141 : const XMMReg2Double &false_expr)
1142 : {
1143 : XMMReg2Double reg;
1144 2 : if (cond.low != 0)
1145 1 : reg.low = true_expr.low;
1146 : else
1147 1 : reg.low = false_expr.low;
1148 2 : if (cond.high != 0)
1149 1 : reg.high = true_expr.high;
1150 : else
1151 1 : reg.high = false_expr.high;
1152 2 : return reg;
1153 : }
1154 :
1155 2 : static inline XMMReg2Double Min(const XMMReg2Double &expr1,
1156 : const XMMReg2Double &expr2)
1157 : {
1158 : XMMReg2Double reg;
1159 2 : reg.low = (expr1.low < expr2.low) ? expr1.low : expr2.low;
1160 2 : reg.high = (expr1.high < expr2.high) ? expr1.high : expr2.high;
1161 2 : return reg;
1162 : }
1163 :
1164 1 : static inline XMMReg2Double Load2Val(const double *ptr)
1165 : {
1166 : XMMReg2Double reg;
1167 1 : reg.nsLoad2Val(ptr);
1168 1 : return reg;
1169 : }
1170 :
1171 : static inline XMMReg2Double Load2ValAligned(const double *ptr)
1172 : {
1173 : XMMReg2Double reg;
1174 : reg.nsLoad2ValAligned(ptr);
1175 : return reg;
1176 : }
1177 :
1178 : static inline XMMReg2Double Load2Val(const float *ptr)
1179 : {
1180 : XMMReg2Double reg;
1181 : reg.nsLoad2Val(ptr);
1182 : return reg;
1183 : }
1184 :
1185 : static inline XMMReg2Double Load2Val(const unsigned char *ptr)
1186 : {
1187 : XMMReg2Double reg;
1188 : reg.nsLoad2Val(ptr);
1189 : return reg;
1190 : }
1191 :
1192 : static inline XMMReg2Double Load2Val(const short *ptr)
1193 : {
1194 : XMMReg2Double reg;
1195 : reg.nsLoad2Val(ptr);
1196 : return reg;
1197 : }
1198 :
1199 : static inline XMMReg2Double Load2Val(const unsigned short *ptr)
1200 : {
1201 : XMMReg2Double reg;
1202 : reg.nsLoad2Val(ptr);
1203 : return reg;
1204 : }
1205 :
1206 : static inline XMMReg2Double Load2Val(const int *ptr)
1207 : {
1208 : XMMReg2Double reg;
1209 : reg.nsLoad2Val(ptr);
1210 : return reg;
1211 : }
1212 :
1213 1 : inline void nsLoad1ValHighAndLow(const double *ptr)
1214 : {
1215 1 : low = ptr[0];
1216 1 : high = ptr[0];
1217 1 : }
1218 :
1219 17 : inline void nsLoad2Val(const double *ptr)
1220 : {
1221 17 : low = ptr[0];
1222 17 : high = ptr[1];
1223 17 : }
1224 :
1225 : inline void nsLoad2ValAligned(const double *ptr)
1226 : {
1227 : low = ptr[0];
1228 : high = ptr[1];
1229 : }
1230 :
1231 2 : inline void nsLoad2Val(const float *ptr)
1232 : {
1233 2 : low = ptr[0];
1234 2 : high = ptr[1];
1235 2 : }
1236 :
1237 : inline void nsLoad2Val(const unsigned char *ptr)
1238 : {
1239 : low = ptr[0];
1240 : high = ptr[1];
1241 : }
1242 :
1243 2 : inline void nsLoad2Val(const short *ptr)
1244 : {
1245 2 : low = ptr[0];
1246 2 : high = ptr[1];
1247 2 : }
1248 :
1249 2 : inline void nsLoad2Val(const unsigned short *ptr)
1250 : {
1251 2 : low = ptr[0];
1252 2 : high = ptr[1];
1253 2 : }
1254 :
1255 : inline void nsLoad2Val(const int *ptr)
1256 : {
1257 : low = ptr[0];
1258 : high = ptr[1];
1259 : }
1260 :
1261 1 : static inline void Load4Val(const unsigned char *ptr, XMMReg2Double &low,
1262 : XMMReg2Double &high)
1263 : {
1264 1 : low.low = ptr[0];
1265 1 : low.high = ptr[1];
1266 1 : high.low = ptr[2];
1267 1 : high.high = ptr[3];
1268 1 : }
1269 :
1270 : static inline void Load4Val(const short *ptr, XMMReg2Double &low,
1271 : XMMReg2Double &high)
1272 : {
1273 : low.nsLoad2Val(ptr);
1274 : high.nsLoad2Val(ptr + 2);
1275 : }
1276 :
1277 : static inline void Load4Val(const unsigned short *ptr, XMMReg2Double &low,
1278 : XMMReg2Double &high)
1279 : {
1280 : low.nsLoad2Val(ptr);
1281 : high.nsLoad2Val(ptr + 2);
1282 : }
1283 :
1284 : static inline void Load4Val(const double *ptr, XMMReg2Double &low,
1285 : XMMReg2Double &high)
1286 : {
1287 : low.nsLoad2Val(ptr);
1288 : high.nsLoad2Val(ptr + 2);
1289 : }
1290 :
1291 1 : static inline void Load4Val(const float *ptr, XMMReg2Double &low,
1292 : XMMReg2Double &high)
1293 : {
1294 1 : low.nsLoad2Val(ptr);
1295 1 : high.nsLoad2Val(ptr + 2);
1296 1 : }
1297 :
1298 : inline void Zeroize()
1299 : {
1300 : low = 0.0;
1301 : high = 0.0;
1302 : }
1303 :
1304 37 : inline XMMReg2Double &operator=(const XMMReg2Double &other)
1305 : {
1306 37 : low = other.low;
1307 37 : high = other.high;
1308 37 : return *this;
1309 : }
1310 :
1311 3 : inline XMMReg2Double &operator+=(const XMMReg2Double &other)
1312 : {
1313 3 : low += other.low;
1314 3 : high += other.high;
1315 3 : return *this;
1316 : }
1317 :
1318 2 : inline XMMReg2Double &operator*=(const XMMReg2Double &other)
1319 : {
1320 2 : low *= other.low;
1321 2 : high *= other.high;
1322 2 : return *this;
1323 : }
1324 :
1325 9 : inline XMMReg2Double operator+(const XMMReg2Double &other) const
1326 : {
1327 : XMMReg2Double ret;
1328 9 : ret.low = low + other.low;
1329 9 : ret.high = high + other.high;
1330 9 : return ret;
1331 : }
1332 :
1333 2 : inline XMMReg2Double operator-(const XMMReg2Double &other) const
1334 : {
1335 : XMMReg2Double ret;
1336 2 : ret.low = low - other.low;
1337 2 : ret.high = high - other.high;
1338 2 : return ret;
1339 : }
1340 :
1341 2 : inline XMMReg2Double operator*(const XMMReg2Double &other) const
1342 : {
1343 : XMMReg2Double ret;
1344 2 : ret.low = low * other.low;
1345 2 : ret.high = high * other.high;
1346 2 : return ret;
1347 : }
1348 :
1349 2 : inline XMMReg2Double operator/(const XMMReg2Double &other) const
1350 : {
1351 : XMMReg2Double ret;
1352 2 : ret.low = low / other.low;
1353 2 : ret.high = high / other.high;
1354 2 : return ret;
1355 : }
1356 :
1357 1 : inline double GetHorizSum() const
1358 : {
1359 1 : return low + high;
1360 : }
1361 :
1362 32 : inline void Store2Val(double *ptr) const
1363 : {
1364 32 : ptr[0] = low;
1365 32 : ptr[1] = high;
1366 32 : }
1367 :
1368 : inline void Store2ValAligned(double *ptr) const
1369 : {
1370 : ptr[0] = low;
1371 : ptr[1] = high;
1372 : }
1373 :
1374 2 : inline void Store2Val(float *ptr) const
1375 : {
1376 2 : ptr[0] = static_cast<float>(low);
1377 2 : ptr[1] = static_cast<float>(high);
1378 2 : }
1379 :
1380 2 : void Store2Val(unsigned char *ptr) const
1381 : {
1382 2 : ptr[0] = static_cast<unsigned char>(low + 0.5);
1383 2 : ptr[1] = static_cast<unsigned char>(high + 0.5);
1384 2 : }
1385 :
1386 2 : void Store2Val(unsigned short *ptr) const
1387 : {
1388 2 : ptr[0] = static_cast<GUInt16>(low + 0.5);
1389 2 : ptr[1] = static_cast<GUInt16>(high + 0.5);
1390 2 : }
1391 :
1392 8 : inline void StoreMask(unsigned char *ptr) const
1393 : {
1394 8 : memcpy(ptr, &low, 8);
1395 8 : memcpy(ptr + 8, &high, 8);
1396 8 : }
1397 :
1398 : inline operator double() const
1399 : {
1400 : return low;
1401 : }
1402 : };
1403 :
1404 : #endif /* defined(__x86_64) || defined(_M_X64) */
1405 :
1406 : #if defined(__AVX__) && !defined(USE_SSE2_EMULATION)
1407 :
1408 : #include <immintrin.h>
1409 :
1410 : class XMMReg4Double
1411 : {
1412 : public:
1413 : __m256d ymm;
1414 :
1415 : XMMReg4Double() : ymm(_mm256_setzero_pd())
1416 : {
1417 : }
1418 :
1419 : XMMReg4Double(const XMMReg4Double &other) : ymm(other.ymm)
1420 : {
1421 : }
1422 :
1423 : static inline XMMReg4Double Zero()
1424 : {
1425 : XMMReg4Double reg;
1426 : reg.Zeroize();
1427 : return reg;
1428 : }
1429 :
1430 : static inline XMMReg4Double Set1(double d)
1431 : {
1432 : XMMReg4Double reg;
1433 : reg.ymm = _mm256_set1_pd(d);
1434 : return reg;
1435 : }
1436 :
1437 : inline void Zeroize()
1438 : {
1439 : ymm = _mm256_setzero_pd();
1440 : }
1441 :
1442 : static inline XMMReg4Double Load1ValHighAndLow(const double *ptr)
1443 : {
1444 : XMMReg4Double reg;
1445 : reg.nsLoad1ValHighAndLow(ptr);
1446 : return reg;
1447 : }
1448 :
1449 : inline void nsLoad1ValHighAndLow(const double *ptr)
1450 : {
1451 : ymm = _mm256_set1_pd(*ptr);
1452 : }
1453 :
1454 : static inline XMMReg4Double Load4Val(const unsigned char *ptr)
1455 : {
1456 : XMMReg4Double reg;
1457 : reg.nsLoad4Val(ptr);
1458 : return reg;
1459 : }
1460 :
1461 : inline void nsLoad4Val(const unsigned char *ptr)
1462 : {
1463 : __m128i xmm_i = GDALCopyInt32ToXMM(ptr);
1464 : xmm_i = _mm_cvtepu8_epi32(xmm_i);
1465 : ymm = _mm256_cvtepi32_pd(xmm_i);
1466 : }
1467 :
1468 : static inline void Load8Val(const unsigned char *ptr, XMMReg4Double &low,
1469 : XMMReg4Double &high)
1470 : {
1471 : const __m128i xmm_i = GDALCopyInt64ToXMM(ptr);
1472 : const __m128i xmm_i_low = _mm_cvtepu8_epi32(xmm_i);
1473 : low.ymm = _mm256_cvtepi32_pd(xmm_i_low);
1474 : const __m128i xmm_i_high = _mm_cvtepu8_epi32(_mm_srli_si128(xmm_i, 4));
1475 : high.ymm = _mm256_cvtepi32_pd(xmm_i_high);
1476 : }
1477 :
1478 : static inline XMMReg4Double Load4Val(const short *ptr)
1479 : {
1480 : XMMReg4Double reg;
1481 : reg.nsLoad4Val(ptr);
1482 : return reg;
1483 : }
1484 :
1485 : inline void nsLoad4Val(const short *ptr)
1486 : {
1487 : __m128i xmm_i = GDALCopyInt64ToXMM(ptr);
1488 : xmm_i = _mm_cvtepi16_epi32(xmm_i);
1489 : ymm = _mm256_cvtepi32_pd(xmm_i);
1490 : }
1491 :
1492 : static inline void Load8Val(const short *ptr, XMMReg4Double &low,
1493 : XMMReg4Double &high)
1494 : {
1495 : low.nsLoad4Val(ptr);
1496 : high.nsLoad4Val(ptr + 4);
1497 : }
1498 :
1499 : static inline XMMReg4Double Load4Val(const unsigned short *ptr)
1500 : {
1501 : XMMReg4Double reg;
1502 : reg.nsLoad4Val(ptr);
1503 : return reg;
1504 : }
1505 :
1506 : inline void nsLoad4Val(const unsigned short *ptr)
1507 : {
1508 : __m128i xmm_i = GDALCopyInt64ToXMM(ptr);
1509 : xmm_i = _mm_cvtepu16_epi32(xmm_i);
1510 : ymm = _mm256_cvtepi32_pd(
1511 : xmm_i); // ok to use signed conversion since we are in the ushort
1512 : // range, so cannot be interpreted as negative int32
1513 : }
1514 :
1515 : static inline void Load8Val(const unsigned short *ptr, XMMReg4Double &low,
1516 : XMMReg4Double &high)
1517 : {
1518 : low.nsLoad4Val(ptr);
1519 : high.nsLoad4Val(ptr + 4);
1520 : }
1521 :
1522 : static inline XMMReg4Double Load4Val(const double *ptr)
1523 : {
1524 : XMMReg4Double reg;
1525 : reg.nsLoad4Val(ptr);
1526 : return reg;
1527 : }
1528 :
1529 : inline void nsLoad4Val(const double *ptr)
1530 : {
1531 : ymm = _mm256_loadu_pd(ptr);
1532 : }
1533 :
1534 : static inline void Load8Val(const double *ptr, XMMReg4Double &low,
1535 : XMMReg4Double &high)
1536 : {
1537 : low.nsLoad4Val(ptr);
1538 : high.nsLoad4Val(ptr + 4);
1539 : }
1540 :
1541 : static inline XMMReg4Double Load4ValAligned(const double *ptr)
1542 : {
1543 : XMMReg4Double reg;
1544 : reg.nsLoad4ValAligned(ptr);
1545 : return reg;
1546 : }
1547 :
1548 : inline void nsLoad4ValAligned(const double *ptr)
1549 : {
1550 : ymm = _mm256_load_pd(ptr);
1551 : }
1552 :
1553 : static inline XMMReg4Double Load4Val(const float *ptr)
1554 : {
1555 : XMMReg4Double reg;
1556 : reg.nsLoad4Val(ptr);
1557 : return reg;
1558 : }
1559 :
1560 : inline void nsLoad4Val(const float *ptr)
1561 : {
1562 : ymm = _mm256_cvtps_pd(_mm_loadu_ps(ptr));
1563 : }
1564 :
1565 : static inline void Load8Val(const float *ptr, XMMReg4Double &low,
1566 : XMMReg4Double &high)
1567 : {
1568 : low.nsLoad4Val(ptr);
1569 : high.nsLoad4Val(ptr + 4);
1570 : }
1571 :
1572 : static inline XMMReg4Double Load4Val(const int *ptr)
1573 : {
1574 : XMMReg4Double reg;
1575 : reg.nsLoad4Val(ptr);
1576 : return reg;
1577 : }
1578 :
1579 : inline void nsLoad4Val(const int *ptr)
1580 : {
1581 : ymm = _mm256_cvtepi32_pd(
1582 : _mm_loadu_si128(reinterpret_cast<const __m128i *>(ptr)));
1583 : }
1584 :
1585 : static inline void Load8Val(const int *ptr, XMMReg4Double &low,
1586 : XMMReg4Double &high)
1587 : {
1588 : low.nsLoad4Val(ptr);
1589 : high.nsLoad4Val(ptr + 4);
1590 : }
1591 :
1592 : static inline XMMReg4Double Equals(const XMMReg4Double &expr1,
1593 : const XMMReg4Double &expr2)
1594 : {
1595 : XMMReg4Double reg;
1596 : reg.ymm = _mm256_cmp_pd(expr1.ymm, expr2.ymm, _CMP_EQ_OQ);
1597 : return reg;
1598 : }
1599 :
1600 : static inline XMMReg4Double NotEquals(const XMMReg4Double &expr1,
1601 : const XMMReg4Double &expr2)
1602 : {
1603 : XMMReg4Double reg;
1604 : reg.ymm = _mm256_cmp_pd(expr1.ymm, expr2.ymm, _CMP_NEQ_OQ);
1605 : return reg;
1606 : }
1607 :
1608 : static inline XMMReg4Double Greater(const XMMReg4Double &expr1,
1609 : const XMMReg4Double &expr2)
1610 : {
1611 : XMMReg4Double reg;
1612 : reg.ymm = _mm256_cmp_pd(expr1.ymm, expr2.ymm, _CMP_GT_OQ);
1613 : return reg;
1614 : }
1615 :
1616 : static inline XMMReg4Double And(const XMMReg4Double &expr1,
1617 : const XMMReg4Double &expr2)
1618 : {
1619 : XMMReg4Double reg;
1620 : reg.ymm = _mm256_and_pd(expr1.ymm, expr2.ymm);
1621 : return reg;
1622 : }
1623 :
1624 : static inline XMMReg4Double Ternary(const XMMReg4Double &cond,
1625 : const XMMReg4Double &true_expr,
1626 : const XMMReg4Double &false_expr)
1627 : {
1628 : XMMReg4Double reg;
1629 : reg.ymm = _mm256_or_pd(_mm256_and_pd(cond.ymm, true_expr.ymm),
1630 : _mm256_andnot_pd(cond.ymm, false_expr.ymm));
1631 : return reg;
1632 : }
1633 :
1634 : static inline XMMReg4Double Min(const XMMReg4Double &expr1,
1635 : const XMMReg4Double &expr2)
1636 : {
1637 : XMMReg4Double reg;
1638 : reg.ymm = _mm256_min_pd(expr1.ymm, expr2.ymm);
1639 : return reg;
1640 : }
1641 :
1642 : inline XMMReg4Double &operator=(const XMMReg4Double &other)
1643 : {
1644 : ymm = other.ymm;
1645 : return *this;
1646 : }
1647 :
1648 : inline XMMReg4Double &operator+=(const XMMReg4Double &other)
1649 : {
1650 : ymm = _mm256_add_pd(ymm, other.ymm);
1651 : return *this;
1652 : }
1653 :
1654 : inline XMMReg4Double &operator*=(const XMMReg4Double &other)
1655 : {
1656 : ymm = _mm256_mul_pd(ymm, other.ymm);
1657 : return *this;
1658 : }
1659 :
1660 : inline XMMReg4Double operator+(const XMMReg4Double &other) const
1661 : {
1662 : XMMReg4Double ret;
1663 : ret.ymm = _mm256_add_pd(ymm, other.ymm);
1664 : return ret;
1665 : }
1666 :
1667 : inline XMMReg4Double operator-(const XMMReg4Double &other) const
1668 : {
1669 : XMMReg4Double ret;
1670 : ret.ymm = _mm256_sub_pd(ymm, other.ymm);
1671 : return ret;
1672 : }
1673 :
1674 : inline XMMReg4Double operator*(const XMMReg4Double &other) const
1675 : {
1676 : XMMReg4Double ret;
1677 : ret.ymm = _mm256_mul_pd(ymm, other.ymm);
1678 : return ret;
1679 : }
1680 :
1681 : inline XMMReg4Double operator/(const XMMReg4Double &other) const
1682 : {
1683 : XMMReg4Double ret;
1684 : ret.ymm = _mm256_div_pd(ymm, other.ymm);
1685 : return ret;
1686 : }
1687 :
1688 : void AddToLow(const XMMReg2Double &other)
1689 : {
1690 : __m256d ymm2 = _mm256_setzero_pd();
1691 : ymm2 = _mm256_insertf128_pd(ymm2, other.xmm, 0);
1692 : ymm = _mm256_add_pd(ymm, ymm2);
1693 : }
1694 :
1695 : inline double GetHorizSum() const
1696 : {
1697 : __m256d ymm_tmp1, ymm_tmp2;
1698 : ymm_tmp2 = _mm256_hadd_pd(ymm, ymm);
1699 : ymm_tmp1 = _mm256_permute2f128_pd(ymm_tmp2, ymm_tmp2, 1);
1700 : ymm_tmp1 = _mm256_add_pd(ymm_tmp1, ymm_tmp2);
1701 : return _mm_cvtsd_f64(_mm256_castpd256_pd128(ymm_tmp1));
1702 : }
1703 :
1704 : inline XMMReg4Double approx_inv_sqrt(const XMMReg4Double &one,
1705 : const XMMReg4Double &half) const
1706 : {
1707 : __m256d reg = ymm;
1708 : __m256d reg_half = _mm256_mul_pd(reg, half.ymm);
1709 : // Compute rough approximation of 1 / sqrt(b) with _mm_rsqrt_ps
1710 : reg = _mm256_cvtps_pd(_mm_rsqrt_ps(_mm256_cvtpd_ps(reg)));
1711 : // And perform one step of Newton-Raphson approximation to improve it
1712 : // approx_inv_sqrt_x = approx_inv_sqrt_x*(1.5 -
1713 : // 0.5*x*approx_inv_sqrt_x*approx_inv_sqrt_x);
1714 : const __m256d one_and_a_half = _mm256_add_pd(one.ymm, half.ymm);
1715 : reg = _mm256_mul_pd(
1716 : reg,
1717 : _mm256_sub_pd(one_and_a_half,
1718 : _mm256_mul_pd(reg_half, _mm256_mul_pd(reg, reg))));
1719 : XMMReg4Double ret;
1720 : ret.ymm = reg;
1721 : return ret;
1722 : }
1723 :
1724 : inline XMMReg4Float cast_to_float() const
1725 : {
1726 : XMMReg4Float ret;
1727 : ret.xmm = _mm256_cvtpd_ps(ymm);
1728 : return ret;
1729 : }
1730 :
1731 : inline void Store4Val(unsigned char *ptr) const
1732 : {
1733 : __m128i xmm_i =
1734 : _mm256_cvttpd_epi32(_mm256_add_pd(ymm, _mm256_set1_pd(0.5)));
1735 : // xmm_i = _mm_packs_epi32(xmm_i, xmm_i); // Pack int32 to int16
1736 : // xmm_i = _mm_packus_epi16(xmm_i, xmm_i); // Pack int16 to uint8
1737 : xmm_i =
1738 : _mm_shuffle_epi8(xmm_i, _mm_cvtsi32_si128(0 | (4 << 8) | (8 << 16) |
1739 : (12 << 24))); // SSSE3
1740 : GDALCopyXMMToInt32(xmm_i, reinterpret_cast<GInt32 *>(ptr));
1741 : }
1742 :
1743 : inline void Store4Val(unsigned short *ptr) const
1744 : {
1745 : __m128i xmm_i =
1746 : _mm256_cvttpd_epi32(_mm256_add_pd(ymm, _mm256_set1_pd(0.5)));
1747 : xmm_i = _mm_packus_epi32(xmm_i, xmm_i); // Pack uint32 to uint16
1748 : GDALCopyXMMToInt64(xmm_i, reinterpret_cast<GInt64 *>(ptr));
1749 : }
1750 :
1751 : inline void Store4Val(float *ptr) const
1752 : {
1753 : _mm_storeu_ps(ptr, _mm256_cvtpd_ps(ymm));
1754 : }
1755 :
1756 : inline void Store4Val(double *ptr) const
1757 : {
1758 : _mm256_storeu_pd(ptr, ymm);
1759 : }
1760 :
1761 : inline void StoreMask(unsigned char *ptr) const
1762 : {
1763 : _mm256_storeu_si256(reinterpret_cast<__m256i *>(ptr),
1764 : _mm256_castpd_si256(ymm));
1765 : }
1766 : };
1767 :
1768 : inline XMMReg4Double XMMReg4Float::cast_to_double() const
1769 : {
1770 : XMMReg4Double ret;
1771 : ret.ymm = _mm256_cvtps_pd(xmm);
1772 : return ret;
1773 : }
1774 :
1775 : inline XMMReg4Double XMMReg4Int::cast_to_double() const
1776 : {
1777 : XMMReg4Double ret;
1778 : ret.ymm = _mm256_cvtepi32_pd(xmm);
1779 : return ret;
1780 : }
1781 :
1782 : class XMMReg8Float
1783 : {
1784 : public:
1785 : __m256 ymm;
1786 :
1787 : XMMReg8Float()
1788 : #if !defined(_MSC_VER)
1789 : : ymm(_mm256_undefined_ps())
1790 : #endif
1791 : {
1792 : }
1793 :
1794 : XMMReg8Float(const XMMReg8Float &other) : ymm(other.ymm)
1795 : {
1796 : }
1797 :
1798 : static inline XMMReg8Float Set1(float f)
1799 : {
1800 : XMMReg8Float reg;
1801 : reg.ymm = _mm256_set1_ps(f);
1802 : return reg;
1803 : }
1804 :
1805 : static inline XMMReg8Float LoadAllVal(const float *ptr)
1806 : {
1807 : return Load8Val(ptr);
1808 : }
1809 :
1810 : static inline XMMReg8Float Load8Val(const float *ptr)
1811 : {
1812 : XMMReg8Float reg;
1813 : reg.nsLoad8Val(ptr);
1814 : return reg;
1815 : }
1816 :
1817 : static inline XMMReg8Float Load8Val(const int *ptr)
1818 : {
1819 : XMMReg8Float reg;
1820 : reg.nsLoad8Val(ptr);
1821 : return reg;
1822 : }
1823 :
1824 : static inline XMMReg8Float Max(const XMMReg8Float &expr1,
1825 : const XMMReg8Float &expr2)
1826 : {
1827 : XMMReg8Float reg;
1828 : reg.ymm = _mm256_max_ps(expr1.ymm, expr2.ymm);
1829 : return reg;
1830 : }
1831 :
1832 : inline void nsLoad8Val(const float *ptr)
1833 : {
1834 : ymm = _mm256_loadu_ps(ptr);
1835 : }
1836 :
1837 : inline void nsLoad8Val(const int *ptr)
1838 : {
1839 : const __m256i ymm_i =
1840 : _mm256_loadu_si256(reinterpret_cast<const __m256i *>(ptr));
1841 : ymm = _mm256_cvtepi32_ps(ymm_i);
1842 : }
1843 :
1844 : inline XMMReg8Float &operator=(const XMMReg8Float &other)
1845 : {
1846 : ymm = other.ymm;
1847 : return *this;
1848 : }
1849 :
1850 : inline XMMReg8Float &operator+=(const XMMReg8Float &other)
1851 : {
1852 : ymm = _mm256_add_ps(ymm, other.ymm);
1853 : return *this;
1854 : }
1855 :
1856 : inline XMMReg8Float &operator-=(const XMMReg8Float &other)
1857 : {
1858 : ymm = _mm256_sub_ps(ymm, other.ymm);
1859 : return *this;
1860 : }
1861 :
1862 : inline XMMReg8Float &operator*=(const XMMReg8Float &other)
1863 : {
1864 : ymm = _mm256_mul_ps(ymm, other.ymm);
1865 : return *this;
1866 : }
1867 :
1868 : inline XMMReg8Float operator+(const XMMReg8Float &other) const
1869 : {
1870 : XMMReg8Float ret;
1871 : ret.ymm = _mm256_add_ps(ymm, other.ymm);
1872 : return ret;
1873 : }
1874 :
1875 : inline XMMReg8Float operator-(const XMMReg8Float &other) const
1876 : {
1877 : XMMReg8Float ret;
1878 : ret.ymm = _mm256_sub_ps(ymm, other.ymm);
1879 : return ret;
1880 : }
1881 :
1882 : inline XMMReg8Float operator*(const XMMReg8Float &other) const
1883 : {
1884 : XMMReg8Float ret;
1885 : ret.ymm = _mm256_mul_ps(ymm, other.ymm);
1886 : return ret;
1887 : }
1888 :
1889 : inline XMMReg8Float operator/(const XMMReg8Float &other) const
1890 : {
1891 : XMMReg8Float ret;
1892 : ret.ymm = _mm256_div_ps(ymm, other.ymm);
1893 : return ret;
1894 : }
1895 :
1896 : inline XMMReg8Float inverse() const
1897 : {
1898 : XMMReg8Float ret;
1899 : ret.ymm = _mm256_div_ps(_mm256_set1_ps(1.0f), ymm);
1900 : return ret;
1901 : }
1902 :
1903 : inline XMMReg8Float approx_inv_sqrt(const XMMReg8Float &one,
1904 : const XMMReg8Float &half) const
1905 : {
1906 : __m256 reg = ymm;
1907 : __m256 reg_half = _mm256_mul_ps(reg, half.ymm);
1908 : // Compute rough approximation of 1 / sqrt(b) with _mm256_rsqrt_ps
1909 : reg = _mm256_rsqrt_ps(reg);
1910 : // And perform one step of Newton-Raphson approximation to improve it
1911 : // approx_inv_sqrt_x = approx_inv_sqrt_x*(1.5 -
1912 : // 0.5*x*approx_inv_sqrt_x*approx_inv_sqrt_x);
1913 : const __m256 one_and_a_half = _mm256_add_ps(one.ymm, half.ymm);
1914 : reg = _mm256_mul_ps(
1915 : reg,
1916 : _mm256_sub_ps(one_and_a_half,
1917 : _mm256_mul_ps(reg_half, _mm256_mul_ps(reg, reg))));
1918 : XMMReg8Float ret;
1919 : ret.ymm = reg;
1920 : return ret;
1921 : }
1922 :
1923 : inline void StoreAllVal(float *ptr) const
1924 : {
1925 : Store8Val(ptr);
1926 : }
1927 :
1928 : inline void Store8Val(float *ptr) const
1929 : {
1930 : _mm256_storeu_ps(ptr, ymm);
1931 : }
1932 :
1933 : XMMReg8Float cast_to_float() const
1934 : {
1935 : return *this;
1936 : }
1937 : };
1938 :
1939 : #if defined(__AVX2__)
1940 :
1941 : class XMMReg8Int
1942 : {
1943 : public:
1944 : __m256i ymm;
1945 :
1946 : XMMReg8Int()
1947 : #if !defined(_MSC_VER)
1948 : : ymm(_mm256_undefined_si256())
1949 : #endif
1950 : {
1951 : }
1952 :
1953 : XMMReg8Int(const XMMReg8Int &other) : ymm(other.ymm)
1954 : {
1955 : }
1956 :
1957 : inline XMMReg8Int &operator=(const XMMReg8Int &other)
1958 : {
1959 : ymm = other.ymm;
1960 : return *this;
1961 : }
1962 :
1963 : static inline XMMReg8Int Zero()
1964 : {
1965 : XMMReg8Int reg;
1966 : reg.ymm = _mm256_setzero_si256();
1967 : return reg;
1968 : }
1969 :
1970 : static inline XMMReg8Int Set1(int i)
1971 : {
1972 : XMMReg8Int reg;
1973 : reg.ymm = _mm256_set1_epi32(i);
1974 : return reg;
1975 : }
1976 :
1977 : static inline XMMReg8Int LoadAllVal(const int *ptr)
1978 : {
1979 : return Load8Val(ptr);
1980 : }
1981 :
1982 : static inline XMMReg8Int Load8Val(const int *ptr)
1983 : {
1984 : XMMReg8Int reg;
1985 : reg.nsLoad8Val(ptr);
1986 : return reg;
1987 : }
1988 :
1989 : inline void nsLoad8Val(const int *ptr)
1990 : {
1991 : ymm = _mm256_loadu_si256(reinterpret_cast<const __m256i *>(ptr));
1992 : }
1993 :
1994 : static inline XMMReg8Int Equals(const XMMReg8Int &expr1,
1995 : const XMMReg8Int &expr2)
1996 : {
1997 : XMMReg8Int reg;
1998 : reg.ymm = _mm256_cmpeq_epi32(expr1.ymm, expr2.ymm);
1999 : return reg;
2000 : }
2001 :
2002 : static inline XMMReg8Int Ternary(const XMMReg8Int &cond,
2003 : const XMMReg8Int &true_expr,
2004 : const XMMReg8Int &false_expr)
2005 : {
2006 : XMMReg8Int reg;
2007 : reg.ymm =
2008 : _mm256_or_si256(_mm256_and_si256(cond.ymm, true_expr.ymm),
2009 : _mm256_andnot_si256(cond.ymm, false_expr.ymm));
2010 : return reg;
2011 : }
2012 :
2013 : inline XMMReg8Int &operator+=(const XMMReg8Int &other)
2014 : {
2015 : ymm = _mm256_add_epi32(ymm, other.ymm);
2016 : return *this;
2017 : }
2018 :
2019 : inline XMMReg8Int &operator-=(const XMMReg8Int &other)
2020 : {
2021 : ymm = _mm256_sub_epi32(ymm, other.ymm);
2022 : return *this;
2023 : }
2024 :
2025 : inline XMMReg8Int operator+(const XMMReg8Int &other) const
2026 : {
2027 : XMMReg8Int ret;
2028 : ret.ymm = _mm256_add_epi32(ymm, other.ymm);
2029 : return ret;
2030 : }
2031 :
2032 : inline XMMReg8Int operator-(const XMMReg8Int &other) const
2033 : {
2034 : XMMReg8Int ret;
2035 : ret.ymm = _mm256_sub_epi32(ymm, other.ymm);
2036 : return ret;
2037 : }
2038 :
2039 : XMMReg8Float cast_to_float() const
2040 : {
2041 : XMMReg8Float ret;
2042 : ret.ymm = _mm256_cvtepi32_ps(ymm);
2043 : return ret;
2044 : }
2045 : };
2046 :
2047 : #endif
2048 :
2049 : #else
2050 :
2051 : class XMMReg4Double
2052 : {
2053 : public:
2054 : XMMReg2Double low, high;
2055 :
2056 1723620054 : XMMReg4Double() : low(XMMReg2Double()), high(XMMReg2Double())
2057 : {
2058 1722050054 : }
2059 :
2060 1344000 : XMMReg4Double(const XMMReg4Double &other) : low(other.low), high(other.high)
2061 : {
2062 1344600 : }
2063 :
2064 253001000 : static inline XMMReg4Double Zero()
2065 : {
2066 253001000 : XMMReg4Double reg;
2067 252970000 : reg.low.Zeroize();
2068 252800000 : reg.high.Zeroize();
2069 252956000 : return reg;
2070 : }
2071 :
2072 111 : static inline XMMReg4Double Set1(double d)
2073 : {
2074 111 : XMMReg4Double reg;
2075 111 : reg.low = XMMReg2Double::Set1(d);
2076 111 : reg.high = XMMReg2Double::Set1(d);
2077 111 : return reg;
2078 : }
2079 :
2080 121304002 : static inline XMMReg4Double Load1ValHighAndLow(const double *ptr)
2081 : {
2082 121304002 : XMMReg4Double reg;
2083 121356002 : reg.low.nsLoad1ValHighAndLow(ptr);
2084 121456002 : reg.high = reg.low;
2085 121304002 : return reg;
2086 : }
2087 :
2088 371272002 : static inline XMMReg4Double Load4Val(const unsigned char *ptr)
2089 : {
2090 371272002 : XMMReg4Double reg;
2091 371007002 : XMMReg2Double::Load4Val(ptr, reg.low, reg.high);
2092 372368002 : return reg;
2093 : }
2094 :
2095 512 : static inline void Load8Val(const unsigned char *ptr, XMMReg4Double &low,
2096 : XMMReg4Double &high)
2097 : {
2098 512 : low = Load4Val(ptr);
2099 512 : high = Load4Val(ptr + 4);
2100 512 : }
2101 :
2102 1014862 : static inline XMMReg4Double Load4Val(const short *ptr)
2103 : {
2104 1014862 : XMMReg4Double reg;
2105 1014862 : reg.low.nsLoad2Val(ptr);
2106 1014862 : reg.high.nsLoad2Val(ptr + 2);
2107 1014862 : return reg;
2108 : }
2109 :
2110 512 : static inline void Load8Val(const short *ptr, XMMReg4Double &low,
2111 : XMMReg4Double &high)
2112 : {
2113 512 : low = Load4Val(ptr);
2114 512 : high = Load4Val(ptr + 4);
2115 512 : }
2116 :
2117 22641102 : static inline XMMReg4Double Load4Val(const unsigned short *ptr)
2118 : {
2119 22641102 : XMMReg4Double reg;
2120 22554202 : reg.low.nsLoad2Val(ptr);
2121 22665002 : reg.high.nsLoad2Val(ptr + 2);
2122 22729202 : return reg;
2123 : }
2124 :
2125 512 : static inline void Load8Val(const unsigned short *ptr, XMMReg4Double &low,
2126 : XMMReg4Double &high)
2127 : {
2128 512 : low = Load4Val(ptr);
2129 512 : high = Load4Val(ptr + 4);
2130 512 : }
2131 :
2132 1024 : static inline XMMReg4Double Load4Val(const int *ptr)
2133 : {
2134 1024 : XMMReg4Double reg;
2135 1024 : reg.low.nsLoad2Val(ptr);
2136 1024 : reg.high.nsLoad2Val(ptr + 2);
2137 1024 : return reg;
2138 : }
2139 :
2140 512 : static inline void Load8Val(const int *ptr, XMMReg4Double &low,
2141 : XMMReg4Double &high)
2142 : {
2143 512 : low = Load4Val(ptr);
2144 512 : high = Load4Val(ptr + 4);
2145 512 : }
2146 :
2147 269856016 : static inline XMMReg4Double Load4Val(const double *ptr)
2148 : {
2149 269856016 : XMMReg4Double reg;
2150 270613016 : reg.low.nsLoad2Val(ptr);
2151 270524016 : reg.high.nsLoad2Val(ptr + 2);
2152 271305016 : return reg;
2153 : }
2154 :
2155 944 : static inline void Load8Val(const double *ptr, XMMReg4Double &low,
2156 : XMMReg4Double &high)
2157 : {
2158 944 : low = Load4Val(ptr);
2159 944 : high = Load4Val(ptr + 4);
2160 944 : }
2161 :
2162 106840000 : static inline XMMReg4Double Load4ValAligned(const double *ptr)
2163 : {
2164 106840000 : XMMReg4Double reg;
2165 106786000 : reg.low.nsLoad2ValAligned(ptr);
2166 106717000 : reg.high.nsLoad2ValAligned(ptr + 2);
2167 106820000 : return reg;
2168 : }
2169 :
2170 38658 : static inline XMMReg4Double Load4Val(const float *ptr)
2171 : {
2172 38658 : XMMReg4Double reg;
2173 38658 : XMMReg2Double::Load4Val(ptr, reg.low, reg.high);
2174 38658 : return reg;
2175 : }
2176 :
2177 512 : static inline void Load8Val(const float *ptr, XMMReg4Double &low,
2178 : XMMReg4Double &high)
2179 : {
2180 512 : low = Load4Val(ptr);
2181 512 : high = Load4Val(ptr + 4);
2182 512 : }
2183 :
2184 2 : static inline XMMReg4Double Equals(const XMMReg4Double &expr1,
2185 : const XMMReg4Double &expr2)
2186 : {
2187 2 : XMMReg4Double reg;
2188 2 : reg.low = XMMReg2Double::Equals(expr1.low, expr2.low);
2189 2 : reg.high = XMMReg2Double::Equals(expr1.high, expr2.high);
2190 2 : return reg;
2191 : }
2192 :
2193 1336732 : static inline XMMReg4Double NotEquals(const XMMReg4Double &expr1,
2194 : const XMMReg4Double &expr2)
2195 : {
2196 1336732 : XMMReg4Double reg;
2197 1338782 : reg.low = XMMReg2Double::NotEquals(expr1.low, expr2.low);
2198 1336522 : reg.high = XMMReg2Double::NotEquals(expr1.high, expr2.high);
2199 1336442 : return reg;
2200 : }
2201 :
2202 6 : static inline XMMReg4Double Greater(const XMMReg4Double &expr1,
2203 : const XMMReg4Double &expr2)
2204 : {
2205 6 : XMMReg4Double reg;
2206 6 : reg.low = XMMReg2Double::Greater(expr1.low, expr2.low);
2207 6 : reg.high = XMMReg2Double::Greater(expr1.high, expr2.high);
2208 6 : return reg;
2209 : }
2210 :
2211 1333960 : static inline XMMReg4Double And(const XMMReg4Double &expr1,
2212 : const XMMReg4Double &expr2)
2213 : {
2214 1333960 : XMMReg4Double reg;
2215 1338550 : reg.low = XMMReg2Double::And(expr1.low, expr2.low);
2216 1338480 : reg.high = XMMReg2Double::And(expr1.high, expr2.high);
2217 1338840 : return reg;
2218 : }
2219 :
2220 2 : static inline XMMReg4Double Ternary(const XMMReg4Double &cond,
2221 : const XMMReg4Double &true_expr,
2222 : const XMMReg4Double &false_expr)
2223 : {
2224 2 : XMMReg4Double reg;
2225 : reg.low =
2226 2 : XMMReg2Double::Ternary(cond.low, true_expr.low, false_expr.low);
2227 : reg.high =
2228 2 : XMMReg2Double::Ternary(cond.high, true_expr.high, false_expr.high);
2229 2 : return reg;
2230 : }
2231 :
2232 4234382 : static inline XMMReg4Double Min(const XMMReg4Double &expr1,
2233 : const XMMReg4Double &expr2)
2234 : {
2235 4234382 : XMMReg4Double reg;
2236 4238692 : reg.low = XMMReg2Double::Min(expr1.low, expr2.low);
2237 4203722 : reg.high = XMMReg2Double::Min(expr1.high, expr2.high);
2238 4231262 : return reg;
2239 : }
2240 :
2241 143412008 : inline XMMReg4Double &operator=(const XMMReg4Double &other)
2242 : {
2243 143412008 : low = other.low;
2244 143466008 : high = other.high;
2245 143340008 : return *this;
2246 : }
2247 :
2248 585750002 : inline XMMReg4Double &operator+=(const XMMReg4Double &other)
2249 : {
2250 585750002 : low += other.low;
2251 588046002 : high += other.high;
2252 590735002 : return *this;
2253 : }
2254 :
2255 11357002 : inline XMMReg4Double &operator*=(const XMMReg4Double &other)
2256 : {
2257 11357002 : low *= other.low;
2258 11365402 : high *= other.high;
2259 11365302 : return *this;
2260 : }
2261 :
2262 8 : inline XMMReg4Double operator+(const XMMReg4Double &other) const
2263 : {
2264 8 : XMMReg4Double ret;
2265 8 : ret.low = low + other.low;
2266 8 : ret.high = high + other.high;
2267 8 : return ret;
2268 : }
2269 :
2270 2 : inline XMMReg4Double operator-(const XMMReg4Double &other) const
2271 : {
2272 2 : XMMReg4Double ret;
2273 2 : ret.low = low - other.low;
2274 2 : ret.high = high - other.high;
2275 2 : return ret;
2276 : }
2277 :
2278 620387002 : inline XMMReg4Double operator*(const XMMReg4Double &other) const
2279 : {
2280 620387002 : XMMReg4Double ret;
2281 619584002 : ret.low = low * other.low;
2282 618505002 : ret.high = high * other.high;
2283 617653002 : return ret;
2284 : }
2285 :
2286 1338232 : inline XMMReg4Double operator/(const XMMReg4Double &other) const
2287 : {
2288 1338232 : XMMReg4Double ret;
2289 1339052 : ret.low = low / other.low;
2290 1338702 : ret.high = high / other.high;
2291 1338432 : return ret;
2292 : }
2293 :
2294 582967 : void AddToLow(const XMMReg2Double &other)
2295 : {
2296 582967 : low += other;
2297 582966 : }
2298 :
2299 242647002 : inline double GetHorizSum() const
2300 : {
2301 242647002 : return (low + high).GetHorizSum();
2302 : }
2303 :
2304 : #if !defined(USE_SSE2_EMULATION)
2305 : inline XMMReg4Double approx_inv_sqrt(const XMMReg4Double &one,
2306 : const XMMReg4Double &half) const
2307 : {
2308 : __m128d reg0 = low.xmm;
2309 : __m128d reg1 = high.xmm;
2310 : __m128d reg0_half = _mm_mul_pd(reg0, half.low.xmm);
2311 : __m128d reg1_half = _mm_mul_pd(reg1, half.low.xmm);
2312 : // Compute rough approximation of 1 / sqrt(b) with _mm_rsqrt_ps
2313 : reg0 = _mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(reg0)));
2314 : reg1 = _mm_cvtps_pd(_mm_rsqrt_ps(_mm_cvtpd_ps(reg1)));
2315 : // And perform one step of Newton-Raphson approximation to improve it
2316 : // approx_inv_sqrt_x = approx_inv_sqrt_x*(1.5 -
2317 : // 0.5*x*approx_inv_sqrt_x*approx_inv_sqrt_x);
2318 : const __m128d one_and_a_half = _mm_add_pd(one.low.xmm, half.low.xmm);
2319 : reg0 = _mm_mul_pd(
2320 : reg0, _mm_sub_pd(one_and_a_half,
2321 : _mm_mul_pd(reg0_half, _mm_mul_pd(reg0, reg0))));
2322 : reg1 = _mm_mul_pd(
2323 : reg1, _mm_sub_pd(one_and_a_half,
2324 : _mm_mul_pd(reg1_half, _mm_mul_pd(reg1, reg1))));
2325 : XMMReg4Double ret;
2326 : ret.low.xmm = reg0;
2327 : ret.high.xmm = reg1;
2328 : return ret;
2329 : }
2330 :
2331 : inline XMMReg4Float cast_to_float() const
2332 : {
2333 : XMMReg4Float ret;
2334 : ret.xmm = _mm_castsi128_ps(
2335 : _mm_unpacklo_epi64(_mm_castps_si128(_mm_cvtpd_ps(low.xmm)),
2336 : _mm_castps_si128(_mm_cvtpd_ps(high.xmm))));
2337 : return ret;
2338 : }
2339 : #endif
2340 :
2341 1672502 : inline void Store4Val(unsigned char *ptr) const
2342 : {
2343 : #ifdef USE_SSE2_EMULATION
2344 1 : low.Store2Val(ptr);
2345 1 : high.Store2Val(ptr + 2);
2346 : #else
2347 3338672 : __m128i tmpLow = _mm_cvttpd_epi32(_mm_add_pd(
2348 1672501 : low.xmm,
2349 : _mm_set1_pd(0.5))); /* Convert the 2 double values to 2 integers */
2350 3334302 : __m128i tmpHigh = _mm_cvttpd_epi32(_mm_add_pd(
2351 1666171 : high.xmm,
2352 : _mm_set1_pd(0.5))); /* Convert the 2 double values to 2 integers */
2353 5011683 : auto tmp = _mm_castps_si128(_mm_shuffle_ps(_mm_castsi128_ps(tmpLow),
2354 : _mm_castsi128_ps(tmpHigh),
2355 : _MM_SHUFFLE(1, 0, 1, 0)));
2356 1674591 : tmp = _mm_packs_epi32(tmp, tmp);
2357 1675511 : tmp = _mm_packus_epi16(tmp, tmp);
2358 1675511 : GDALCopyXMMToInt32(tmp, reinterpret_cast<GInt32 *>(ptr));
2359 : #endif
2360 1664442 : }
2361 :
2362 2570502 : inline void Store4Val(unsigned short *ptr) const
2363 : {
2364 : #if 1
2365 2570502 : low.Store2Val(ptr);
2366 2568312 : high.Store2Val(ptr + 2);
2367 : #else
2368 : __m128i xmm0 = _mm_cvtpd_epi32(low.xmm);
2369 : __m128i xmm1 = _mm_cvtpd_epi32(high.xmm);
2370 : xmm0 = _mm_or_si128(xmm0, _mm_slli_si128(xmm1, 8));
2371 : #if __SSE4_1__
2372 : xmm0 = _mm_packus_epi32(xmm0, xmm0); // Pack uint32 to uint16
2373 : #else
2374 : xmm0 = _mm_add_epi32(xmm0, _mm_set1_epi32(-32768));
2375 : xmm0 = _mm_packs_epi32(xmm0, xmm0);
2376 : xmm0 = _mm_sub_epi16(xmm0, _mm_set1_epi16(-32768));
2377 : #endif
2378 : GDALCopyXMMToInt64(xmm0, (GInt64 *)ptr);
2379 : #endif
2380 2585382 : }
2381 :
2382 50900902 : inline void Store4Val(float *ptr) const
2383 : {
2384 50900902 : low.Store2Val(ptr);
2385 50913102 : high.Store2Val(ptr + 2);
2386 50938802 : }
2387 :
2388 3616 : inline void Store4Val(double *ptr) const
2389 : {
2390 3616 : low.Store2Val(ptr);
2391 3616 : high.Store2Val(ptr + 2);
2392 3616 : }
2393 :
2394 8 : inline void StoreMask(unsigned char *ptr) const
2395 : {
2396 8 : low.StoreMask(ptr);
2397 8 : high.StoreMask(ptr + 16);
2398 8 : }
2399 : };
2400 :
2401 : #if !defined(USE_SSE2_EMULATION)
2402 : inline XMMReg4Double XMMReg4Float::cast_to_double() const
2403 : {
2404 : XMMReg4Double ret;
2405 : ret.low.xmm = _mm_cvtps_pd(xmm);
2406 : ret.high.xmm = _mm_cvtps_pd(
2407 : _mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(xmm), 8)));
2408 : return ret;
2409 : }
2410 :
2411 : inline XMMReg4Double XMMReg4Int::cast_to_double() const
2412 : {
2413 : XMMReg4Double ret;
2414 : ret.low.xmm = _mm_cvtepi32_pd(xmm);
2415 : ret.high.xmm = _mm_cvtepi32_pd(_mm_srli_si128(xmm, 8));
2416 : return ret;
2417 : }
2418 : #endif
2419 :
2420 : #endif
2421 :
2422 : #endif /* #ifndef DOXYGEN_SKIP */
2423 :
2424 : #endif /* GDALSSE_PRIV_H_INCLUDED */
|
