include/flint/longlong.h

0001 /*
0002     Copyright 1991, 1992, 1993, 1994, 1996, 1997, 1999, 2000, 2001, 2002, 2003,
0003     2004, 2005 Free Software Foundation, Inc.
0004
0005     Copyright (C) 2009, 2015, 2016 William Hart
0006     Copyright (C) 2011 Fredrik Johansson
0007     Copyright (C) 2023 Albin Ahlbäck
0008
0009     This file is part of FLINT.
0010
0011     Contains code from GNU MP Library.
0012
0013     FLINT is free software: you can redistribute it and/or modify it under
0014     the terms of the GNU Lesser General Public License (LGPL) as published
0015     by the Free Software Foundation; either version 3 of the License, or
0016     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0017 */
0018
0019 #ifndef FLINT_LONGLONG_H
0020 #define FLINT_LONGLONG_H
0021
0022 #include "flint.h"
0023
0024 #ifdef __cplusplus
0025 extern "C" {
0026 #endif
0027
0028 /* FIXME: Bug in Clang 16+ */
0029 #if __clang_major__ >= 16
0030 # define CLANG_volatile volatile
0031 #else
0032 # define CLANG_volatile
0033 #endif
0034
0035 #if defined(__GNUC__)
0036
0037 /* Trailing and leading zeros */
0038 # if FLINT_LONG_LONG
0039 #  define flint_clz __builtin_clzll
0040 #  define flint_ctz __builtin_ctzll
0041 # else
0042 #  define flint_clz __builtin_clzl
0043 #  define flint_ctz __builtin_ctzl
0044 # endif
0045
0046 /* Addition, subtraction and multiplication */
0047 # if defined(__clang__)
0048 #  include "longlong_asm_clang.h"
0049 # elif FLINT_WANT_ASSEMBLY
0050 #  include "longlong_asm_gcc.h"
0051 # endif
0052
0053 # if FLINT_WANT_ASSEMBLY
0054 #  include "longlong_asm_gnu.h"
0055 #  include "longlong_div_gnu.h"
0056 # endif
0057
0058 #elif defined(_MSC_VER)
0059
0060 # if defined(_M_X64) || defined(_M_IX86)
0061 #  include "longlong_msc_x86.h"
0062 # elif defined(_M_ARM64)
0063 #  include "longlong_msc_arm64.h"
0064 # endif
0065
0066 #endif
0067
0068 /* Generics ******************************************************************/
0069
0070 #define __ll_B ((ulong) 1 << (FLINT_BITS / 2))
0071 #define __ll_lowpart(t) ((ulong) (t) & (__ll_B - 1))
0072 #define __ll_highpart(t) ((ulong) (t) >> (FLINT_BITS / 2))
0073
0074 /* Trailing and leading zeros */
0075 #if !defined(flint_ctz)
0076 # define NEED_CLZ_TAB 1
0077 FLINT_DLL extern const unsigned char __flint_clz_tab[128];
0078
0079 # define flint_clz flint_clz
0080 static inline int flint_clz(ulong x)
0081 {
0082     ulong a, xr = x;
0083     const unsigned int bits4 = FLINT_BITS / 4;
0084     if (FLINT_BITS == 32)
0085         a = xr < ((ulong) 1 << 2 * bits4)
0086         ? (xr < ((ulong) 1 << bits4) ? 1 : bits4 + 1)
0087         : (xr < ((ulong) 1 << 3 * bits4) ? 2 * bits4 + 1 : 3 * bits4 + 1);
0088     else
0089     {
0090         for (a = FLINT_BITS - 8; a > 0; a -= 8)
0091             if (((xr >> a) & 0xff) != 0)
0092                 break;
0093         ++a;
0094     }
0095     return FLINT_BITS + 1 - a - __flint_clz_tab[xr >> a];
0096 }
0097
0098 # define flint_ctz flint_ctz
0099 static inline int flint_ctz(ulong x)
0100 {
0101     return FLINT_BITS - 1 - flint_clz(x & -x);
0102 }
0103 #endif
0104
0105 /* Beware when using the unsigned return value in signed arithmetic */
0106 FLINT_FORCE_INLINE
0107 flint_bitcnt_t FLINT_BIT_COUNT(ulong x)
0108 {
0109     flint_bitcnt_t zeros = FLINT_BITS;
0110     if (x) zeros = flint_clz(x);
0111     return FLINT_BITS - zeros;
0112 }
0113
0114 #define FLINT_FLOG2(k) (FLINT_BIT_COUNT(k) - 1)
0115 #define FLINT_CLOG2(k) FLINT_BIT_COUNT((k) - 1)
0116
0117 /* Addition and subtraction */
0118 #if !defined(add_ssaaaa)
0119 # define add_ssaaaa(s1, s0, a1, a0, b1, b0) \
0120   do { \
0121     ulong __t0 = (a0); \
0122     (s0) = (a0) + (b0); \
0123     (s1) = (a1) + (b1) + ((ulong) (s0) < __t0); \
0124   } while (0)
0125
0126 # define sub_ddmmss(s1, s0, a1, a0, b1, b0) \
0127   do { \
0128     ulong __t0 = (a0); \
0129     (s0) = (a0) - (b0); \
0130     (s1) = (a1) - (b1) - ((ulong) (s0) > __t0); \
0131   } while (0)
0132 #endif
0133
0134 #if !defined(add_sssaaaaaa)
0135 # define add_sssaaaaaa(s2, s1, s0, a2, a1, a0, b2, b1, b0) \
0136   do { \
0137     ulong __t1, __t2; \
0138     add_ssaaaa(__t1, s0, (ulong) 0, a0, (ulong) 0, b0); \
0139     add_ssaaaa(__t2, s1, (ulong) 0, a1, (ulong) 0, b1); \
0140     add_ssaaaa(s2, s1, (a2) + (b2), s1, __t2, __t1); \
0141   } while (0)
0142
0143 # define sub_dddmmmsss(d2, d1, d0, m2, m1, m0, s2, s1, s0) \
0144   do { \
0145     ulong __t1, __t2; \
0146     sub_ddmmss(__t1, d0, (ulong) 0, m0, (ulong) 0, s0); \
0147     sub_ddmmss(__t2, d1, (ulong) 0, m1, (ulong) 0, s1); \
0148     sub_ddmmss(d2, d1, (m2) - (s2), d1, -__t2, -__t1); \
0149   } while (0)
0150 #endif
0151
0152 #if !defined(add_ssssaaaaaaaa)
0153 # define add_ssssaaaaaaaa(s3, s2, s1, s0, a3, a2, a1, a0, b3, b2, b1, b0) \
0154   do { \
0155     ulong __t3; \
0156     add_sssaaaaaa(__t3, s1, s0, (ulong) 0, a1, a0, (ulong) 0, b1, b0); \
0157     add_ssaaaa(s3, s2, a3, a2, b3, b2); \
0158     add_ssaaaa(s3, s2, s3, s2, (ulong) 0, __t3); \
0159   } while (0)
0160
0161 # define sub_ddddmmmmssss(s3, s2, s1, s0, a3, a2, a1, a0, b3, b2, b1, b0) \
0162   do { \
0163     ulong __t3, __t4; \
0164     sub_dddmmmsss(__t3, s1, s0, (ulong) 0, a1, a0, (ulong) 0, b1, b0); \
0165     sub_ddmmss(__t4, s2, (ulong) 0, a2, (ulong) 0, b2); \
0166     sub_ddmmss(s3, s2, (a3) - (b3), s2, -__t4, -__t3); \
0167   } while (0)
0168 #endif
0169
0170 #if !defined(add_sssssaaaaaaaaaa)
0171 # define add_sssssaaaaaaaaaa(s4, s3, s2, s1, s0, a4, a3, a2, a1, a0, b4, b3, b2, b1, b0) \
0172   do { \
0173     ulong __t4 = 0; \
0174     add_ssssaaaaaaaa(__t4, s2, s1, s0, (ulong) 0, a2, a1, a0, (ulong) 0, b2, b1, b0); \
0175     add_ssaaaa(s4, s3, a4, a3, b4, b3); \
0176     add_ssaaaa(s4, s3, s4, s3, (ulong) 0, __t4); \
0177   } while (0)
0178
0179 # define sub_dddddmmmmmsssss(s4, s3, s2, s1, s0, a4, a3, a2, a1, a0, b4, b3, b2, b1, b0) \
0180   do { \
0181     ulong __t5, __t6; \
0182     sub_ddddmmmmssss(__t5, s2, s1, s0, (ulong) 0, a2, a1, a0, (ulong) 0, b2, b1, b0); \
0183     sub_ddmmss(__t6, s3, (ulong) 0, a3, (ulong) 0, b3); \
0184     sub_ddmmss(s4, s3, (a4) - (b4), s3, -__t6, -__t5); \
0185   } while (0)
0186 #endif
0187
0188 #if !defined(add_ssssssaaaaaaaaaaaa)
0189 # define add_ssssssaaaaaaaaaaaa(s5, s4, s3, s2, s1, s0, a5, a4, a3, a2, a1, a0, b5, b4, b3, b2, b1, b0) \
0190   do { \
0191     ulong __t5 = 0; \
0192     add_sssssaaaaaaaaaa(__t5, s3, s2, s1, s0, (ulong) 0, a3, a2, a1, a0, (ulong) 0, b3, b2, b1, b0); \
0193     add_ssaaaa(s5, s4, a5, a4, b5, b4); \
0194     add_ssaaaa(s5, s4, s5, s4, (ulong) 0, __t5); \
0195   } while (0)
0196
0197 # define sub_ddddddmmmmmmssssss(s5, s4, s3, s2, s1, s0, a5, a4, a3, a2, a1, a0, b5, b4, b3, b2, b1, b0) \
0198   do { \
0199     ulong __t7, __t8; \
0200     sub_dddddmmmmmsssss(__t7, s3, s2, s1, s0, (ulong) 0, a3, a2, a1, a0, (ulong) 0, b3, b2, b1, b0);\
0201     sub_ddmmss(__t8, s4, (ulong) 0, a4, (ulong) 0, b4); \
0202     sub_ddmmss(s5, s4, (a5) - (b5), s4, -__t8, -__t7); \
0203   } while (0)
0204 #endif
0205
0206 #if !defined(add_sssssssaaaaaaaaaaaaaa)
0207 # define add_sssssssaaaaaaaaaaaaaa(s6, s5, s4, s3, s2, s1, s0, a6, a5, a4, a3, a2, a1, a0, b6, b5, b4, b3, b2, b1, b0) \
0208   do { \
0209     ulong __t6 = 0; \
0210     add_ssssssaaaaaaaaaaaa(__t6, s4, s3, s2, s1, s0, (ulong) 0, a4, a3, a2, a1, a0, (ulong) 0, b4, b3, b2, b1, b0); \
0211     add_ssaaaa(s6, s5, a6, a5, b6, b5); \
0212     add_ssaaaa(s6, s5, s6, s5, (ulong) 0, __t6); \
0213   } while (0)
0214
0215 # define sub_dddddddmmmmmmmsssssss(s6, s5, s4, s3, s2, s1, s0, a6, a5, a4, a3, a2, a1, a0, b6, b5, b4, b3, b2, b1, b0) \
0216   do { \
0217     ulong __t9, __t10; \
0218     sub_ddddddmmmmmmssssss(__t9, s4, s3, s2, s1, s0, (ulong) 0, a4, a3, a2, a1, a0, (ulong) 0, b4, b3, b2, b1, b0); \
0219     sub_ddmmss(__t10, s5, (ulong) 0, a5, (ulong) 0, b5); \
0220     sub_ddmmss(s6, s5, (a6) - (b6), s5, -__t10, -__t9); \
0221   } while (0)
0222 #endif
0223
0224 #if !defined(add_ssssssssaaaaaaaaaaaaaaaa)
0225 # define add_ssssssssaaaaaaaaaaaaaaaa(s7, s6, s5, s4, s3, s2, s1, s0, a7, a6, a5, a4, a3, a2, a1, a0, b7, b6, b5, b4, b3, b2, b1, b0) \
0226   do { \
0227     ulong __t7 = 0; \
0228     add_sssssssaaaaaaaaaaaaaa(__t7, s5, s4, s3, s2, s1, s0, (ulong) 0, a5, a4, a3, a2, a1, a0, (ulong) 0, b5, b4, b3, b2, b1, b0); \
0229     add_ssaaaa(s7, s6, a7, a6, b7, b6); \
0230     add_ssaaaa(s7, s6, s7, s6, (ulong) 0, __t7); \
0231   } while (0)
0232
0233 # define sub_ddddddddmmmmmmmmssssssss(s7, s6, s5, s4, s3, s2, s1, s0, a7, a6, a5, a4, a3, a2, a1, a0, b7, b6, b5, b4, b3, b2, b1, b0) \
0234   do { \
0235     ulong __t11, __t12; \
0236     sub_dddddddmmmmmmmsssssss(__t11, s5, s4, s3, s2, s1, s0, (ulong) 0, a5, a4, a3, a2, a1, a0, (ulong) 0, b5, b4, b3, b2, b1, b0); \
0237     sub_ddmmss(__t12, s6, (ulong) 0, a6, (ulong) 0, b6); \
0238     sub_ddmmss(s7, s6, (a7) - (b7), s6, -__t12, -__t11); \
0239   } while (0)
0240 #endif
0241
0242 #if !defined(MPN_INCR_U)
0243 # if FLINT_WANT_ASSERT
0244 #  define MPN_INCR_U(ptr, size, incr) \
0245   do { \
0246     ulong __cy = mpn_add_1(ptr, ptr, size, incr); \
0247     FLINT_ASSERT(__cy == 0); \
0248   } while (0)
0249 #  define MPN_DECR_U(ptr, size, incr) \
0250   do { \
0251     ulong __cy = mpn_sub_1(ptr, ptr, size, incr); \
0252     FLINT_ASSERT(__cy == 0); \
0253   } while (0)
0254 # else
0255 #  define MPN_INCR_U(ptr, size, incr) mpn_add_1(ptr, ptr, size, incr)
0256 #  define MPN_DECR_U(ptr, size, incr) mpn_sub_1(ptr, ptr, size, incr)
0257 # endif
0258 #endif
0259
0260 /* Multiplication */
0261 #if !defined(umul_ppmm)
0262 # define umul_ppmm(w1, w0, u, v) \
0263   do { \
0264     ulong __x0, __x1, __x2, __x3; \
0265     ulong __ul, __vl, __uh, __vh; \
0266     __ul = __ll_lowpart(u); \
0267     __uh = __ll_highpart(u); \
0268     __vl = __ll_lowpart(v); \
0269     __vh = __ll_highpart(v); \
0270     __x0 = __ul * __vl; \
0271     __x1 = __ul * __vh; \
0272     __x2 = __uh * __vl; \
0273     __x3 = __uh * __vh; \
0274     __x1 += __ll_highpart(__x0); /* this can't give carry */ \
0275     __x1 += __x2; /* but this indeed can */ \
0276     if (__x1 < __x2) \
0277       __x3 += __ll_B; \
0278     (w1) = __x3 + __ll_highpart(__x1); \
0279     (w0) = (__x1 << (FLINT_BITS / 2)) + __ll_lowpart(__x0); \
0280   } while (0)
0281
0282 # define smul_ppmm(w1, w0, u, v) \
0283   do { \
0284     ulong __w1, __u = (u), __v = (v); \
0285     umul_ppmm(__w1, w0, u, v); \
0286     (w1) = __w1 - (-(__u >> (FLINT_BITS - 1)) & __v) - (-(__v >> (FLINT_BITS - 1)) & __u); \
0287   } while (0)
0288 #endif
0289
0290 /* Division */
0291 #if !defined(udiv_qrnnd)
0292 # define udiv_qrnnd_int(q, r, n1, n0, d) \
0293   do { \
0294     ulong __d1, __d0, __q1, __q0, __r1, __r0, __m; \
0295     FLINT_ASSERT((d) != 0); \
0296     FLINT_ASSERT((n1) < (d)); \
0297     __d1 = __ll_highpart (d); \
0298     __d0 = __ll_lowpart (d); \
0299     __q1 = (n1) / __d1; \
0300     __r1 = (n1) - __q1 * __d1; \
0301     __m = __q1 * __d0; \
0302     __r1 = __r1 * __ll_B | __ll_highpart (n0); \
0303     if (__r1 < __m) \
0304     { \
0305       __q1--, __r1 += (d); \
0306       if (__r1 >= (d)) /* i.e. we didn't get carry when adding to __r1 */ \
0307         if (__r1 < __m) \
0308           __q1--, __r1 += (d); \
0309     } \
0310     __r1 -= __m; \
0311     __q0 = __r1 / __d1; \
0312     __r0 = __r1  - __q0 * __d1; \
0313     __m = __q0 * __d0; \
0314     __r0 = __r0 * __ll_B | __ll_lowpart (n0); \
0315     if (__r0 < __m) \
0316     { \
0317       __q0--, __r0 += (d); \
0318       if (__r0 >= (d)) \
0319         if (__r0 < __m) \
0320           __q0--, __r0 += (d); \
0321     } \
0322     __r0 -= __m; \
0323     (q) = __q1 * __ll_B | __q0; \
0324     (r) = __r0; \
0325   } while (0)
0326
0327 # define udiv_qrnnd(q, r, n1, n0, d) \
0328   do { \
0329     CLANG_volatile int __norm = flint_clz(d); \
0330     if (__norm != 0) \
0331     { \
0332       udiv_qrnnd_int((q), (r), ((n1) << __norm) + ((n0) >> (FLINT_BITS - __norm)), (n0) << __norm, (d) << __norm); \
0333       (r) = ((ulong) (r) >> __norm); \
0334     } \
0335     else \
0336        udiv_qrnnd_int((q), (r), (n1), (n0), (d)); \
0337   } while (0)
0338
0339 # define __highbit (~(ulong) 0 ^ ((~(ulong) 0) >> 1))
0340
0341 # define sdiv_qrnnd(q, r, n1, n0, d) \
0342   do { \
0343     ulong __n1, __n0, __d; \
0344     ulong __q, __r; \
0345     unsigned int __sgn_n = 0, __sgn_d = 0; \
0346     if ((n1) & __highbit) \
0347     { \
0348       __n0 = -(n0); \
0349       __n1 = ~(n1) + (__n0 == 0); \
0350       __sgn_n = ~__sgn_n; \
0351     } \
0352     else \
0353     { \
0354       __n0 = (n0); \
0355       __n1 = (n1); \
0356     } \
0357     if ((d) & __highbit) \
0358     { \
0359         __d = -(d); \
0360         __sgn_d = ~__sgn_d; \
0361     } \
0362     else \
0363     { \
0364         __d = (d); \
0365     } \
0366     udiv_qrnnd(__q, __r, __n1, __n0, __d); \
0367     q = (__sgn_n == __sgn_d) ? __q : -__q; \
0368     r = (__sgn_n == 0) ? __r : -__r; \
0369   } while (0)
0370 #endif
0371
0372 #define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
0373   do { \
0374     ulong _n2, _n10, _nmask, _nadj, _q1; \
0375     ulong _xh, _xl; \
0376     _n2 = (nh); \
0377     _n10 = (nl); \
0378     _nmask = (slong) (_n10) >> (FLINT_BITS - 1); \
0379     _nadj = _n10 + (_nmask & (d)); \
0380     umul_ppmm(_xh, _xl, di, _n2 - _nmask); \
0381     add_ssaaaa(_xh, _xl, _xh, _xl, _n2, _nadj); \
0382     _q1 = ~_xh; \
0383     umul_ppmm(_xh, _xl, _q1, d); \
0384     add_ssaaaa(_xh, _xl, _xh, _xl, nh, nl); \
0385     _xh -= (d); /* xh = 0 or -1 */ \
0386     (r) = _xl + ((d) & _xh); \
0387     (q) = _xh - _q1; \
0388   } while (0)
0389
0390 #ifdef __cplusplus
0391 }
0392 #endif
0393
0394 #endif