Back to home page

EIC code displayed by LXR
 
 

    

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

0001 /*
0002     Copyright (C) 2010 William Hart
0003     Copyright (C) 2021 Fredrik Johansson
0004     Copyright (C) 2024 Vincent Neiger
0005 
0006     This file is part of FLINT.
0007 
0008     FLINT is free software: you can redistribute it and/or modify it under
0009     the terms of the GNU Lesser General Public License (LGPL) as published
0010     by the Free Software Foundation; either version 3 of the License, or
0011     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0012 */
0013 
0014 #ifndef NMOD_VEC_H
0015 #define NMOD_VEC_H
0016 
0017 #ifdef NMOD_VEC_INLINES_C
0018 #define NMOD_VEC_INLINE
0019 #else
0020 #define NMOD_VEC_INLINE static inline
0021 #endif
0022 
0023 #include "flint.h"
0024 #include "nmod.h"  // nmod_mul, nmod_fmma
0025 
0026 #ifdef __cplusplus
0027 extern "C" {
0028 #endif
0029 
0030 #define NMOD_VEC_NORM(vec, i)                   \
0031 do {                                            \
0032     while ((i) && vec[(i) - 1] == UWORD(0))     \
0033         (i)--;                                  \
0034 } while (0)
0035 
0036 NMOD_VEC_INLINE
0037 nn_ptr _nmod_vec_init(slong len)
0038 {
0039     return (nn_ptr) flint_malloc(len * sizeof(ulong));
0040 }
0041 
0042 NMOD_VEC_INLINE
0043 void _nmod_vec_clear(nn_ptr vec)
0044 {
0045    flint_free(vec);
0046 }
0047 
0048 void _nmod_vec_randtest(nn_ptr vec, flint_rand_t state, slong len, nmod_t mod);
0049 
0050 void _nmod_vec_rand(nn_ptr vec, flint_rand_t state, slong len, nmod_t mod);
0051 
0052 NMOD_VEC_INLINE
0053 void _nmod_vec_zero(nn_ptr vec, slong len)
0054 {
0055     slong i;
0056     for (i = 0; i < len; i++)
0057         vec[i] = 0;
0058 }
0059 
0060 flint_bitcnt_t _nmod_vec_max_bits(nn_srcptr vec, slong len);
0061 
0062 NMOD_VEC_INLINE
0063 void _nmod_vec_set(nn_ptr res, nn_srcptr vec, slong len)
0064 {
0065     slong i;
0066     for (i = 0; i < len; i++)
0067         res[i] = vec[i];
0068 }
0069 
0070 NMOD_VEC_INLINE
0071 void _nmod_vec_swap(nn_ptr a, nn_ptr b, slong length)
0072 {
0073     slong i;
0074     for (i = 0; i < length; i++)
0075     {
0076         ulong t = a[i];
0077         a[i] = b[i];
0078         b[i] = t;
0079     }
0080 }
0081 
0082 NMOD_VEC_INLINE
0083 int _nmod_vec_equal(nn_srcptr vec, nn_srcptr vec2, slong len)
0084 {
0085     slong i;
0086 
0087     for (i = 0; i < len; i++)
0088         if (vec[i] != vec2[i]) return 0;
0089 
0090     return 1;
0091 }
0092 
0093 NMOD_VEC_INLINE
0094 int _nmod_vec_is_zero(nn_srcptr vec, slong len)
0095 {
0096     slong i;
0097 
0098     for (i = 0; i < len; i++)
0099         if (vec[i] != 0) return 0;
0100 
0101     return 1;
0102 }
0103 
0104 /* some IO functions */
0105 #ifdef FLINT_HAVE_FILE
0106 int _nmod_vec_fprint_pretty(FILE * file, nn_srcptr vec, slong len, nmod_t mod);
0107 int _nmod_vec_fprint(FILE * f, nn_srcptr vec, slong len, nmod_t mod);
0108 #endif
0109 
0110 void _nmod_vec_print_pretty(nn_srcptr vec, slong len, nmod_t mod);
0111 int _nmod_vec_print(nn_srcptr vec, slong len, nmod_t mod);
0112 
0113 /* reduce, add, scalar mul */
0114 void _nmod_vec_reduce(nn_ptr res, nn_srcptr vec, slong len, nmod_t mod);
0115 
0116 void _nmod_vec_add(nn_ptr res, nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0117 void _nmod_vec_sub(nn_ptr res, nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0118 void _nmod_vec_neg(nn_ptr res, nn_srcptr vec, slong len, nmod_t mod);
0119 
0120 void _nmod_vec_scalar_mul_nmod(nn_ptr res, nn_srcptr vec, slong len, ulong c, nmod_t mod);
0121 void _nmod_vec_scalar_mul_nmod_redc(nn_ptr res, nn_srcptr vec, slong len, ulong c, nmod_t mod);
0122 void _nmod_vec_scalar_mul_nmod_generic(nn_ptr res, nn_srcptr vec, slong len, ulong c, nmod_t mod);
0123 void _nmod_vec_scalar_mul_nmod_shoup(nn_ptr res, nn_srcptr vec, slong len, ulong c, nmod_t mod);
0124 
0125 void _nmod_vec_scalar_addmul_nmod(nn_ptr res, nn_srcptr vec, slong len, ulong c, nmod_t mod);
0126 void _nmod_vec_scalar_addmul_nmod_generic(nn_ptr res, nn_srcptr vec, slong len, ulong c, nmod_t mod);
0127 void _nmod_vec_scalar_addmul_nmod_shoup(nn_ptr res, nn_srcptr vec, slong len, ulong c, nmod_t mod);
0128 
0129 #if FLINT_BITS == 64 && defined(__AVX2__)
0130 void _nmod_vec_nored_scalar_addmul_halflimb_avx2(nn_ptr res, nn_srcptr vec, slong len, ulong c);
0131 #endif
0132 
0133 NMOD_VEC_INLINE void
0134 _nmod_vec_nored_scalar_addmul_halflimb(nn_ptr res, nn_srcptr vec, slong len, ulong c)
0135 {
0136 #if FLINT_BITS == 64 && defined(__AVX2__)
0137     if (len >= 16)
0138     {
0139         _nmod_vec_nored_scalar_addmul_halflimb_avx2(res, vec, len, c);
0140         return;
0141     }
0142 #endif
0143 
0144     slong i;
0145     for (i = 0; i < len; i++)
0146         res[i] += vec[i] * c;
0147 }
0148 
0149 void _nmod_vec_nored_scalar_addmul_halflimb(nn_ptr res, nn_srcptr vec, slong len, ulong c);
0150 void _nmod_vec_nored_ll_scalar_addmul_halflimb(nn_ptr res, nn_srcptr vec, slong len, ulong c);
0151 void _nmod_vec_nored_ll_scalar_addmul(nn_ptr res, nn_srcptr vec, slong len, ulong c);
0152 void _nmod_vec_nored_lll_scalar_addmul(nn_ptr res, nn_srcptr vec, slong len, ulong c);
0153 
0154 /* invert each coefficient */
0155 void _nmod_vec_invert(nn_ptr res, nn_srcptr vec, ulong len, nmod_t mod);
0156 void _nmod_vec_invert_naive(nn_ptr res, nn_srcptr vec, ulong len, nmod_t mod);
0157 void _nmod_vec_invert_generic(nn_ptr res, nn_srcptr vec, ulong len, nmod_t mod);
0158 void _nmod_vec_invert_shoup(nn_ptr res, nn_srcptr vec, ulong len, nmod_t mod);
0159 
0160 /* ---- compute dot parameters ---- */
0161 
0162 // for _DOT2_SPLIT
0163 #if (FLINT_BITS == 64)
0164 #   define DOT_SPLIT_BITS 56
0165 #   define DOT_SPLIT_MASK UWORD(72057594037927935) // (1L << DOT_SPLIT_BITS) - 1
0166 #endif // FLINT_BITS == 64
0167 
0168 #define _FIXED_LEN_MOD_BOUNDS(fixedlen, onelimb_bnd, twolimb_bnd) \
0169         if (len == fixedlen)                                      \
0170         {                                                         \
0171             if (mod.n <= UWORD(onelimb_bnd))                      \
0172             {                                                     \
0173                 dot_params_t params = {_DOT1, UWORD(0)};          \
0174                 return params;                                    \
0175             }                                                     \
0176             if (mod.n <= UWORD(twolimb_bnd))                      \
0177             {                                                     \
0178                 dot_params_t params = {_DOT2, UWORD(0)};          \
0179                 return params;                                    \
0180             }                                                     \
0181             dot_params_t params = {_DOT3, UWORD(0)};              \
0182             return params;                                        \
0183         }
0184 
0185 FLINT_FORCE_INLINE dot_params_t _nmod_vec_dot_params(ulong len, nmod_t mod)
0186 {
0187     if (len == 0 || mod.n == 1)
0188     {
0189         dot_params_t params = {_DOT0, UWORD(0)};
0190         return params;
0191     }
0192     if ((mod.n & (mod.n - 1)) == 0)
0193     {
0194         dot_params_t params = {_DOT_POW2, UWORD(0)};
0195         return params;
0196     }
0197     // from here on len >= 1, n > 1 not power of 2
0198 
0199     // short dot products: we use only _DOT1, _DOT2, _DOT3 in that case
0200     if (len <= 11)
0201     {
0202 #if FLINT_BITS == 64
0203         // 64 bits:  k limbs  <=>  n <= ceil(2**(32*k) / sqrt(len))
0204         _FIXED_LEN_MOD_BOUNDS(11, 1294981365,  5561902608746059656);
0205         _FIXED_LEN_MOD_BOUNDS(10, 1358187914,  5833372668713515885);
0206         _FIXED_LEN_MOD_BOUNDS( 9, 1431655766,  6148914691236517206);
0207         _FIXED_LEN_MOD_BOUNDS( 8, 1518500250,  6521908912666391107);
0208         _FIXED_LEN_MOD_BOUNDS( 7, 1623345051,  6972213902555716131);
0209         _FIXED_LEN_MOD_BOUNDS( 6, 1753413057,  7530851732716320753);
0210         _FIXED_LEN_MOD_BOUNDS( 5, 1920767767,  8249634742471189718);
0211         _FIXED_LEN_MOD_BOUNDS( 4, 2147483648,  9223372036854775808);
0212         _FIXED_LEN_MOD_BOUNDS( 3, 2479700525, 10650232656628343402);
0213         _FIXED_LEN_MOD_BOUNDS( 2, 3037000500, 13043817825332782213);
0214 #else  // FLINT_BITS == 64
0215         // 32 bits: k limbs  <=>  n <= ceil(2**(16*k) / sqrt(len))
0216         _FIXED_LEN_MOD_BOUNDS(11, 19760, 1294981365);
0217         _FIXED_LEN_MOD_BOUNDS(10, 20725, 1358187914);
0218         _FIXED_LEN_MOD_BOUNDS( 9, 21846, 1431655766);
0219         _FIXED_LEN_MOD_BOUNDS( 8, 23171, 1518500250);
0220         _FIXED_LEN_MOD_BOUNDS( 7, 24771, 1623345051);
0221         _FIXED_LEN_MOD_BOUNDS( 6, 26755, 1753413057);
0222         _FIXED_LEN_MOD_BOUNDS( 5, 29309, 1920767767);
0223         _FIXED_LEN_MOD_BOUNDS( 4, 32768, 2147483648);
0224         _FIXED_LEN_MOD_BOUNDS( 3, 37838, 2479700525);
0225         _FIXED_LEN_MOD_BOUNDS( 2, 46341, 3037000500);
0226 #endif  // FLINT_BITS == 64
0227         // remains len == 1
0228         if (mod.n <= (UWORD(1) << FLINT_BITS / 2))
0229         {
0230             dot_params_t params = {_DOT1, UWORD(0)};
0231             return params;
0232         }
0233         dot_params_t params = {_DOT2, UWORD(0)};
0234         return params;
0235     }
0236 
0237     if (mod.n <= UWORD(1) << (FLINT_BITS / 2)) // implies <= 2 limbs
0238     {
0239         const ulong t0 = (mod.n - 1) * (mod.n - 1);
0240         ulong u1 = n_mulhi(t0, len);
0241         if (u1 == 0)  // 1 limb
0242         {
0243             dot_params_t params = {_DOT1, UWORD(0)};
0244             return params;
0245         }
0246 
0247         // u1 != 0 <=> 2 limbs
0248 #if (FLINT_BITS == 64) // _SPLIT: see end of file for these constraints
0249         if (mod.n <= UWORD(1515531528) && len <= WORD(380368697))
0250         {
0251             ulong pow2_precomp;
0252             NMOD_RED(pow2_precomp, (UWORD(1) << DOT_SPLIT_BITS), mod);
0253             dot_params_t params = {_DOT2_SPLIT, pow2_precomp};
0254             return params;
0255         }
0256 #endif // FLINT_BITS == 64
0257 #if (FLINT_BITS == 64) && defined(__AVX2__)
0258         ulong pow2_precomp;
0259         NMOD_RED(pow2_precomp, (UWORD(1) << DOT_SPLIT_BITS), mod);
0260         dot_params_t params = {_DOT2_HALF, pow2_precomp};
0261         return params;
0262 #else // (FLINT_BITS == 64) && defined(__AVX2__)
0263         dot_params_t params = {_DOT2_HALF, UWORD(0)};
0264         return params;
0265 #endif // (FLINT_BITS == 64) && defined(__AVX2__)
0266     }
0267     // from here on, mod.n > 2**(FLINT_BITS / 2)
0268     // --> unreduced dot cannot fit in 1 limb
0269 
0270     ulong t2, t1, t0, u1;
0271     umul_ppmm(t1, t0, mod.n - 1, mod.n - 1);
0272     umul_ppmm(t2, t1, t1, len);
0273     u1 = n_mulhi(t0, len);
0274     add_ssaaaa(t2, t1, t2, t1, UWORD(0), u1);
0275 
0276     if (t2 == 0) // 2 limbs
0277     {
0278         dot_params_t params = {_DOT2, UWORD(0)};
0279         return params;
0280     }
0281 
0282     // 3 limbs:
0283 #if (FLINT_BITS == 64)
0284     if (mod.n <= UWORD(6521908912666391107))  // room for accumulating 8 terms
0285 #else
0286     if (mod.n <= UWORD(1518500250))           // room for accumulating 8 terms
0287 #endif
0288     {
0289         dot_params_t params = {_DOT3_ACC, UWORD(0)};
0290         return params;
0291     }
0292 
0293     dot_params_t params = {_DOT3, UWORD(0)};
0294     return params;
0295 }
0296 
0297 #undef _FIXED_LEN_MOD_BOUNDS
0298 
0299 int _nmod_vec_dot_bound_limbs(slong len, nmod_t mod);
0300 int _nmod_vec_dot_bound_limbs_from_params(slong len, nmod_t mod, dot_params_t params);
0301 
0302 
0303 /* ------ dot product, specific algorithms ------ */
0304 
0305 /* vec1[i] * vec2[i] */
0306 ulong _nmod_vec_dot_pow2(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0307 ulong _nmod_vec_dot1(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0308 ulong _nmod_vec_dot2_half(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod, ulong pow2_precomp);
0309 ulong _nmod_vec_dot2(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0310 ulong _nmod_vec_dot3_acc(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0311 ulong _nmod_vec_dot3(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0312 #if FLINT_BITS == 64
0313 ulong _nmod_vec_dot2_split(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod, ulong pow2_precomp);
0314 #endif  // FLINT_BITS == 64
0315 
0316 /* vec1[i] * vec2[len-1-i] */
0317 ulong _nmod_vec_dot_pow2_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0318 ulong _nmod_vec_dot1_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0319 ulong _nmod_vec_dot2_half_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0320 ulong _nmod_vec_dot2_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0321 ulong _nmod_vec_dot3_acc_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0322 ulong _nmod_vec_dot3_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod);
0323 #if FLINT_BITS == 64
0324 ulong _nmod_vec_dot2_split_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod, ulong pow2_precomp);
0325 #endif  // FLINT_BITS == 64
0326 
0327 /* vec1[i] * vec2[i][offset] */
0328 ulong _nmod_vec_dot_pow2_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod);
0329 ulong _nmod_vec_dot1_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod);
0330 ulong _nmod_vec_dot2_half_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod);
0331 ulong _nmod_vec_dot2_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod);
0332 ulong _nmod_vec_dot3_acc_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod);
0333 ulong _nmod_vec_dot3_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod);
0334 #if FLINT_BITS == 64
0335 ulong _nmod_vec_dot2_split_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod, ulong pow2_precomp);
0336 #endif  // FLINT_BITS == 64
0337 
0338 
0339 
0340 
0341 
0342 /* ------------- dot product, general --------------- */
0343 
0344 // auxiliary for short dot products
0345 // (fixedlen small constant: compiler unrolls the loops completely)
0346 #define _NMOD_VEC_DOT_SHORT1(fixedlen,expr1,expr2)       \
0347         {                                                \
0348             ulong res = (expr1) * (expr2); i++;          \
0349             for (slong j = 0; j < fixedlen-1; j++, i++)  \
0350                 res += (expr1) * (expr2);                \
0351             NMOD_RED(res, res, mod);                     \
0352             return res;                                  \
0353         }                                                \
0354 
0355 #define _NMOD_VEC_DOT_SHORT2(fixedlen,expr1,expr2)          \
0356         {                                                   \
0357             ulong s0, s1, u0, u1;                           \
0358             umul_ppmm(u1, u0, (expr1), (expr2)); i++;       \
0359             for (slong j = 0; j < fixedlen-1; j++, i++)     \
0360             {                                               \
0361                 umul_ppmm(s1, s0, (expr1), (expr2));        \
0362                 add_ssaaaa(u1, u0, u1, u0, s1, s0);         \
0363             }                                               \
0364             NMOD2_RED2(s0, u1, u0, mod);                    \
0365             return s0;                                      \
0366         }                                                   \
0367 
0368 #define _NMOD_VEC_DOT_SHORT3(fixedlen,expr1,expr2)           \
0369         {                                                    \
0370             ulong t2 = UWORD(0);                             \
0371             ulong t1, t0;                                    \
0372             umul_ppmm(t1, t0, (expr1), (expr2)); i++;        \
0373             for (slong j = 0; j < fixedlen - 1; j++, i++)    \
0374             {                                                \
0375                 ulong s0, s1;                                \
0376                 umul_ppmm(s1, s0, (expr1), (expr2));         \
0377                 add_sssaaaaaa(t2, t1, t0,                    \
0378                                 t2, t1, t0,                  \
0379                                 UWORD(0), s1, s0);           \
0380             }                                                \
0381                                                              \
0382             NMOD_RED(t2, t2, mod);                           \
0383             ulong res;                                       \
0384             NMOD_RED3(res, t2, t1, t0, mod);                 \
0385             return res;                                      \
0386         }                                                    \
0387 
0388 // * supports 1 <= len <= 11, requires method==DOT1|DOT2|DOT3|DOT_POW2
0389 // * i must be already initialized at the first wanted value
0390 #define _NMOD_VEC_DOT_SHORT(i, expr1, expr2, len, mod, method)          \
0391 {                                                                       \
0392     if (method == _DOT1 || method == _DOT_POW2)                         \
0393     {                                                                   \
0394         if (len ==  1) _NMOD_VEC_DOT_SHORT1( 1, expr1, expr2)           \
0395         if (len ==  2) _NMOD_VEC_DOT_SHORT1( 2, expr1, expr2)           \
0396         if (len ==  3) _NMOD_VEC_DOT_SHORT1( 3, expr1, expr2)           \
0397         if (len ==  4) _NMOD_VEC_DOT_SHORT1( 4, expr1, expr2)           \
0398         if (len ==  5) _NMOD_VEC_DOT_SHORT1( 5, expr1, expr2)           \
0399         if (len ==  6) _NMOD_VEC_DOT_SHORT1( 6, expr1, expr2)           \
0400         if (len ==  7) _NMOD_VEC_DOT_SHORT1( 7, expr1, expr2)           \
0401         if (len ==  8) _NMOD_VEC_DOT_SHORT1( 8, expr1, expr2)           \
0402         if (len ==  9) _NMOD_VEC_DOT_SHORT1( 9, expr1, expr2)           \
0403         if (len == 10) _NMOD_VEC_DOT_SHORT1(10, expr1, expr2)           \
0404         _NMOD_VEC_DOT_SHORT1(11, expr1, expr2)                          \
0405     }                                                                   \
0406                                                                         \
0407     else if (method == _DOT2)                                           \
0408     {                                                                   \
0409         if (len ==  1) return nmod_mul((expr1), (expr2), mod);          \
0410         if (len ==  2) _NMOD_VEC_DOT_SHORT2( 2, expr1, expr2)           \
0411         if (len ==  3) _NMOD_VEC_DOT_SHORT2( 3, expr1, expr2)           \
0412         if (len ==  4) _NMOD_VEC_DOT_SHORT2( 4, expr1, expr2)           \
0413         if (len ==  5) _NMOD_VEC_DOT_SHORT2( 5, expr1, expr2)           \
0414         if (len ==  6) _NMOD_VEC_DOT_SHORT2( 6, expr1, expr2)           \
0415         if (len ==  7) _NMOD_VEC_DOT_SHORT2( 7, expr1, expr2)           \
0416         if (len ==  8) _NMOD_VEC_DOT_SHORT2( 8, expr1, expr2)           \
0417         if (len ==  9) _NMOD_VEC_DOT_SHORT2( 9, expr1, expr2)           \
0418         if (len == 10) _NMOD_VEC_DOT_SHORT2(10, expr1, expr2)           \
0419         _NMOD_VEC_DOT_SHORT2(11, expr1, expr2)                          \
0420     }                                                                   \
0421                                                                         \
0422     else if (method == _DOT3)                                           \
0423     {                                                                   \
0424         if (len ==  1) return nmod_mul((expr1), (expr2), mod);          \
0425         if (len ==  2) _NMOD_VEC_DOT_SHORT3( 2, expr1, expr2)           \
0426         if (len ==  3) _NMOD_VEC_DOT_SHORT3( 3, expr1, expr2)           \
0427         if (len ==  4) _NMOD_VEC_DOT_SHORT3( 4, expr1, expr2)           \
0428         if (len ==  5) _NMOD_VEC_DOT_SHORT3( 5, expr1, expr2)           \
0429         if (len ==  6) _NMOD_VEC_DOT_SHORT3( 6, expr1, expr2)           \
0430         if (len ==  7) _NMOD_VEC_DOT_SHORT3( 7, expr1, expr2)           \
0431         if (len ==  8) _NMOD_VEC_DOT_SHORT3( 8, expr1, expr2)           \
0432         if (len ==  9) _NMOD_VEC_DOT_SHORT3( 9, expr1, expr2)           \
0433         if (len == 10) _NMOD_VEC_DOT_SHORT3(10, expr1, expr2)           \
0434         _NMOD_VEC_DOT_SHORT3(11, expr1, expr2)                          \
0435     }                                                                   \
0436 }  while(0);                                                            \
0437 
0438 FLINT_FORCE_INLINE ulong _nmod_vec_dot(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod, dot_params_t params)
0439 {
0440     if (len <= 11)
0441     {
0442         if (len == 0) return UWORD(0);
0443         slong i = 0;
0444         _NMOD_VEC_DOT_SHORT(i, vec1[i], vec2[i], len, mod, params.method);
0445     }
0446 
0447     if (params.method == _DOT1)
0448         return _nmod_vec_dot1(vec1, vec2, len, mod);
0449 
0450 #if FLINT_BITS == 64
0451     if (params.method == _DOT2_SPLIT)
0452         return _nmod_vec_dot2_split(vec1, vec2, len, mod, params.pow2_precomp);
0453 #endif // FLINT_BITS == 64
0454 
0455     if (params.method == _DOT2)
0456         return _nmod_vec_dot2(vec1, vec2, len, mod);
0457 
0458     if (params.method == _DOT3_ACC)
0459         return _nmod_vec_dot3_acc(vec1, vec2, len, mod);
0460 
0461     if (params.method == _DOT3)
0462         return _nmod_vec_dot3(vec1, vec2, len, mod);
0463 
0464     if (params.method == _DOT2_HALF)
0465         return _nmod_vec_dot2_half(vec1, vec2, len, mod, params.pow2_precomp);
0466 
0467     if (params.method == _DOT_POW2)
0468     {
0469         if (mod.n <= UWORD(1) << (FLINT_BITS / 2))
0470             return _nmod_vec_dot1(vec1, vec2, len, mod);
0471         else
0472             return _nmod_vec_dot_pow2(vec1, vec2, len, mod);
0473     }
0474 
0475     // covers _DOT0 for len > 11 (i.e. mod.n == 1...)
0476     return UWORD(0);
0477 }
0478 
0479 FLINT_FORCE_INLINE ulong _nmod_vec_dot_rev(nn_srcptr vec1, nn_srcptr vec2, slong len, nmod_t mod, dot_params_t params)
0480 {
0481     if (len <= 11)
0482     {
0483         if (len == 0) return UWORD(0);
0484         slong i = 0;
0485         _NMOD_VEC_DOT_SHORT(i, vec1[i], vec2[len-1-i], len, mod, params.method);
0486     }
0487 
0488     if (params.method == _DOT1)
0489         return _nmod_vec_dot1_rev(vec1, vec2, len, mod);
0490 
0491 #if FLINT_BITS == 64
0492     if (params.method == _DOT2_SPLIT)
0493         return _nmod_vec_dot2_split_rev(vec1, vec2, len, mod, params.pow2_precomp);
0494 #endif // FLINT_BITS == 64
0495 
0496     if (params.method == _DOT2)
0497         return _nmod_vec_dot2_rev(vec1, vec2, len, mod);
0498 
0499     if (params.method == _DOT3_ACC)
0500         return _nmod_vec_dot3_acc_rev(vec1, vec2, len, mod);
0501 
0502     if (params.method == _DOT3)
0503         return _nmod_vec_dot3_rev(vec1, vec2, len, mod);
0504 
0505     if (params.method == _DOT2_HALF)
0506         return _nmod_vec_dot2_half_rev(vec1, vec2, len, mod);
0507 
0508     if (params.method == _DOT_POW2)
0509     {
0510         if (mod.n <= UWORD(1) << (FLINT_BITS / 2))
0511             return _nmod_vec_dot1_rev(vec1, vec2, len, mod);
0512         else
0513             return _nmod_vec_dot_pow2_rev(vec1, vec2, len, mod);
0514     }
0515 
0516     // covers _DOT0 for len > 11 (i.e. mod.n == 1...)
0517     return UWORD(0);
0518 }
0519 
0520 FLINT_FORCE_INLINE ulong _nmod_vec_dot_ptr(nn_srcptr vec1, const nn_ptr * vec2, slong offset, slong len, nmod_t mod, dot_params_t params)
0521 {
0522     if (len <= 11)
0523     {
0524         if (len == 0) return UWORD(0);
0525         slong i = 0;
0526         _NMOD_VEC_DOT_SHORT(i, vec1[i], vec2[i][offset], len, mod, params.method);
0527     }
0528 
0529     if (params.method == _DOT1)
0530         return _nmod_vec_dot1_ptr(vec1, vec2, offset, len, mod);
0531 
0532 #if FLINT_BITS == 64
0533     if (params.method == _DOT2_SPLIT)
0534         return _nmod_vec_dot2_split_ptr(vec1, vec2, offset, len, mod, params.pow2_precomp);
0535 #endif // FLINT_BITS == 64
0536 
0537     if (params.method == _DOT2)
0538         return _nmod_vec_dot2_ptr(vec1, vec2, offset, len, mod);
0539 
0540     if (params.method == _DOT3_ACC)
0541         return _nmod_vec_dot3_acc_ptr(vec1, vec2, offset, len, mod);
0542 
0543     if (params.method == _DOT3)
0544         return _nmod_vec_dot3_ptr(vec1, vec2, offset, len, mod);
0545 
0546     if (params.method == _DOT2_HALF)
0547         return _nmod_vec_dot2_half_ptr(vec1, vec2, offset, len, mod);
0548 
0549     if (params.method == _DOT_POW2)
0550     {
0551         if (mod.n <= UWORD(1) << (FLINT_BITS / 2))
0552             return _nmod_vec_dot1_ptr(vec1, vec2, offset, len, mod);
0553         else
0554             return _nmod_vec_dot_pow2_ptr(vec1, vec2, offset, len, mod);
0555     }
0556 
0557     // covers _DOT0 for len > 11 (i.e. mod.n == 1...)
0558     return UWORD(0);
0559 }
0560 
0561 #undef _NMOD_VEC_DOT_SHORT1
0562 #undef _NMOD_VEC_DOT_SHORT2
0563 #undef _NMOD_VEC_DOT_SHORT3
0564 
0565 
0566 /* ---- macros for dot product with expressions, specific algorithms ---- */
0567 
0568 // _DOT1   (1 limb)
0569 #define _NMOD_VEC_DOT1(res, i, len, expr1, expr2, mod)                \
0570 do                                                                    \
0571 {                                                                     \
0572     res = UWORD(0);                                                   \
0573     for (i = 0; i < (len); i++)                                       \
0574         res += (expr1) * (expr2);                                     \
0575     NMOD_RED(res, res, mod);                                          \
0576 } while(0);
0577 
0578 // _DOT2_SPLIT   (2 limbs, splitting at DOT_SPLIT_BITS bits, 8-unrolling)
0579 #if (FLINT_BITS == 64)
0580 #define _NMOD_VEC_DOT2_SPLIT(res, i, len, expr1, expr2, mod, pow2_precomp) \
0581 do                                                    \
0582 {                                                     \
0583     ulong dp_lo = 0;                                  \
0584     ulong dp_hi = 0;                                  \
0585                                                       \
0586     for (i = 0; i+7 < (len); )                        \
0587     {                                                 \
0588         dp_lo += (expr1) * (expr2); i++;              \
0589         dp_lo += (expr1) * (expr2); i++;              \
0590         dp_lo += (expr1) * (expr2); i++;              \
0591         dp_lo += (expr1) * (expr2); i++;              \
0592         dp_lo += (expr1) * (expr2); i++;              \
0593         dp_lo += (expr1) * (expr2); i++;              \
0594         dp_lo += (expr1) * (expr2); i++;              \
0595         dp_lo += (expr1) * (expr2); i++;              \
0596                                                       \
0597         dp_hi += dp_lo >> DOT_SPLIT_BITS;             \
0598         dp_lo &= DOT_SPLIT_MASK;                      \
0599     }                                                 \
0600                                                       \
0601     for ( ; i < (len); i++)                           \
0602         dp_lo += (expr1) * (expr2);                   \
0603                                                       \
0604     res = pow2_precomp * dp_hi + dp_lo;               \
0605     NMOD_RED(res, res, mod);                          \
0606 } while(0);
0607 #endif  // FLINT_BITS == 64
0608 
0609 // _DOT2_HALF   (two limbs, modulus < 2**32)
0610 // mod.n is too close to 2**32 to accumulate in some ulong
0611 // still interesting: a bit faster than _NMOD_VEC_DOT2
0612 #define _NMOD_VEC_DOT2_HALF(res, i, len, expr1, expr2, mod)           \
0613 do                                                                    \
0614 {                                                                     \
0615     ulong s0zz = UWORD(0);                                            \
0616     ulong s1zz = UWORD(0);                                            \
0617     for (i = 0; i < (len); i++)                                       \
0618     {                                                                 \
0619         const ulong prodzz = (expr1) * (expr2);                       \
0620         add_ssaaaa(s1zz, s0zz, s1zz, s0zz, 0, prodzz);                \
0621     }                                                                 \
0622     NMOD2_RED2(res, s1zz, s0zz, mod);                                 \
0623 } while(0);
0624 
0625 // _DOT2   (two limbs, general)
0626 #define _NMOD_VEC_DOT2(res, i, len, expr1, expr2, mod)                \
0627 do                                                                    \
0628 {                                                                     \
0629     ulong u0zz = UWORD(0);                                            \
0630     ulong u1zz = UWORD(0);                                            \
0631                                                                       \
0632     for (i = 0; i+7 < (len); )                                        \
0633     {                                                                 \
0634         ulong s0zz, s1zz;                                             \
0635         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0636         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0637         i++;                                                          \
0638         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0639         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0640         i++;                                                          \
0641         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0642         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0643         i++;                                                          \
0644         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0645         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0646         i++;                                                          \
0647         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0648         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0649         i++;                                                          \
0650         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0651         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0652         i++;                                                          \
0653         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0654         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0655         i++;                                                          \
0656         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0657         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0658         i++;                                                          \
0659     }                                                                 \
0660     for ( ; i < (len); i++)                                           \
0661     {                                                                 \
0662         ulong s0zz, s1zz;                                             \
0663         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0664         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0665     }                                                                 \
0666                                                                       \
0667     NMOD2_RED2(res, u1zz, u0zz, mod);                                 \
0668 } while(0);
0669 
0670 // _DOT3_ACC   (three limbs, delayed accumulations)
0671 // 8-unroll: requires  8 * (mod.n - 1)**2 < 2**128
0672 #define _NMOD_VEC_DOT3_ACC(res, i, len, expr1, expr2, mod)            \
0673 do                                                                    \
0674 {                                                                     \
0675     ulong t2zz = UWORD(0);                                            \
0676     ulong t1zz = UWORD(0);                                            \
0677     ulong t0zz = UWORD(0);                                            \
0678                                                                       \
0679     for (i = 0; i+7 < (len); )                                        \
0680     {                                                                 \
0681         ulong s0zz, s1zz;                                             \
0682         ulong u0zz = UWORD(0);                                        \
0683         ulong u1zz = UWORD(0);                                        \
0684         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0685         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0686         i++;                                                          \
0687         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0688         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0689         i++;                                                          \
0690         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0691         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0692         i++;                                                          \
0693         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0694         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0695         i++;                                                          \
0696         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0697         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0698         i++;                                                          \
0699         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0700         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0701         i++;                                                          \
0702         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0703         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0704         i++;                                                          \
0705         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0706         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0707         i++;                                                          \
0708         add_sssaaaaaa(t2zz, t1zz, t0zz,                               \
0709                         t2zz, t1zz, t0zz,                             \
0710                         UWORD(0), u1zz, u0zz);                        \
0711     }                                                                 \
0712                                                                       \
0713     ulong s0zz, s1zz;                                                 \
0714     ulong u0zz = UWORD(0);                                            \
0715     ulong u1zz = UWORD(0);                                            \
0716     for ( ; i < (len); i++)                                           \
0717     {                                                                 \
0718         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0719         add_ssaaaa(u1zz, u0zz, u1zz, u0zz, s1zz, s0zz);               \
0720     }                                                                 \
0721                                                                       \
0722     add_sssaaaaaa(t2zz, t1zz, t0zz,                                   \
0723                     t2zz, t1zz, t0zz,                                 \
0724                     UWORD(0), u1zz, u0zz);                            \
0725                                                                       \
0726     NMOD_RED(t2zz, t2zz, mod);                                        \
0727     NMOD_RED3(res, t2zz, t1zz, t0zz, mod);                            \
0728 } while(0);
0729 
0730 // _DOT3   (three limbs, general)
0731 // mod.n is too close to 2**64 to accumulate in two words
0732 #define _NMOD_VEC_DOT3(res, i, len, expr1, expr2, mod)                \
0733 do                                                                    \
0734 {                                                                     \
0735     ulong t2zz = UWORD(0);                                            \
0736     ulong t1zz = UWORD(0);                                            \
0737     ulong t0zz = UWORD(0);                                            \
0738     for (i = 0; i < (len); i++)                                       \
0739     {                                                                 \
0740         ulong s0zz, s1zz;                                             \
0741         umul_ppmm(s1zz, s0zz, (expr1), (expr2));                      \
0742         add_sssaaaaaa(t2zz, t1zz, t0zz,                               \
0743                         t2zz, t1zz, t0zz,                             \
0744                         UWORD(0), s1zz, s0zz);                        \
0745     }                                                                 \
0746                                                                       \
0747     NMOD_RED(t2zz, t2zz, mod);                                        \
0748     NMOD_RED3(res, t2zz, t1zz, t0zz, mod);                            \
0749 } while(0);
0750 
0751 /* ---- macros for dot product with expressions, general ---- */
0752 // currently no vectorization here
0753 
0754 #if (FLINT_BITS == 64)
0755 
0756 #define NMOD_VEC_DOT(res, i, len, expr1, expr2, mod, params)          \
0757 do                                                                    \
0758 {                                                                     \
0759     res = UWORD(0);   /* covers _DOT0 */                              \
0760     if (params.method == _DOT1 || params.method == _DOT_POW2)         \
0761         _NMOD_VEC_DOT1(res, i, len, expr1, expr2, mod)                \
0762     else if (params.method == _DOT2_SPLIT)                            \
0763         _NMOD_VEC_DOT2_SPLIT(res, i, len, expr1, expr2, mod,          \
0764                 params.pow2_precomp)                                  \
0765     else if (params.method == _DOT2_HALF)                             \
0766         _NMOD_VEC_DOT2_HALF(res, i, len, expr1, expr2, mod)           \
0767     else if (params.method == _DOT2)                                  \
0768         _NMOD_VEC_DOT2(res, i, len, expr1, expr2, mod)                \
0769     else if (params.method == _DOT3_ACC)                              \
0770         _NMOD_VEC_DOT3_ACC(res, i, len, expr1, expr2, mod)            \
0771     else if (params.method == _DOT3)                                  \
0772         _NMOD_VEC_DOT3(res, i, len, expr1, expr2, mod)                \
0773 } while(0);
0774 
0775 #else  // FLINT_BITS == 64
0776 
0777 #define NMOD_VEC_DOT(res, i, len, expr1, expr2, mod, params)          \
0778 do                                                                    \
0779 {                                                                     \
0780     res = UWORD(0);   /* covers _DOT0 */                              \
0781     if (params.method == _DOT1 || params.method == _DOT_POW2)         \
0782         _NMOD_VEC_DOT1(res, i, len, expr1, expr2, mod)                \
0783     else if (params.method == _DOT2_HALF)                             \
0784         _NMOD_VEC_DOT2_HALF(res, i, len, expr1, expr2, mod)           \
0785     else if (params.method == _DOT2)                                  \
0786         _NMOD_VEC_DOT2(res, i, len, expr1, expr2, mod)                \
0787     else if (params.method == _DOT3_ACC)                              \
0788         _NMOD_VEC_DOT3_ACC(res, i, len, expr1, expr2, mod)            \
0789     else if (params.method == _DOT3)                                  \
0790         _NMOD_VEC_DOT3(res, i, len, expr1, expr2, mod)                \
0791 } while(0);
0792 
0793 #endif  // FLINT_BITS == 64
0794 
0795 /* Note: if de-inlining this, consider specializing for the common use
0796    case of stride1 == 1 or stride2 == 1. */
0797 NMOD_VEC_INLINE ulong _nmod_vec_dot_strided(nn_srcptr vec1, slong stride1,
0798     nn_srcptr vec2, slong stride2, slong len, nmod_t mod, dot_params_t params)
0799 {
0800     slong i;
0801     ulong r;
0802     NMOD_VEC_DOT(r, i, len, vec1[i * stride1], vec2[i * stride2], mod, params)
0803     return r;
0804 }
0805 
0806 #ifdef __cplusplus
0807 }
0808 #endif
0809 
0810 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!