include/flint/crt_helpers.h

0001 /*
0002     Copyright (C) 2022 Daniel Schultz
0003
0004     This file is part of FLINT.
0005
0006     FLINT is free software: you can redistribute it and/or modify it under
0007     the terms of the GNU Lesser General Public License (LGPL) as published
0008     by the Free Software Foundation; either version 3 of the License, or
0009     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0010 */
0011
0012 #ifndef CRT_HELPERS_H
0013 #define CRT_HELPERS_H
0014
0015 #if defined(__GNUC__) && defined(__AVX2__)
0016 # include <immintrin.h>
0017 #elif defined(_MSC_VER) && defined(__AVX2__)
0018 # include <intrin.h>
0019 #endif
0020
0021 #include "longlong.h"
0022 #include "templates.h"
0023
0024 #ifdef __cplusplus
0025 extern "C" {
0026 #endif
0027
0028 #if defined(__AVX2__)
0029
0030 FLINT_FORCE_INLINE unsigned char _addcarry_ulong(unsigned char cf, ulong x, ulong y, ulong* s)
0031 {
0032     long long unsigned int _s;
0033     cf = _addcarry_u64(cf, (long long unsigned int)(x),
0034                            (long long unsigned int)(y),
0035                            &_s);
0036     *s = (ulong)(_s);
0037     return cf;
0038 }
0039
0040 FLINT_FORCE_INLINE unsigned char _subborrow_ulong(unsigned char cf, ulong x, ulong y, ulong* s)
0041 {
0042     long long unsigned int _s;
0043     cf = _subborrow_u64(cf, (long long unsigned int)(x),
0044                             (long long unsigned int)(y),
0045                            &_s);
0046     *s = (ulong)(_s);
0047     return cf;
0048 }
0049
0050 #else
0051
0052 FLINT_FORCE_INLINE unsigned char _addcarry_ulong(unsigned char cf, ulong x, ulong y, ulong* s)
0053 {
0054 #if 0
0055     ulong cf2;
0056     *s = __builtin_addcl(x, y, cf, &cf2);
0057     return cf2;
0058 #else
0059     ulong hi, lo;
0060     add_ssaaaa(hi, lo, 0, x, 0, y);
0061     add_ssaaaa(hi, lo, hi, lo, 0, (ulong) cf);
0062     *s = lo;
0063     return hi;
0064 #endif
0065 }
0066
0067 FLINT_FORCE_INLINE unsigned char _subborrow_ulong(unsigned char cf, ulong x, ulong y, ulong* s)
0068 {
0069 #if 0
0070     ulong cf2;
0071     *s = __builtin_subcl(x, y, cf, &cf2);
0072     return cf2;
0073 #else
0074     ulong hi, lo;
0075     sub_ddmmss(hi, lo, 0, x, 0, y);
0076     sub_ddmmss(hi, lo, hi, lo, 0, (ulong) cf);
0077     *s = lo;
0078     return hi != 0;
0079 #endif
0080 }
0081
0082 #endif
0083
0084
0085 #if 1
0086
0087 #if defined(__GNUC__) && defined(__AVX2__)
0088
0089
0090 #elif defined(__GNUC__) && defined(__ARM_NEON)
0091
0092
0093
0094 #elif defined(_MSC_VER) && (defined(__AVX2__) || defined(_M_ARM64))
0095
0096
0097
0098 #else
0099 # error crt_helpers.h requires AVX2 or Neon instructions
0100 #endif
0101
0102 FLINT_FORCE_INLINE void multi_add_0(ulong FLINT_UNUSED(z[]), const ulong FLINT_UNUSED(a[]))
0103 {
0104 }
0105
0106 FLINT_FORCE_INLINE void multi_add_1(ulong z[], const ulong a[])
0107 {
0108     z[0] += a[0];
0109 }
0110
0111 FLINT_FORCE_INLINE void multi_add_2(ulong z[], const ulong a[])
0112 {
0113     add_ssaaaa(z[1],z[0],
0114                z[1],z[0],
0115                a[1],a[0]);
0116 }
0117
0118 FLINT_FORCE_INLINE void multi_add_3(ulong z[], const ulong a[])
0119 {
0120     add_sssaaaaaa(z[2],z[1],z[0],
0121                   z[2],z[1],z[0],
0122                   a[2],a[1],a[0]);
0123 }
0124
0125 FLINT_FORCE_INLINE void multi_add_4(ulong z[], const ulong a[])
0126 {
0127     add_ssssaaaaaaaa(z[3],z[2],z[1],z[0],
0128                      z[3],z[2],z[1],z[0],
0129                      a[3],a[2],a[1],a[0]);
0130 }
0131
0132 FLINT_FORCE_INLINE void multi_add_5(ulong z[], const ulong a[])
0133 {
0134     add_sssssaaaaaaaaaa(z[4],z[3],z[2],z[1],z[0],
0135                         z[4],z[3],z[2],z[1],z[0],
0136                         a[4],a[3],a[2],a[1],a[0]);
0137 }
0138
0139 FLINT_FORCE_INLINE void multi_add_6(ulong z[], const ulong a[])
0140 {
0141     add_ssssssaaaaaaaaaaaa(z[5],z[4],z[3],z[2],z[1],z[0],
0142                            z[5],z[4],z[3],z[2],z[1],z[0],
0143                            a[5],a[4],a[3],a[2],a[1],a[0]);
0144 }
0145
0146 FLINT_FORCE_INLINE void multi_add_7(ulong z[], const ulong a[])
0147 {
0148     add_sssssssaaaaaaaaaaaaaa(z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0149                               z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0150                               a[6],a[5],a[4],a[3],a[2],a[1],a[0]);
0151 }
0152
0153 FLINT_FORCE_INLINE void multi_add_8(ulong z[], const ulong a[])
0154 {
0155     add_ssssssssaaaaaaaaaaaaaaaa(z[7],z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0156                                  z[7],z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0157                                  a[7],a[6],a[5],a[4],a[3],a[2],a[1],a[0]);
0158 }
0159
0160 FLINT_FORCE_INLINE void multi_sub_0(ulong FLINT_UNUSED(z[]), const ulong FLINT_UNUSED(a[]))
0161 {
0162 }
0163
0164 FLINT_FORCE_INLINE void multi_sub_1(ulong z[], const ulong a[])
0165 {
0166     z[0] -= a[0];
0167 }
0168
0169 FLINT_FORCE_INLINE void multi_sub_2(ulong z[], const ulong a[])
0170 {
0171     sub_ddmmss(z[1],z[0],
0172                z[1],z[0],
0173                a[1],a[0]);
0174 }
0175
0176 FLINT_FORCE_INLINE void multi_sub_3(ulong z[], const ulong a[])
0177 {
0178     sub_dddmmmsss(z[2],z[1],z[0],
0179                   z[2],z[1],z[0],
0180                   a[2],a[1],a[0]);
0181 }
0182
0183 FLINT_FORCE_INLINE void multi_sub_4(ulong z[], const ulong a[])
0184 {
0185     sub_ddddmmmmssss(z[3],z[2],z[1],z[0],
0186                      z[3],z[2],z[1],z[0],
0187                      a[3],a[2],a[1],a[0]);
0188 }
0189
0190 FLINT_FORCE_INLINE void multi_sub_5(ulong z[], const ulong a[])
0191 {
0192     sub_dddddmmmmmsssss(z[4],z[3],z[2],z[1],z[0],
0193                         z[4],z[3],z[2],z[1],z[0],
0194                         a[4],a[3],a[2],a[1],a[0]);
0195 }
0196
0197 FLINT_FORCE_INLINE void multi_sub_6(ulong z[], const ulong a[])
0198 {
0199     sub_ddddddmmmmmmssssss(z[5],z[4],z[3],z[2],z[1],z[0],
0200                            z[5],z[4],z[3],z[2],z[1],z[0],
0201                            a[5],a[4],a[3],a[2],a[1],a[0]);
0202 }
0203
0204 FLINT_FORCE_INLINE void multi_sub_7(ulong z[], const ulong a[])
0205 {
0206     sub_dddddddmmmmmmmsssssss(z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0207                               z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0208                               a[6],a[5],a[4],a[3],a[2],a[1],a[0]);
0209 }
0210
0211 FLINT_FORCE_INLINE void multi_sub_8(ulong z[], const ulong a[])
0212 {
0213     sub_ddddddddmmmmmmmmssssssss(z[7],z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0214                                  z[7],z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0215                                  a[7],a[6],a[5],a[4],a[3],a[2],a[1],a[0]);
0216 }
0217
0218 FLINT_FORCE_INLINE void multi_rsub_0(ulong FLINT_UNUSED(z[]), const ulong FLINT_UNUSED(a[]))
0219 {
0220 }
0221
0222 FLINT_FORCE_INLINE void multi_rsub_1(ulong z[], const ulong a[])
0223 {
0224     z[0] = a[0] - z[0];
0225 }
0226
0227 FLINT_FORCE_INLINE void multi_rsub_2(ulong z[], const ulong a[])
0228 {
0229     sub_ddmmss(z[1],z[0],
0230                a[1],a[0],
0231                z[1],z[0]);
0232 }
0233
0234 FLINT_FORCE_INLINE void multi_rsub_3(ulong z[], const ulong a[])
0235 {
0236     sub_dddmmmsss(z[2],z[1],z[0],
0237                   a[2],a[1],a[0],
0238                   z[2],z[1],z[0]);
0239 }
0240
0241 FLINT_FORCE_INLINE void multi_rsub_4(ulong z[], const ulong a[])
0242 {
0243     sub_ddddmmmmssss(z[3],z[2],z[1],z[0],
0244                      a[3],a[2],a[1],a[0],
0245                      z[3],z[2],z[1],z[0]);
0246 }
0247
0248 FLINT_FORCE_INLINE void multi_rsub_5(ulong z[], const ulong a[])
0249 {
0250     sub_dddddmmmmmsssss(z[4],z[3],z[2],z[1],z[0],
0251                         a[4],a[3],a[2],a[1],a[0],
0252                         z[4],z[3],z[2],z[1],z[0]);
0253 }
0254
0255 FLINT_FORCE_INLINE void multi_rsub_6(ulong z[], const ulong a[])
0256 {
0257     sub_ddddddmmmmmmssssss(z[5],z[4],z[3],z[2],z[1],z[0],
0258                            a[5],a[4],a[3],a[2],a[1],a[0],
0259                            z[5],z[4],z[3],z[2],z[1],z[0]);
0260 }
0261
0262 FLINT_FORCE_INLINE void multi_rsub_7(ulong z[], const ulong a[])
0263 {
0264     sub_dddddddmmmmmmmsssssss(z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0265                               a[6],a[5],a[4],a[3],a[2],a[1],a[0],
0266                               z[6],z[5],z[4],z[3],z[2],z[1],z[0]);
0267 }
0268
0269 FLINT_FORCE_INLINE void multi_rsub_8(ulong z[], const ulong a[])
0270 {
0271     sub_ddddddddmmmmmmmmssssssss(z[7],z[6],z[5],z[4],z[3],z[2],z[1],z[0],
0272                                  a[7],a[6],a[5],a[4],a[3],a[2],a[1],a[0],
0273                                  z[7],z[6],z[5],z[4],z[3],z[2],z[1],z[0]);
0274 }
0275
0276 #else
0277
0278
0279 #define DEFINE_IT(n) \
0280 FLINT_FORCE_INLINE void CAT(multi_add, n)(ulong z[], const ulong a[]) \
0281 { \
0282     unsigned char cf = 0; \
0283     for (ulong i = 0; i < n; i++) \
0284         cf = _addcarry_ulong(cf, z[i], a[i], &z[i]); \
0285 }
0286
0287 DEFINE_IT(0)
0288 DEFINE_IT(1)
0289 DEFINE_IT(2)
0290 DEFINE_IT(3)
0291 DEFINE_IT(4)
0292 DEFINE_IT(5)
0293 DEFINE_IT(6)
0294 DEFINE_IT(7)
0295 DEFINE_IT(8)
0296 #undef DEFINE_IT
0297
0298 #define DEFINE_IT(n) \
0299 FLINT_FORCE_INLINE void CAT(multi_sub, n)(ulong z[], const ulong a[]) \
0300 { \
0301     unsigned char cf = 0; \
0302     for (ulong i = 0; i < n; i++) \
0303         cf = _subborrow_ulong(cf, z[i], a[i], &z[i]); \
0304 }
0305
0306 DEFINE_IT(0)
0307 DEFINE_IT(1)
0308 DEFINE_IT(2)
0309 DEFINE_IT(3)
0310 DEFINE_IT(4)
0311 DEFINE_IT(5)
0312 DEFINE_IT(6)
0313 DEFINE_IT(7)
0314 #undef DEFINE_IT
0315
0316 #endif
0317
0318 #ifdef __GNUC__
0319 FLINT_FORCE_INLINE void _mul(ulong* hi, ulong* lo, ulong y, ulong x)
0320 {
0321     __uint128_t p = ((__uint128_t) x) * ((__uint128_t) y);
0322     *lo = (ulong) (p);
0323     *hi = (ulong) (p >> 64);
0324 }
0325
0326 FLINT_FORCE_INLINE void _madd(ulong* hi, ulong* lo, ulong y, ulong x)
0327 {
0328     __uint128_t p = ((__uint128_t) *lo) | (((__uint128_t) *hi) << 64);
0329     p += ((__uint128_t) x) * ((__uint128_t) y);
0330     *lo = (ulong) (p);
0331     *hi = (ulong) (p >> 64);
0332 }
0333 #else
0334 FLINT_FORCE_INLINE void _mul(ulong* hi, ulong* lo, ulong y, ulong x)
0335 {
0336     ulong r1, r0;
0337     umul_ppmm(r1, r0, x, y);
0338     *lo = r0;
0339     *hi = r1;
0340 }
0341
0342 FLINT_FORCE_INLINE void _madd(ulong* hi, ulong* lo, ulong y, ulong x)
0343 {
0344     ulong r1, r0;
0345     umul_ppmm(r1, r0, x, y);
0346     add_ssaaaa(*hi, *lo, r1, r0, *hi, *lo);
0347 }
0348 #endif
0349
0350 /* NOTE: Define these manually to avoid compiler warnings or errors */
0351 FLINT_FORCE_INLINE void CAT3(_big_mul, 1, 0)(ulong r[], ulong FLINT_UNUSED(t[]), ulong FLINT_UNUSED(C[]), ulong FLINT_UNUSED(y))
0352 {
0353     r[0] = 0;
0354 }
0355
0356 FLINT_FORCE_INLINE void CAT3(_big_addmul, 1, 0)(ulong FLINT_UNUSED(r[]), ulong FLINT_UNUSED(t[]), ulong FLINT_UNUSED(C[]), ulong FLINT_UNUSED(y))
0357 {
0358 }
0359
0360 #define DEFINE_IT(n, m) \
0361 FLINT_FORCE_INLINE void CAT3(_big_mul, n, m)(ulong r[], ulong t[], ulong C[], ulong y) \
0362 { \
0363     for (ulong k = 0; k < n; k += 2) \
0364     { \
0365         if (k + 1 < n) \
0366         { \
0367             FLINT_ASSERT(k < m); \
0368             _mul(&r[k+1],&r[k+0], C[k+0], y); \
0369         } \
0370         else \
0371         { \
0372             FLINT_ASSERT(k + 1 == n); \
0373             if (k < m) \
0374                 r[k+0] = C[k+0]*y; \
0375             else \
0376                 r[k+0] = 0; \
0377         } \
0378  \
0379         if (k + 2 < n) \
0380         { \
0381             FLINT_ASSERT(k + 1 < m); \
0382             _mul(&t[k+2],&t[k+1], C[k+1], y); \
0383         } \
0384         else if (k + 1 < n) \
0385         { \
0386             if (k + 1 < m) \
0387                 t[k+1] = C[k+1]*y; \
0388             else \
0389                 t[k+1] = 0; \
0390         } \
0391     } \
0392 } \
0393 FLINT_FORCE_INLINE void CAT3(_big_addmul, n, m)(ulong r[], ulong t[], ulong C[], ulong y) \
0394 { \
0395     for (ulong k = 0; k < n; k += 2) \
0396     { \
0397         if (k + 1 < n) \
0398         { \
0399             FLINT_ASSERT(k < m); \
0400             _madd(&r[k+1],&r[k+0], C[k+0], y); \
0401         } \
0402         else \
0403         { \
0404             FLINT_ASSERT(k + 1 == n); \
0405             if (k < m) \
0406                 r[k+0] += C[k+0]*y; \
0407         } \
0408  \
0409         if (k + 2 < n) \
0410         { \
0411             FLINT_ASSERT(k + 1 < m); \
0412             _madd(&t[k+2],&t[k+1], C[k+1], y); \
0413         } \
0414         else if (k + 1 < n) \
0415         { \
0416             if (k + 1 < m) \
0417                 t[k+1] += C[k+1]*y; \
0418         } \
0419     } \
0420 }
0421
0422 DEFINE_IT(2, 1)
0423 DEFINE_IT(3, 2)
0424 DEFINE_IT(4, 3)
0425 DEFINE_IT(4, 4)
0426 DEFINE_IT(5, 4)
0427 DEFINE_IT(6, 5)
0428 DEFINE_IT(7, 6)
0429 #undef DEFINE_IT
0430
0431
0432
0433 #define DEFINE_IT(n, n_minus_1) \
0434 FLINT_FORCE_INLINE void CAT(_reduce_big_sum, n)(ulong r[], ulong t[], const ulong* limit) \
0435 { \
0436     CAT(multi_add, n_minus_1)(r+1, t+1); \
0437 check: \
0438     for (ulong k = n; k > 1; k--) \
0439     { \
0440         if (FLINT_LIKELY(r[k-1] > limit[k-1])) \
0441             goto sub; \
0442         if (r[k-1] < limit[k-1]) \
0443             return; \
0444     } \
0445     if (r[0] < limit[0]) \
0446         return; \
0447 sub: \
0448     CAT(multi_sub, n)(r, limit); \
0449     goto check; \
0450 }
0451
0452 DEFINE_IT(1, 0)
0453 DEFINE_IT(2, 1)
0454 DEFINE_IT(3, 2)
0455 DEFINE_IT(4, 3)
0456 DEFINE_IT(5, 4)
0457 DEFINE_IT(6, 5)
0458 DEFINE_IT(7, 6)
0459 #undef DEFINE_IT
0460
0461 #ifdef __cplusplus
0462 }
0463 #endif
0464
0465 #endif /* CRT_HELPERS_H */