EIC code displayed by LXR master/​include/​flint/​machine_vectors.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-11-07 10:15:27

 /*
     Copyright (C) 2022 Daniel Schultz
     Copyright (C) 2023 Mathieu Gouttenoire
 
     This file is part of FLINT.
 
     FLINT is free software: you can redistribute it and/or modify it under
     the terms of the GNU Lesser General Public License (LGPL) as published
     by the Free Software Foundation; either version 3 of the License, or
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
 #ifndef MACHINE_VECTORS_H
 #define MACHINE_VECTORS_H
 
 #define ALIGN_STRUCT(x) __attribute__((aligned(x)))
 
 #include <math.h>
 
 #if defined(__GNUC__)
 # if defined(__AVX2__)
 #  include <immintrin.h>
 # elif defined(__ARM_NEON)
 #  include <arm_neon.h>
 # endif
 #elif defined(_MSC_VER)
 # if defined(__AVX2__)
 #  include <intrin.h>
 # elif defined(_M_ARM64)
 #  include <arm_neon.h>
 # endif
 #endif
 
 #include "flint.h"
 #include "templates.h"
 
 #ifdef __cplusplus
 extern "C" {
 #endif
 
 #if defined(__AVX2__)
 /*
     In general the machine vector types should either be passed by const ref or
     the whole function should be forced inline as some platforms have buggy
     pass by value.
 */
 
 typedef ulong vec1n;
 typedef __m128i vec2n;
 typedef __m256i vec4n;
 typedef struct {__m256i e1, e2;} vec8n;
 
 typedef double vec1d;
 typedef __m128d vec2d;
 typedef __m256d vec4d;
 typedef struct {__m256d e1, e2;} vec8d;
 
 
 
 /* vec1d -- AVX2 ***********************************************************/
 
 FLINT_FORCE_INLINE vec1d vec1d_load(const double* a) {
     return a[0];
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_load_aligned(const double* a) {
     return a[0];
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_load_unaligned(const double* a) {
     return a[0];
 }
 
 FLINT_FORCE_INLINE void vec1d_store(double* z, vec1d a) {
     z[0] = a;
 }
 
 FLINT_FORCE_INLINE void vec1d_store_aligned(double* z, vec1d a) {
     z[0] = a;
 }
 
 FLINT_FORCE_INLINE void vec1d_store_unaligned(double* z, vec1d a) {
     z[0] = a;
 }
 
 FLINT_FORCE_INLINE int vec1d_same(double a, double b) {
     return a == b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_set_d(double a) {
     return a;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_round(vec1d a) {
     return rint(a);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_zero(void) {
     return 0.0;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_one(void) {
     return 1.0;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_add(vec1d a, vec1d b) {
     return a + b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_sub(vec1d a, vec1d b) {
     return a - b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_addsub(vec1d a, vec1d b) {
     return a - b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_neg(vec1d a) {
     return -a;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_abs(vec1d a) {
     return fabs(a);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_max(vec1d a, vec1d b) {
     return fmax(a, b);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_min(vec1d a, vec1d b) {
     return fmin(a, b);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_mul(vec1d a, vec1d b) {
     return a*b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_half(vec1d a) {
     return a*0.5;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_div(vec1d a, vec1d b) {
     return a/b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fmadd(vec1d a, vec1d b, vec1d c) {
     return fma(a, b, c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fmsub(vec1d a, vec1d b, vec1d c) {
     return fma(a, b, -c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fnmadd(vec1d a, vec1d b, vec1d c) {
     return fma(-a, b, c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fnmsub(vec1d a, vec1d b, vec1d c) {
     return fma(-a, b, -c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_blendv(vec1d a, vec1d b, vec1d c) {
     return c >= 0 ? a : b;
 }
 
 /* [0,n] -> [-n/2, n/2] */
 FLINT_FORCE_INLINE vec1d vec1d_reduce_0n_to_pmhn(vec1d a, vec1d n) {
     vec1d halfn = 0.5*n;
     return a > halfn ? a - n : a;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_reduce_pm1n_to_pmhn(vec1d a, vec1d n) {
     vec1d t = a + n;
     vec1d halfn = 0.5*n;
     if (a > halfn)
         return a - n;
     else if (t < halfn)
         return t;
     else
         return a;
 }
 
 /* [0,2n) to [0,n) */
 FLINT_FORCE_INLINE vec1d vec1d_reduce_2n_to_n(vec1d a, vec1d n) {
     return a-n >= 0 ? a-n : a;
 }
 
 
 /* vec4d -- AVX2 ***********************************************************/
 
 FLINT_FORCE_INLINE void vec4d_print(vec4d a)
 {
 #ifdef _MSC_VER
     double as[4];
     _mm256_storeu_pd(as, a);
     flint_printf("{%f, %f, %f, %f}", as[0], as[1], as[2], as[3]);
 #else
     flint_printf("{%f, %f, %f, %f}", a[0], a[1], a[2], a[3]);
 #endif
 }
 
 FLINT_FORCE_INLINE double vec4d_get_index(vec4d a, const int i) {
 #ifdef _MSC_VER
     double as[4];
     _mm256_storeu_pd(as, a);
     return as[i];
 #else
     return a[i];
 #endif
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_load(const double* a) {
     return _mm256_load_pd(a);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_load_aligned(const double* a) {
     return _mm256_load_pd(a);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_load_unaligned(const double* a) {
     return _mm256_loadu_pd(a);
 }
 
 FLINT_FORCE_INLINE void vec4d_store(double* z, vec4d a) {
     _mm256_store_pd(z, a);
 }
 
 FLINT_FORCE_INLINE void vec4d_store_aligned(double* z, vec4d a) {
     _mm256_store_pd(z, a);
 }
 
 FLINT_FORCE_INLINE void vec4d_store_unaligned(double* z, vec4d a) {
     _mm256_storeu_pd(z, a);
 }
 
 FLINT_FORCE_INLINE int vec4d_same(vec4d a, vec4d b) {
 #ifdef _MSC_VER
     double as[4], bs[4];
     _mm256_storeu_pd(as, a);
     _mm256_storeu_pd(bs, b);
     return as[0] == bs[0] && as[1] == bs[1] && as[2] == bs[2] && as[3] == bs[3];
 #else
     return a[0] == b[0] && a[1] == b[1] && a[2] == b[2] && a[3] == b[3];
 #endif
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_set_d(double a) {
     return _mm256_set1_pd(a);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_set_d4(double a0, double a1, double a2, double a3) {
     return _mm256_set_pd(a3, a2, a1, a0);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_round(vec4d a) {
     return _mm256_round_pd(a, 4);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_zero(void) {
     return _mm256_setzero_pd();
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_one(void) {
     return vec4d_set_d(1);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_add(vec4d a, vec4d b) {
     return _mm256_add_pd(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_sub(vec4d a, vec4d b) {
     return _mm256_sub_pd(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_addsub(vec4d a, vec4d b) {
     return _mm256_addsub_pd(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_neg(vec4d a) {
     __m256d mask = _mm256_castsi256_pd(_mm256_set1_epi64x(UWORD(0x8000000000000000)));
     return _mm256_xor_pd(a, mask);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_abs(vec4d a) {
     __m256d mask = _mm256_castsi256_pd(_mm256_set1_epi64x(UWORD(0x7fffffffffffffff)));
     return _mm256_and_pd(a, mask);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_max(vec4d a, vec4d b) {
     return _mm256_max_pd(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_min(vec4d a, vec4d b) {
     return _mm256_min_pd(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_mul(vec4d a, vec4d b) {
     return _mm256_mul_pd(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_half(vec4d a) {
     return vec4d_mul(a, vec4d_set_d(0.5));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_div(vec4d a, vec4d b) {
     return _mm256_div_pd(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_fmadd(vec4d a, vec4d b, vec4d c) {
 #ifndef AVOID_AVX2
     return _mm256_fmadd_pd(a, b, c);
 #else
     return _mm256_macc_pd(a, b, c);
 #endif
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_fmsub(vec4d a, vec4d b, vec4d c) {
 #ifndef AVOID_AVX2
     return _mm256_fmsub_pd(a, b, c);
 #else
     return _mm256_msub_pd(a, b, c);
 #endif
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_fnmadd(vec4d a, vec4d b, vec4d c) {
 #ifndef AVOID_AVX2
     return _mm256_fnmadd_pd(a, b, c);
 #else
     return _mm256_nmacc_pd(a, b, c);
 #endif
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_fnmsub(vec4d a, vec4d b, vec4d c) {
 #ifndef AVOID_AVX2
     return _mm256_fnmsub_pd(a, b, c);
 #else
     return _mm256_nmsub_pd(a, b, c);
 #endif
 }
 
 
 FLINT_FORCE_INLINE vec4d vec4d_blendv(vec4d a, vec4d b, vec4d c) {
     return _mm256_blendv_pd(a, b, c);
 }
 
 // return {a[0], b[0], a[2], b[2]}
 FLINT_FORCE_INLINE vec4d vec4d_unpacklo(vec4d a, vec4d b) {
     return _mm256_unpacklo_pd(a, b);
 }
 
 // return {a[1], b[1], a[3], b[3]}
 FLINT_FORCE_INLINE vec4d vec4d_unpackhi(vec4d a, vec4d b) {
     return _mm256_unpackhi_pd(a, b);
 }
 
 /* permute_i0_i1_i2_i3(a): return {a[i0], a[i1], a[i2], a[i3]} */
 #ifndef AVOID_AVX2
 #define DEFINE_IT(i0, i1, i2, i3) \
 FLINT_FORCE_INLINE vec4d CAT6(vec4d, permute, i0, i1, i2, i3)(vec4d a) { \
     return _mm256_permute4x64_pd(a, i0 + 4*(i1 + 4*(i2 + 4*i3))); \
 }
 #else
 #define DEFINE_IT(i0, i1, i2, i3) \
 FLINT_FORCE_INLINE vec4d CAT6(vec4d, permute, i0, i1, i2, i3)(vec4d a) { \
     return vec4d_set_d4(a[i0], a[i1], a[i2], a[i3]); \
 }
 #endif
 DEFINE_IT(0,2,1,3)
 DEFINE_IT(3,1,2,0)
 DEFINE_IT(3,2,1,0)
 #undef DEFINE_IT
 
 /* permute2_i0_i1(a):  return {v[i0], v[i1]}
                     |   v[0]     |    v[1]    |    v[2]    |   v[3]     |
                       a[0], a[1]   a[2], a[3]   b[0], b[1]   b[2], b[3]
 */
 #define DEFINE_IT(i0, i1) \
 FLINT_FORCE_INLINE vec4d CAT4(vec4d, permute2, i0, i1)(vec4d a, vec4d b) { \
     return _mm256_permute2f128_pd(a, b, i0 + 16*i1); \
 }
 DEFINE_IT(0,2)
 DEFINE_IT(1,3)
 #undef DEFINE_IT
 
 FLINT_FORCE_INLINE vec4d vec4d_unpack_lo_permute_0_2_1_3(vec4d u, vec4d v) {
     return vec4d_permute_0_2_1_3(vec4d_unpacklo(u, v));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_unpack_hi_permute_0_2_1_3(vec4d u, vec4d v) {
     return vec4d_permute_0_2_1_3(vec4d_unpackhi(u, v));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_unpackhi_permute_3_1_2_0(vec4d u, vec4d v) {
     return vec4d_permute_3_1_2_0(vec4d_unpackhi(u, v));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_unpacklo_permute_3_1_2_0(vec4d u, vec4d v) {
     return vec4d_permute_3_1_2_0(vec4d_unpacklo(u, v));
 }
 
 /* view the 4 vectors as the rows of a 4x4 matrix */
 #define VEC4D_TRANSPOSE(z0, z1, z2, z3, a0, a1, a2, a3) \
 { \
     vec4d _t0, _t1, _t2, _t3; \
     _t0 = vec4d_unpacklo(a0, a1); \
     _t1 = vec4d_unpackhi(a0, a1); \
     _t2 = vec4d_unpacklo(a2, a3); \
     _t3 = vec4d_unpackhi(a2, a3); \
     z0 = vec4d_permute2_0_2(_t0, _t2); \
     z1 = vec4d_permute2_0_2(_t1, _t3); \
     z2 = vec4d_permute2_1_3(_t0, _t2); \
     z3 = vec4d_permute2_1_3(_t1, _t3); \
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_cmp_ge(vec4d a, vec4d b) {
     return _mm256_cmp_pd(a, b, _CMP_GE_OQ);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_cmp_gt(vec4d a, vec4d b) {
     return _mm256_cmp_pd(a, b, _CMP_GT_OQ);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_reduce_0n_to_pmhn(vec4d a, vec4d n) {
     vec4d halfn = vec4d_half(n);
     return vec4d_blendv(a, vec4d_sub(a, n), vec4d_cmp_gt(a, halfn));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_reduce_pm1n_to_pmhn(vec4d a, vec4d n) {
     vec4d halfn = vec4d_half(n);
     vec4d t = vec4d_blendv(n, vec4d_neg(n), a);
     a = vec4d_blendv(a, vec4d_sub(a, t), vec4d_cmp_gt(vec4d_abs(a), halfn));
     return a;
 }
 
 /* [0,2n) to [0,n) */
 FLINT_FORCE_INLINE vec4d vec4d_reduce_2n_to_n(vec4d a, vec4d n) {
     vec4d s = vec4d_sub(a, n);
     return vec4d_blendv(s, a, s);
 }
 
 FLINT_FORCE_INLINE vec4n vec4d_convert_limited_vec4n(vec4d a) {
     __m256d t = _mm256_set1_pd(0x1.0p52);
     return _mm256_castpd_si256(_mm256_xor_pd(_mm256_add_pd(a, t), t));
 }
 
 /* vec4n -- AVX2 ***********************************************************/
 
 FLINT_FORCE_INLINE void vec4n_print(vec4n a)
 {
     flint_printf("[hi %016llx_%016llx_%016llx_%016llx lo]",
         _mm256_extract_epi64(a, 3),
         _mm256_extract_epi64(a, 2),
         _mm256_extract_epi64(a, 1),
         _mm256_extract_epi64(a, 0));
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_load_unaligned(const ulong* a) {
     return _mm256_loadu_si256((__m256i*) a);
 }
 
 FLINT_FORCE_INLINE void vec4n_store_unaligned(ulong* z, vec4n a) {
     _mm256_storeu_si256((__m256i*) z, a);
 }
 
 /* permute_i0_i1_i2_i3(a): return {a[i0], a[i1], a[i2], a[i3]} */
 #ifndef AVOID_AVX2
 #define DEFINE_IT(i0, i1, i2, i3)                                        \
 FLINT_FORCE_INLINE vec4n CAT6(vec4n, permute, i0, i1, i2, i3)(vec4n a) { \
     return _mm256_permute4x64_epi64(a, i0 + 4*(i1 + 4*(i2 + 4*i3)));     \
 }
 #else
 #define DEFINE_IT(i0, i1, i2, i3)                                        \
 FLINT_FORCE_INLINE vec4n CAT6(vec4n, permute, i0, i1, i2, i3)(vec4n a) { \
     return vec4n_set_n4(a[i0], a[i1], a[i2], a[i3]);                     \
 }
 #endif
 DEFINE_IT(3,2,1,0)
 #undef DEFINE_IT
 
 FLINT_FORCE_INLINE vec4n vec4n_zero(void)
 {
     return _mm256_setzero_si256();
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_set_n4(ulong a0, ulong a1, ulong a2, ulong a3) {
     return _mm256_set_epi64x(a3, a2, a1, a0);
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_set_n(ulong a) {
   return _mm256_set1_epi64x(a);
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_add(vec4n a, vec4n b)
 {
     return _mm256_add_epi64(a, b);
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_sub(vec4n a, vec4n b)
 {
     return _mm256_sub_epi64(a, b);
 }
 
 /* for n < 2^63 */
 FLINT_FORCE_INLINE vec4n vec4n_addmod_limited(vec4n a, vec4n b, vec4n n)
 {
     vec4n s = vec4n_add(a, b);
     vec4n t = vec4n_sub(s, n);
     vec4n m = _mm256_srai_epi32(t, 31);
           m = _mm256_shuffle_epi32(m, 1 + 4*(1 + 4*(3 + 4*(3))));
     return _mm256_blendv_epi8(t, s, m);
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_addmod(vec4n a, vec4n b, vec4n n)
 {
     vec4n tt = vec4n_set_n(0x8000000000000000);
     vec4n s = vec4n_add(a, b);
 #if 0
     vec4n  m = vec4n_sub(n, tt);
     vec4n t0 = vec4n_sub(m, a);
     vec4n t1 = vec4n_sub(b, tt);
     vec4n t2 = vec4n_sub(t1, t0);
     return _mm256_blendv_epi8(s, t2, _mm256_cmpgt_epi64(t1, t0));
 #else
     vec4n t0 = vec4n_sub(s, n);
     vec4n t1 = vec4n_sub(a, tt);
     vec4n t2 = vec4n_sub(t0, tt);
     return _mm256_blendv_epi8(t0, s, _mm256_cmpgt_epi64(t2, t1));
 #endif
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_mul(vec4n u, vec4n v)
 {
     return _mm256_mul_epu32(u, v);
 }
 
 // compilers accept non-(compile-time-)constant b with srli
 // (still, faster if constant)
 FLINT_FORCE_INLINE vec4n vec4n_bit_shift_right(vec4n a, ulong b) {
     return _mm256_srli_epi64(a, b);
 }
 
 #define vec4n_bit_shift_right_32(a) vec4n_bit_shift_right((a), 32)
 
 FLINT_FORCE_INLINE vec4n vec4n_bit_and(vec4n a, vec4n b) {
     return _mm256_and_si256(a, b);
 }
 
 FLINT_FORCE_INLINE vec4d vec4n_convert_limited_vec4d(vec4n a) {
     __m256d t = _mm256_set1_pd(0x1.0p52);
     return _mm256_sub_pd(_mm256_or_pd(_mm256_castsi256_pd(a), t), t);
 }
 
 // horizontal sum
 FLINT_FORCE_INLINE ulong vec4n_horizontal_sum(vec4n a) {
     vec4n a_hi = _mm256_shuffle_epi32(a, 14);  // 14 == 0b00001110
     vec4n sum_lo = _mm256_add_epi64(a, a_hi);
     vec2n sum_hi = _mm256_extracti128_si256(sum_lo, 1);
     vec2n sum = _mm_add_epi64(_mm256_castsi256_si128(sum_lo), sum_hi);
     return (ulong) _mm_cvtsi128_si64(sum);
 }
 
 
 /* vec8d -- AVX2 ***********************************************************/
 
 FLINT_FORCE_INLINE double vec8d_get_index(vec8d a, int i) {
     return i < 4 ? vec4d_get_index(a.e1, i) : vec4d_get_index(a.e2, i - 4);
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_set_d(double a) {
     vec4d z1 = vec4d_set_d(a);
     vec8d z = {z1, z1};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_set_d8(double a0, double a1, double a2, double a3, double a4, double a5, double a6, double a7) {
     vec4d z1 = vec4d_set_d4(a0, a1, a2, a3);
     vec4d z2 = vec4d_set_d4(a4, a5, a6, a7);
     vec8d z = {z1, z2};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_load(const double* a) {
     vec8d z = {vec4d_load(a+0), vec4d_load(a+4)};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_load_aligned(const double* a) {
     vec8d z = {vec4d_load_aligned(a+0), vec4d_load_aligned(a+4)};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_load_unaligned(const double* a) {
     vec8d z = {vec4d_load_unaligned(a+0), vec4d_load_unaligned(a+4)};
     return z;
 }
 
 FLINT_FORCE_INLINE void vec8d_store(double* z, vec8d a) {
     vec4d_store(z+0, a.e1);
     vec4d_store(z+4, a.e2);
 }
 
 FLINT_FORCE_INLINE void vec8d_store_aligned(double* z, vec8d a) {
     vec4d_store_aligned(z+0, a.e1);
     vec4d_store_aligned(z+4, a.e2);
 }
 
 FLINT_FORCE_INLINE void vec8d_store_unaligned(double* z, vec8d a) {
     vec4d_store_unaligned(z+0, a.e1);
     vec4d_store_unaligned(z+4, a.e2);
 }
 
 FLINT_FORCE_INLINE int vec8d_same(vec8d a, vec8d b) {
     return vec4d_same(a.e1, b.e1) && vec4d_same(a.e2, b.e2);
 }
 
 FLINT_FORCE_INLINE vec8d vec8n_convert_limited_vec8d(vec8n a) {
     vec8d z = {vec4n_convert_limited_vec4d(a.e1), vec4n_convert_limited_vec4d(a.e2)};
     return z;
 }
 
 /* !!! the outputs are also permuted !!! */
 FLINT_FORCE_INLINE vec8d _vec8i32_convert_vec8d(__m256i a)
 {
     __m256i mask = _mm256_set1_epi32(0x43300000);
     __m256i ak0 = _mm256_unpacklo_epi32(a, mask);
     __m256i ak1 = _mm256_unpackhi_epi32(a, mask);
     __m256d t = _mm256_set1_pd(0x1.0p52);
     vec8d z;
     z.e1 = _mm256_sub_pd(_mm256_castsi256_pd(ak0), t);
     z.e2 = _mm256_sub_pd(_mm256_castsi256_pd(ak1), t);
     return z;
 }
 
 
 /* reduce_pm1no_to_0n(a, n): return a mod n in [0,n) assuming a in (-n,n) */
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE V V##_reduce_pm1no_to_0n(V a, V n) { \
     return V##_blendv(a, V##_add(a, n), a); \
 }
 DEFINE_IT(vec1d)
 DEFINE_IT(vec4d)
 #undef DEFINE_IT
 
 /* reduce_to_pm1n(a, n, ninv): return a mod n in [-n,n] */
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE V V##_reduce_to_pm1n(V a, V n, V ninv) { \
     return V##_fnmadd(V##_round(V##_mul(a, ninv)), n, a); \
 }
 DEFINE_IT(vec1d)
 DEFINE_IT(vec4d)
 #undef DEFINE_IT
 
 /* reduce_to_pm1n(a, n, ninv): return a mod n in (-n,n) */
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE V V##_reduce_to_pm1no(V a, V n, V ninv) { \
     return V##_fnmadd(V##_round(V##_mul(a, ninv)), n, a); \
 }
 DEFINE_IT(vec1d)
 DEFINE_IT(vec4d)
 #undef DEFINE_IT
 
 
 /* reduce_to_0n(a, n, ninv): return a mod n in [0,n) */
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE V V##_reduce_to_0n(V a, V n, V ninv) { \
     return V##_reduce_pm1no_to_0n(V##_reduce_to_pm1no(a, n, ninv), n); \
 }
 DEFINE_IT(vec1d)
 DEFINE_IT(vec4d)
 #undef DEFINE_IT
 
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE int V##_same_mod(V a, V b, V n, V ninv) { \
     return V##_same(V##_reduce_to_0n(a, n, ninv), V##_reduce_to_0n(b, n, ninv)); \
 }
 DEFINE_IT(vec1d)
 DEFINE_IT(vec4d)
 #undef DEFINE_IT
 
 /*
     mulmod(a, b, n, 1.0/n): return a*b mod n
     assuming fft_small_mulmod_satisfies_bounds(n):
      Products a*b in the range (-2*n^2, 2*n^2) are reduced to the range (-1*n, 1*n).
      Products a*b in the range (-4*n^2, 4*n^2) are reduced to the range (-3/2*n, 3/2*n).
     otherwise the more pessimistic:
      Products a*b in the range (-2*n^2, 2*n^2) are reduced to the range (-9/8*n, 9/8*n).
      Products a*b in the range (-4*n^2, 4*n^2) are reduced to the range (-7/4*n, 7/4*n).
     floating point round mode is always assumed to be set to nearest
 */
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE V V##_mulmod(V a, V b, V n, V ninv) { \
     V h = V##_mul(a, b); \
     V q = V##_round(V##_mul(h, ninv)); \
     V l = V##_fmsub(a, b, h); \
     return V##_add(V##_fnmadd(q, n, h), l); \
 } \
  \
 FLINT_FORCE_INLINE V V##_nmulmod(V a, V b, V n, V ninv) { \
     V h = V##_mul(a, b); \
     V q = V##_round(V##_mul(h, ninv)); \
     V l = V##_fnmadd(a, b, h); \
     return V##_sub(l, V##_fnmadd(q, n, h)); \
 }
 
 DEFINE_IT(vec1d)
 DEFINE_IT(vec4d)
 #undef DEFINE_IT
 
 
 
 
 #define EXTEND_VEC_DEF0(U, V, f) \
 FLINT_FORCE_INLINE V V##f(void) { \
     U z1 = U##f(); \
     U z2 = U##f(); \
     V z = {z1, z2}; \
     return z; \
 }
 
 #define EXTEND_VEC_DEF1(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a) { \
     U z1 = U##f(a.e1); \
     U z2 = U##f(a.e2); \
     V z = {z1, z2}; \
     return z; \
 }
 
 #define EXTEND_VEC_DEF2(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a, V b) { \
     U z1 = U##f(a.e1, b.e1); \
     U z2 = U##f(a.e2, b.e2); \
     V z = {z1, z2}; \
     return z; \
 }
 
 #define EXTEND_VEC_DEF3(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a, V b, V c) { \
     U z1 = U##f(a.e1, b.e1, c.e1); \
     U z2 = U##f(a.e2, b.e2, c.e2); \
     V z = {z1, z2}; \
     return z; \
 }
 
 #define EXTEND_VEC_DEF4(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a, V b, V c, V d) { \
     U z1 = U##f(a.e1, b.e1, c.e1, d.e1); \
     U z2 = U##f(a.e2, b.e2, c.e2, d.e2); \
     V z = {z1, z2}; \
     return z; \
 }
 
 EXTEND_VEC_DEF0(vec4d, vec8d, _zero)
 EXTEND_VEC_DEF1(vec4d, vec8d, _neg)
 EXTEND_VEC_DEF1(vec4d, vec8d, _round)
 EXTEND_VEC_DEF2(vec4d, vec8d, _add)
 EXTEND_VEC_DEF2(vec4d, vec8d, _sub)
 EXTEND_VEC_DEF2(vec4d, vec8d, _min)
 EXTEND_VEC_DEF2(vec4d, vec8d, _max)
 EXTEND_VEC_DEF2(vec4d, vec8d, _mul)
 EXTEND_VEC_DEF2(vec4d, vec8d, _div)
 EXTEND_VEC_DEF2(vec4d, vec8d, _reduce_pm1n_to_pmhn)
 EXTEND_VEC_DEF2(vec4d, vec8d, _reduce_pm1no_to_0n)
 EXTEND_VEC_DEF2(vec4d, vec8d, _reduce_2n_to_n)
 EXTEND_VEC_DEF3(vec4d, vec8d, _reduce_to_pm1n)
 EXTEND_VEC_DEF3(vec4d, vec8d, _reduce_to_pm1no)
 EXTEND_VEC_DEF3(vec4d, vec8d, _reduce_to_0n)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fmadd)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fmsub)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fnmadd)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fnmsub)
 EXTEND_VEC_DEF3(vec4d, vec8d, _blendv)
 EXTEND_VEC_DEF3(vec4n, vec8n, _addmod)
 EXTEND_VEC_DEF3(vec4n, vec8n, _addmod_limited)
 EXTEND_VEC_DEF4(vec4d, vec8d, _mulmod)
 EXTEND_VEC_DEF4(vec4d, vec8d, _nmulmod)
 
 #undef EXTEND_VEC_DEF4
 #undef EXTEND_VEC_DEF3
 #undef EXTEND_VEC_DEF2
 #undef EXTEND_VEC_DEF1
 #undef EXTEND_VEC_DEF0
 
 
 
 /* vec8n -- AVX2 ***********************************************************/
 
 FLINT_FORCE_INLINE vec8n vec8n_load_unaligned(const ulong* a) {
     vec8n z = {vec4n_load_unaligned(a+0), vec4n_load_unaligned(a+4)};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8n vec8n_set_n(ulong a) {
     vec4n x = vec4n_set_n(a);
     vec8n z = {x, x};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8n vec8n_bit_shift_right(vec8n a, ulong b) {
     vec8n z = {vec4n_bit_shift_right(a.e1, b), vec4n_bit_shift_right(a.e2, b)};
     return z;
 }
 
 #define vec8n_bit_shift_right_32(a) vec8n_bit_shift_right((a), 32)
 
 FLINT_FORCE_INLINE vec8n vec8n_bit_and(vec8n a, vec8n b) {
     vec8n z = {vec4n_bit_and(a.e1, b.e1), vec4n_bit_and(a.e2, b.e2)};
     return z;
 }
 
 
 
 #elif defined(__ARM_NEON) || defined(_M_ARM64)
 
 typedef ulong vec1n;
 typedef uint64x2_t vec2n;
 typedef struct {vec2n e1, e2;} vec4n;
 typedef struct {vec4n e1, e2;} vec8n;
 
 typedef double vec1d;
 typedef float64x2_t vec2d;
 typedef struct {vec2d e1, e2;} vec4d;
 typedef struct {vec4d e1, e2;} vec8d;
 
 #define EXTEND_VEC_DEF0(U, V, f) \
 FLINT_FORCE_INLINE V V##f(void) { \
     U z1 = U##f(); \
     U z2 = U##f(); \
     V z = {z1, z2}; return z; \
 }
 
 #define EXTEND_VEC_DEF1(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a) { \
     U z1 = U##f(a.e1); \
     U z2 = U##f(a.e2); \
     V z = {z1, z2}; return z; \
 }
 
 #define EXTEND_VEC_DEF2(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a, V b) { \
     U z1 = U##f(a.e1, b.e1); \
     U z2 = U##f(a.e2, b.e2); \
     V z = {z1, z2}; return z; \
 }
 
 #define EXTEND_VEC_DEF3(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a, V b, V c) { \
     U z1 = U##f(a.e1, b.e1, c.e1); \
     U z2 = U##f(a.e2, b.e2, c.e2); \
     V z = {z1, z2}; return z; \
 }
 
 #define EXTEND_VEC_DEF4(U, V, f) \
 FLINT_FORCE_INLINE V V##f(V a, V b, V c, V d) { \
     U z1 = U##f(a.e1, b.e1, c.e1, d.e1); \
     U z2 = U##f(a.e2, b.e2, c.e2, d.e2); \
     V z = {z1, z2}; return z; \
 }
 
 
 /* floating point stuff ******************************************************/
 
 /* vec1d -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE vec1d vec1d_load(const double* a) {
     return a[0];
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_load_aligned(const double* a) {
     return a[0];
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_load_unaligned(const double* a) {
     return a[0];
 }
 
 FLINT_FORCE_INLINE void vec1d_store(double* z, vec1d a) {
     z[0] = a;
 }
 
 FLINT_FORCE_INLINE void vec1d_store_aligned(double* z, vec1d a) {
     z[0] = a;
 }
 
 FLINT_FORCE_INLINE void vec1d_store_unaligned(double* z, vec1d a) {
     z[0] = a;
 }
 
 FLINT_FORCE_INLINE int vec1d_same(double a, double b) {
     return a == b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_set_d(double a) {
     return a;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_round(vec1d a) {
     return rint(a);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_zero(void) {
     return 0.0;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_one(void) {
     return 1.0;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_add(vec1d a, vec1d b) {
     return a + b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_sub(vec1d a, vec1d b) {
     return a - b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_addsub(vec1d a, vec1d b) {
     return a - b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_neg(vec1d a) {
     return -a;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_abs(vec1d a) {
     return fabs(a);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_max(vec1d a, vec1d b) {
     return fmax(a, b);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_min(vec1d a, vec1d b) {
     return fmin(a, b);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_mul(vec1d a, vec1d b) {
     return a*b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_half(vec1d a) {
     return a*0.5;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_div(vec1d a, vec1d b) {
     return a/b;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fmadd(vec1d a, vec1d b, vec1d c) {
     return fma(a, b, c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fmsub(vec1d a, vec1d b, vec1d c) {
     return fma(a, b, -c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fnmadd(vec1d a, vec1d b, vec1d c) {
     return fma(-a, b, c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_fnmsub(vec1d a, vec1d b, vec1d c) {
     return fma(-a, b, -c);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_blendv(vec1d a, vec1d b, vec1d c) {
     return c >= 0 ? a : b;
 }
 
 /* [0,n] -> [-n/2, n/2] */
 FLINT_FORCE_INLINE vec1d vec1d_reduce_0n_to_pmhn(vec1d a, vec1d n) {
     vec1d halfn = 0.5*n;
     return a > halfn ? a - n : a;
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_reduce_pm1n_to_pmhn(vec1d a, vec1d n) {
     vec1d t = a + n;
     vec1d halfn = 0.5*n;
     if (a > halfn)
         return a - n;
     else if (t < halfn)
         return t;
     else
         return a;
 }
 
 /* vec2d -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE double vec2d_get_index(vec2d a, int i) {
     return a[i];
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_set_d2(double a0, double a1)
 {
     double ALIGN_STRUCT(16) data[2] = {a0, a1};
     return vld1q_f64((float64_t *) data);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_set_d(double a)
 {
     return vdupq_n_f64(a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_load(const double* a)
 {
     return vld1q_f64(a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_load_unaligned(const double* a)
 {
     return vld1q_f64(a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_load_aligned(const double* a)
 {
     return vld1q_f64(a);
 }
 
 FLINT_FORCE_INLINE void vec2d_store(double* z, vec2d a)
 {
     vst1q_f64(z, a);
 }
 
 FLINT_FORCE_INLINE void vec2d_store_unaligned(double* z, vec2d a)
 {
     vst1q_f64(z, a);
 }
 
 FLINT_FORCE_INLINE void vec2d_store_aligned(double* z, vec2d a)
 {
     vst1q_f64(z, a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_zero(void)
 {
     return vdupq_n_f64(0.0);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_one(void)
 {
     return vdupq_n_f64(1.0);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_neg(vec2d a)
 {
     return vnegq_f64(a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_abs(vec2d a)
 {
     return vabsq_f64(a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_min(vec2d a, vec2d b)
 {
     return vminq_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_max(vec2d a, vec2d b)
 {
     return vmaxq_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_cmp_lt(vec2d a, vec2d b) {
     return vcvtq_f64_u64(vcltq_f64(a, b));
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_cmp_gt(vec2d a, vec2d b) {
     return vcvtq_f64_u64(vcgtq_f64(a, b));
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_blendv(vec2d a, vec2d b, vec2d c) {
     return vbslq_f64(vcvtq_u64_f64(c), b, a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_add(vec2d a, vec2d b)
 {
     return vaddq_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_sub(vec2d a, vec2d b)
 {
     return vsubq_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_mul(vec2d a, vec2d b)
 {
     return vmulq_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_half(vec2d a)
 {
     return vmulq_n_f64(a, 0.5);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_fmadd(vec2d a, vec2d b, vec2d c) {
     return vfmaq_f64(c, a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_fnmadd(vec2d a, vec2d b, vec2d c) {
     return vfmsq_f64(c, a, b);
 }
 
 /* these two might be slow due to the extra negation */
 FLINT_FORCE_INLINE vec2d vec2d_fmsub(vec2d a, vec2d b, vec2d c) {
     return -vfmsq_f64(c, a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_fnmsub(vec2d a, vec2d b, vec2d c) {
     return -vfmaq_f64(c, a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_div(vec2d a, vec2d b)
 {
     return vdivq_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_round(vec2d a)
 {
     return vrndnq_f64(a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_unpacklo(vec2d a, vec2d b) {
     return vtrn1q_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_unpackhi(vec2d a, vec2d b) {
     return vtrn2q_f64(a, b);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_reduce_to_pm1no(vec1d a, vec1d n, vec1d ninv) {
     return vec1d_fnmadd(vec1d_round(vec1d_mul(a, ninv)), n, a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_reduce_to_pm1no(vec2d a, vec2d n, vec2d ninv) {
     return vec2d_fnmadd(vec2d_round(vec2d_mul(a, ninv)), n, a);
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_reduce_to_pm1n(vec1d a, vec1d n, vec1d ninv) {
     return vec1d_fnmadd(vec1d_round(vec1d_mul(a, ninv)), n, a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_reduce_to_pm1n(vec2d a, vec2d n, vec2d ninv) {
     return vec2d_fnmadd(vec2d_round(vec2d_mul(a, ninv)), n, a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_reduce_0n_to_pmhn(vec2d a, vec2d n) {
     vec2d halfn = vec2d_half(n);
     return vbslq_f64(vcgtq_f64(a, halfn), vec2d_sub(a, n), a);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_reduce_pm1n_to_pmhn(vec2d a, vec2d n) {
     vec2d halfn = vec2d_half(n);
     vec2d t = vec2d_add(a, n);
 
     vec2n condition_a = vcgtq_f64(a, halfn);
     vec2n condition_t = vcltq_f64(t, halfn);
 
     return vbslq_f64(condition_a, vec2d_sub(a, n), vbslq_f64(condition_t, t, a));
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_reduce_pm1no_to_0n(vec1d a, vec1d n) {
     return a >= 0 ? a : a + n;
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_reduce_pm1no_to_0n(vec2d a, vec2d n) {
     return vbslq_f64(vcgeq_f64(a, vec2d_zero()), a, vaddq_f64(a, n));
 }
 
 FLINT_FORCE_INLINE vec1d vec1d_reduce_to_0n(vec1d a, vec1d n, vec1d ninv) {
     return vec1d_reduce_pm1no_to_0n(vec1d_reduce_to_pm1no(a, n, ninv), n);
 }
 
 FLINT_FORCE_INLINE vec2d vec2d_reduce_to_0n(vec2d a, vec2d n, vec2d ninv) {
     return vec2d_reduce_pm1no_to_0n(vec2d_reduce_to_pm1no(a, n, ninv), n);
 }
 
 
 /* mulmod(a, b, n, ninv): return a*b mod n in [-n,n] with assumptions */
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE V V##_mulmod(V a, V b, V n, V ninv) { \
     V h = V##_mul(a, b); \
     V q = V##_round(V##_mul(h, ninv)); \
     V l = V##_fnmadd(a, b, h); \
     return V##_sub(V##_fnmadd(q, n, h), l); \
 } \
  \
 FLINT_FORCE_INLINE V V##_nmulmod(V a, V b, V n, V ninv) { \
     V h = V##_mul(a, b); \
     V q = V##_round(V##_mul(h, ninv)); \
     V l = V##_fnmadd(a, b, h); \
     return V##_sub(l, V##_fnmadd(q, n, h)); \
 }
 
 DEFINE_IT(vec1d)
 DEFINE_IT(vec2d)
 #undef DEFINE_IT
 
 
 
 /* vec4d -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE vec4d vec4d_set_vec2d2(vec2d a, vec2d b) {
     vec4d z = {a, b}; return z;
 }
 
 FLINT_FORCE_INLINE double vec4d_get_index(vec4d a, int i) {
     return i < 2 ? vec2d_get_index(a.e1, i) : vec2d_get_index(a.e2, i - 2);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_set_d4(double a0, double a1, double a2, double a3)
 {
     return vec4d_set_vec2d2(vec2d_set_d2(a0, a1), vec2d_set_d2(a2, a3));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_set_d(double a)
 {
     vec2d z1 = vec2d_set_d(a);
     vec4d z = {z1, z1}; return z;
 }
 
 FLINT_FORCE_INLINE void vec4d_store(double* z, vec4d a)
 {
     vec2d_store(z+0, a.e1);
     vec2d_store(z+2, a.e2);
 }
 
 FLINT_FORCE_INLINE void vec4d_store_aligned(double* z, vec4d a)
 {
     vec2d_store_aligned(z+0, a.e1);
     vec2d_store_aligned(z+2, a.e2);
 }
 
 FLINT_FORCE_INLINE void vec4d_store_unaligned(double* z, vec4d a)
 {
     vec2d_store_unaligned(z+0, a.e1);
     vec2d_store_unaligned(z+2, a.e2);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_load(const double* a)
 {
     return vec4d_set_vec2d2(vec2d_load(a+0),
                             vec2d_load(a+2));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_load_aligned(const double* a)
 {
     return vec4d_set_vec2d2(vec2d_load_aligned(a+0),
                             vec2d_load_aligned(a+2));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_load_unaligned(const double* a)
 {
     return vec4d_set_vec2d2(vec2d_load_unaligned(a+0),
                             vec2d_load_unaligned(a+2));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_permute_0_2_1_3(vec4d a)
 {
     return vec4d_set_vec2d2(vec2d_unpacklo(a.e1, a.e2),
                             vec2d_unpackhi(a.e1, a.e2));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_permute_3_2_1_0(vec4d a) {
     return vec4d_set_vec2d2(vextq_f64(a.e2, a.e2, 1),
                             vextq_f64(a.e1, a.e1, 1));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_permute_3_1_2_0(vec4d a) {
     vec2d z1 = vzip2q_f64(a.e2, a.e1);
     vec2d z2 = vzip1q_f64(a.e2, a.e1);
     return vec4d_set_vec2d2(z1, z2);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_permute2_0_2(vec4d a, vec4d b) {
     return vec4d_set_vec2d2(a.e1, b.e1);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_permute2_1_3(vec4d a, vec4d b) {
     return vec4d_set_vec2d2(a.e2, b.e2);
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_unpack_lo_permute_0_2_1_3(vec4d u, vec4d v) {
     return vec4d_set_vec2d2(vec2d_unpacklo(u.e1, u.e2),
                             vec2d_unpacklo(v.e1, v.e2));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_unpack_hi_permute_0_2_1_3(vec4d u, vec4d v) {
     return vec4d_set_vec2d2(vec2d_unpackhi(u.e1, u.e2),
                             vec2d_unpackhi(v.e1, v.e2));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_unpackhi_permute_3_1_2_0(vec4d u, vec4d v) {
     return vec4d_set_vec2d2(vec2d_unpackhi(v.e2, v.e1),
                             vec2d_unpackhi(u.e2, u.e1));
 }
 
 FLINT_FORCE_INLINE vec4d vec4d_unpacklo_permute_3_1_2_0(vec4d u, vec4d v) {
     return vec4d_set_vec2d2(vec2d_unpacklo(v.e2, v.e1),
                             vec2d_unpacklo(u.e2, u.e1));
 }
 
 #define VEC4D_TRANSPOSE(z0, z1, z2, z3, a0, a1, a2, a3) \
 { \
     vec4d _t0, _t1, _t2, _t3; \
     _t0 = vec4d_unpacklo(a0, a1); \
     _t1 = vec4d_unpackhi(a0, a1); \
     _t2 = vec4d_unpacklo(a2, a3); \
     _t3 = vec4d_unpackhi(a2, a3); \
     z0 = vec4d_permute2_0_2(_t0, _t2); \
     z1 = vec4d_permute2_0_2(_t1, _t3); \
     z2 = vec4d_permute2_1_3(_t0, _t2); \
     z3 = vec4d_permute2_1_3(_t1, _t3); \
 }
 
 EXTEND_VEC_DEF0(vec2d, vec4d, _zero)
 EXTEND_VEC_DEF0(vec2d, vec4d, _one)
 EXTEND_VEC_DEF1(vec2d, vec4d, _neg)
 EXTEND_VEC_DEF1(vec2d, vec4d, _round)
 EXTEND_VEC_DEF2(vec2d, vec4d, _add)
 EXTEND_VEC_DEF2(vec2d, vec4d, _sub)
 EXTEND_VEC_DEF2(vec2d, vec4d, _min)
 EXTEND_VEC_DEF2(vec2d, vec4d, _max)
 EXTEND_VEC_DEF2(vec2d, vec4d, _mul)
 EXTEND_VEC_DEF2(vec2d, vec4d, _div)
 EXTEND_VEC_DEF2(vec2d, vec4d, _reduce_pm1n_to_pmhn)
 EXTEND_VEC_DEF2(vec2d, vec4d, _reduce_pm1no_to_0n)
 EXTEND_VEC_DEF2(vec2d, vec4d, _unpacklo)
 EXTEND_VEC_DEF2(vec2d, vec4d, _unpackhi)
 EXTEND_VEC_DEF3(vec2d, vec4d, _reduce_to_pm1n)
 EXTEND_VEC_DEF3(vec2d, vec4d, _reduce_to_pm1no)
 EXTEND_VEC_DEF3(vec2d, vec4d, _reduce_to_0n)
 EXTEND_VEC_DEF3(vec2d, vec4d, _fmadd)
 EXTEND_VEC_DEF3(vec2d, vec4d, _fmsub)
 EXTEND_VEC_DEF3(vec2d, vec4d, _fnmadd)
 EXTEND_VEC_DEF3(vec2d, vec4d, _fnmsub)
 EXTEND_VEC_DEF3(vec2d, vec4d, _blendv)
 EXTEND_VEC_DEF4(vec2d, vec4d, _mulmod)
 EXTEND_VEC_DEF4(vec2d, vec4d, _nmulmod)
 
 
 /* vec8d -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE vec8d vec8d_set_vec4d2(vec4d a, vec4d b) {
     vec8d z = {a, b}; return z;
 }
 
 FLINT_FORCE_INLINE double vec8d_get_index(vec8d a, int i) {
     return i < 4 ? vec4d_get_index(a.e1, i) : vec4d_get_index(a.e2, i - 4);
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_set_d8(double a0, double a1, double a2, double a3,
                                       double a4, double a5, double a6, double a7)
 {
     vec4d z1 = vec4d_set_d4(a0, a1, a2, a3);
     vec4d z2 = vec4d_set_d4(a4, a5, a6, a7);
     vec8d z = {z1, z2}; return z;
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_set_d(double a)
 {
     vec4d z1 = vec4d_set_d(a);
     vec8d z = {z1, z1}; return z;
 }
 
 FLINT_FORCE_INLINE void vec8d_store(double* z, vec8d a)
 {
     vec4d_store(z+0, a.e1);
     vec4d_store(z+4, a.e2);
 }
 
 FLINT_FORCE_INLINE void vec8d_store_aligned(double* z, vec8d a)
 {
     vec4d_store_aligned(z+0, a.e1);
     vec4d_store_aligned(z+4, a.e2);
 }
 
 FLINT_FORCE_INLINE void vec8d_store_unaligned(double* z, vec8d a)
 {
     vec4d_store_unaligned(z+0, a.e1);
     vec4d_store_unaligned(z+4, a.e2);
 }
 
 
 FLINT_FORCE_INLINE vec8d vec8d_load(const double* a)
 {
     vec4d z1 = vec4d_load(a+0);
     vec4d z2 = vec4d_load(a+4);
     vec8d z = {z1, z2}; return z;
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_load_aligned(const double* a)
 {
     vec4d z1 = vec4d_load_aligned(a+0);
     vec4d z2 = vec4d_load_aligned(a+4);
     vec8d z = {z1, z2}; return z;
 }
 
 FLINT_FORCE_INLINE vec8d vec8d_load_unaligned(const double* a)
 {
     vec4d z1 = vec4d_load_unaligned(a+0);
     vec4d z2 = vec4d_load_unaligned(a+4);
     vec8d z = {z1, z2}; return z;
 }
 
 
 EXTEND_VEC_DEF0(vec4d, vec8d, _zero)
 EXTEND_VEC_DEF0(vec4d, vec8d, _one)
 EXTEND_VEC_DEF1(vec4d, vec8d, _neg)
 EXTEND_VEC_DEF1(vec4d, vec8d, _round)
 EXTEND_VEC_DEF2(vec4d, vec8d, _add)
 EXTEND_VEC_DEF2(vec4d, vec8d, _sub)
 EXTEND_VEC_DEF2(vec4d, vec8d, _min)
 EXTEND_VEC_DEF2(vec4d, vec8d, _max)
 EXTEND_VEC_DEF2(vec4d, vec8d, _mul)
 EXTEND_VEC_DEF2(vec4d, vec8d, _div)
 EXTEND_VEC_DEF2(vec4d, vec8d, _reduce_pm1n_to_pmhn)
 EXTEND_VEC_DEF2(vec4d, vec8d, _reduce_pm1no_to_0n)
 EXTEND_VEC_DEF2(vec4d, vec8d, _unpacklo)
 EXTEND_VEC_DEF2(vec4d, vec8d, _unpackhi)
 EXTEND_VEC_DEF3(vec4d, vec8d, _reduce_to_pm1n)
 EXTEND_VEC_DEF3(vec4d, vec8d, _reduce_to_pm1no)
 EXTEND_VEC_DEF3(vec4d, vec8d, _reduce_to_0n)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fmadd)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fmsub)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fnmadd)
 EXTEND_VEC_DEF3(vec4d, vec8d, _fnmsub)
 EXTEND_VEC_DEF3(vec4d, vec8d, _blendv)
 EXTEND_VEC_DEF4(vec4d, vec8d, _mulmod)
 EXTEND_VEC_DEF4(vec4d, vec8d, _nmulmod)
 
 
 
 
 FLINT_FORCE_INLINE int vec2d_same(vec2d a, vec2d b) {
     return a[0] == b[0] && a[1] == b[1];
 }
 
 
 #define DEFINE_IT(V) \
 FLINT_FORCE_INLINE int V##_same_mod(V a, V b, V n, V ninv) { \
     return V##_same(V##_reduce_to_0n(a, n, ninv), V##_reduce_to_0n(b, n, ninv)); \
 }
 DEFINE_IT(vec1d)
 DEFINE_IT(vec2d)
 #undef DEFINE_IT
 
 /* integer stuff *************************************************************/
 
 /* vec1n -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE void vec1n_store_unaligned(ulong* z, vec1n a) {
     z[0] = a;
 }
 
 FLINT_FORCE_INLINE vec1n vec1d_convert_limited_vec1n(vec1d a) {
     return (slong)a;
 }
 
 // (a + b) % n
 FLINT_FORCE_INLINE vec1n vec1n_addmod(vec1n a, vec1n b, vec1n n) {
     vec1n nmb = n - b;
     return nmb > a ? a + b : a - nmb;
 }
 
 /* vec2n -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE vec2d vec2n_convert_limited_vec2d(vec2n a) {
     float64x2_t t = vdupq_n_f64(0x1.0p52);
     return vsubq_f64(vreinterpretq_f64_u64(vorrq_u64(a, vreinterpretq_u64_f64(t))), t);
 }
 
 FLINT_FORCE_INLINE void vec2n_store_unaligned(ulong* z, vec2n a) {
    vst1q_u64((uint64_t *) z, a);
 }
 
 FLINT_FORCE_INLINE vec2n vec2d_convert_limited_vec2n(vec2d a) {
     return vcvtnq_u64_f64(a);
 }
 
 FLINT_FORCE_INLINE vec2n vec2n_set_n(ulong a) {
     vec2n x = vdupq_n_u64(a);
     return x;
 }
 
 FLINT_FORCE_INLINE vec2n vec2n_load_unaligned(const ulong* a) {
     return vld1q_u64((const uint64_t *) a);
 }
 
 // Right shift 32bits
 FLINT_FORCE_INLINE vec2n vec2n_bit_shift_right_32(vec2n a) {
     return vshrq_n_u64(a, 32);
 }
 
 // AND operation
 FLINT_FORCE_INLINE vec2n vec2n_bit_and(vec2n a, vec2n b) {
     return vandq_u64(a, b);
 }
 
 // Addition
 FLINT_FORCE_INLINE vec2n vec2n_add(vec2n a, vec2n b) {
     return vaddq_u64(a, b);
 }
 
 // Substraction
 FLINT_FORCE_INLINE vec2n vec2n_sub(vec2n a, vec2n b) {
     return vsubq_u64(a, b);
 }
 
 // (a + b) % n
 FLINT_FORCE_INLINE vec2n vec2n_addmod(vec2n a, vec2n b, vec2n n) {
     vec2n nmb = vec2n_sub(n, b);
     vec2n sum = vec2n_sub(a, nmb);
 
     vec2n mask = vcgtq_u64(nmb, a);
 
     return vec2n_add(sum, vandq_u64(n, mask));
 }
 
 // (a + b) % n for n < 2^63
 FLINT_FORCE_INLINE vec2n vec2n_addmod_limited(vec2n a, vec2n b, vec2n n) {
     vec2n s = vec2n_add(a, b);
 
     vec2n mask = vcgeq_u64(s, n);
 
     return vec2n_sub(s, vandq_u64(n, mask));
 }
 
 /* vec4n -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE vec4d vec4n_convert_limited_vec4d(vec4n a) {
     vec2d z1 = vec2n_convert_limited_vec2d(a.e1);
     vec2d z2 = vec2n_convert_limited_vec2d(a.e2);
     vec4d z = {z1, z2}; return z;
 }
 
 FLINT_FORCE_INLINE void vec4n_store_unaligned(ulong* z, vec4n a) {
     vec2n_store_unaligned(z+0, a.e1);
     vec2n_store_unaligned(z+2, a.e2);
 }
 
 FLINT_FORCE_INLINE vec4n vec4d_convert_limited_vec4n(vec4d a) {
     vec2n z1 = vec2d_convert_limited_vec2n(a.e1);
     vec2n z2 = vec2d_convert_limited_vec2n(a.e2);
     vec4n z = {z1, z2}; return z;
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_set_n(ulong a) {
     vec2n x = vec2n_set_n(a);
     vec4n z = {x, x};
     return z;
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_load_unaligned(const ulong* a) {
     vec2n z1 = vec2n_load_unaligned(a+0);
     vec2n z2 = vec2n_load_unaligned(a+2);
     vec4n z = {z1, z2}; return z;
 }
 
 EXTEND_VEC_DEF1(vec2n, vec4n, _bit_shift_right_32)
 EXTEND_VEC_DEF2(vec2n, vec4n, _bit_and)
 EXTEND_VEC_DEF2(vec2n, vec4n, _add)
 EXTEND_VEC_DEF2(vec2n, vec4n, _sub)
 EXTEND_VEC_DEF3(vec2n, vec4n, _addmod)
 EXTEND_VEC_DEF3(vec2n, vec4n, _addmod_limited)
 
 /* vec8n -- NEON/ARM64 *********************************************/
 
 FLINT_FORCE_INLINE vec8d vec8n_convert_limited_vec8d(vec8n a) {
     vec4d z1 = vec4n_convert_limited_vec4d(a.e1);
     vec4d z2 = vec4n_convert_limited_vec4d(a.e2);
     vec8d z = {z1, z2}; return z;
 }
 
 FLINT_FORCE_INLINE vec8n vec8n_set_n(ulong a) {
     vec4n x = vec4n_set_n(a);
     vec8n z = {x, x};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8n vec8n_load_unaligned(const ulong* a)
 {
     vec4n z1 = vec4n_load_unaligned(a+0);
     vec4n z2 = vec4n_load_unaligned(a+4);
     vec8n z = {z1, z2}; return z;
 }
 
 
 EXTEND_VEC_DEF1(vec4n, vec8n, _bit_shift_right_32)
 EXTEND_VEC_DEF2(vec4n, vec8n, _bit_and)
 EXTEND_VEC_DEF2(vec4n, vec8n, _add)
 EXTEND_VEC_DEF2(vec4n, vec8n, _sub)
 EXTEND_VEC_DEF3(vec4n, vec8n, _addmod)
 EXTEND_VEC_DEF3(vec4n, vec8n, _addmod_limited)
 
 
 #if 0
 /*
     vshrq_n_u64(a, n) cannot be used because n must be a compile-time
     constant, and the compiler doesn't see that n is constant
     even if the function is forced inline.
     And vshlq_s64(a, vdupq_n_s64(-(slong) n)) cannot be used to emulate
     vshrq_n_u64(a, n) as it propagates the sign bit.
 */
 FLINT_FORCE_INLINE vec2n vec2n_bit_shift_right(vec2n a, ulong n)
 {
     return vshrq_n_u64(a, n);
 }
 
 FLINT_FORCE_INLINE vec4n vec4n_bit_shift_right(vec4n a, ulong n)
 {
     vec4n z = {vec2n_bit_shift_right(a.e1, n), vec2n_bit_shift_right(a.e2, n)};
     return z;
 }
 
 FLINT_FORCE_INLINE vec8n vec8n_bit_shift_right(vec8n a, ulong n)
 {
     vec8n z = {vec4n_bit_shift_right(a.e1, n), vec4n_bit_shift_right(a.e2, n)};
     return z;
 }
 #endif
 
 
 #undef EXTEND_VEC_DEF4
 #undef EXTEND_VEC_DEF3
 #undef EXTEND_VEC_DEF2
 #undef EXTEND_VEC_DEF1
 #undef EXTEND_VEC_DEF0
 
 
 #else
 
 #error machine_vector.h not implemented
 
 #endif
 
 
 #ifdef __cplusplus
 }
 #endif
 
 #endif

This page was automatically generated by the 2.3.7 LXR engine. The LXR team