include/flint/radix.h

0001 /*
0002     Copyright (C) 2026 Fredrik Johansson
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 RADIX_H
0013 #define RADIX_H
0014
0015 #ifdef RADIX_INLINES_C
0016 #define RADIX_INLINE
0017 #else
0018 #define RADIX_INLINE static inline
0019 #endif
0020
0021 #include "ulong_extras.h"
0022 #include "mpn_extras.h"
0023 #include "nmod.h"
0024 #include "fmpz_types.h"
0025
0026 #ifdef __cplusplus
0027 extern "C" {
0028 #endif
0029
0030 typedef struct
0031 {
0032     ulong m;
0033     unsigned int e;
0034     unsigned int c;
0035 }
0036 n_div_precomp_struct;
0037
0038 typedef n_div_precomp_struct n_div_precomp_t[1];
0039
0040 typedef struct
0041 {
0042     ulong a;
0043     ulong b;
0044 }
0045 n_pair_struct;
0046
0047 typedef struct
0048 {
0049     nmod_t b;           /* Digit radix */
0050     unsigned int exp;   /* B = b^exp */
0051     nmod_t B;           /* Limb radix */
0052     ulong * bpow;       /* Precomputed powers of b */
0053     n_div_precomp_struct * bpow_div;  /* Division by powers of b */
0054     unsigned char bits_to_digit_size[FLINT_BITS];
0055     unsigned int bval;  /* Binary valuation of b */
0056     n_pair_struct * bpow_oddinv;    /* Inverses of (b>>bval)^k */
0057     ulong (*val_func)(ulong, const void *);  /* Function to compute digit valuation */
0058     /* Extended data which may be used by some ring contexts */
0059     slong trunc_limbs;
0060     slong trunc_digits;
0061 }
0062 radix_struct;
0063
0064 typedef radix_struct radix_t[1];
0065
0066 #define DIGIT_RADIX(radix) ((radix)->b.n)
0067 #define LIMB_RADIX(radix) ((radix)->B.n)
0068
0069 void radix_init(radix_t radix, ulong b, unsigned int exp);
0070 void radix_clear(radix_t radix);
0071 void radix_init_randtest(radix_t radix, flint_rand_t state);
0072 void radix_init_randtest_prime(radix_t radix, flint_rand_t state);
0073
0074 RADIX_INLINE ulong radix_digit_radix(const radix_t radix) { return radix->b.n; }
0075 RADIX_INLINE ulong radix_limb_radix(const radix_t radix) { return radix->B.n; }
0076 RADIX_INLINE ulong radix_limb_exponent(const radix_t radix) { return radix->exp; }
0077
0078 RADIX_INLINE ulong
0079 _radix_size_digits_1(ulong c, const radix_t radix)
0080 {
0081     FLINT_ASSERT(c != 0);
0082     FLINT_ASSERT(c < LIMB_RADIX(radix));
0083
0084     unsigned int bc = FLINT_BIT_COUNT(c);
0085
0086     ulong n = radix->bits_to_digit_size[bc - 1];
0087     n += (c >= radix->bpow[n]);
0088     return n;
0089 }
0090
0091 RADIX_INLINE ulong
0092 radix_size_digits_1(ulong c, const radix_t radix)
0093 {
0094     return (c == 0) ? 0 : _radix_size_digits_1(c, radix);
0095 }
0096
0097 RADIX_INLINE ulong _radix_size_digits(nn_srcptr x, slong n, const radix_t radix)
0098 {
0099     return (n - 1) * radix->exp + _radix_size_digits_1(x[n - 1], radix);
0100 }
0101
0102 RADIX_INLINE ulong radix_size_digits(nn_srcptr x, slong n, const radix_t radix)
0103 {
0104     while (n > 0 && x[n - 1] == 0)
0105         n--;
0106
0107     return (n == 0) ? 0 : _radix_size_digits_1(x[n - 1], radix);
0108 }
0109
0110 RADIX_INLINE ulong _radix_valuation_digits_1(ulong c, const radix_t radix)
0111 {
0112     FLINT_ASSERT(c != 0);
0113     FLINT_ASSERT(c < LIMB_RADIX(radix));
0114
0115     return radix->val_func(c, radix);
0116 }
0117
0118 void radix_rand_limbs(nn_ptr res, flint_rand_t state, slong n, const radix_t radix);
0119 void radix_rand_digits(nn_ptr res, flint_rand_t state, slong n, const radix_t radix);
0120 void radix_randtest_limbs(nn_ptr res, flint_rand_t state, slong n, const radix_t radix);
0121 void radix_randtest_digits(nn_ptr res, flint_rand_t state, slong n, const radix_t radix);
0122
0123 ulong radix_neg(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0124 ulong radix_add(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix);
0125 ulong radix_sub(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix);
0126
0127 /* todo
0128 #define radix_add_n(res, a, b, n, radix) radix_add(res, a, n, b, n, radix)
0129 #define radix_sub_n(res, a, b, n, radix) radix_sub(res, a, n, b, n, radix)
0130
0131 RADIX_INLINE ulong
0132 radix_add_1(nn_ptr res, nn_srcptr a, slong n, ulong c, const radix_t radix)
0133 {
0134     return radix_add(res, a, n, &c, 1, radix);
0135 }
0136
0137 RADIX_INLINE ulong
0138 radix_sub_1(nn_ptr res, nn_srcptr a, slong n, ulong c, const radix_t radix)
0139 {
0140     return radix_sub(res, a, n, &c, 1, radix);
0141 }
0142  */
0143
0144 /* Shift */
0145
0146 ulong radix_lshift_digits(nn_ptr res, nn_srcptr a, slong n, unsigned int e, const radix_t radix);
0147 ulong radix_rshift_digits(nn_ptr res, nn_srcptr a, slong n, unsigned int e, const radix_t radix);
0148
0149 /* fused add/sub with an on-the-fly left shift by sh digits, 1 <= sh < exp */
0150 ulong radix_addlsh(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, unsigned int sh, const radix_t radix);
0151 ulong radix_sublsh(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, unsigned int sh, const radix_t radix);
0152 ulong radix_lshsub(nn_ptr res, nn_srcptr b, slong bn, unsigned int sh, nn_srcptr a, slong an, const radix_t radix);
0153 int radix_cmplsh(nn_srcptr a, slong an, nn_srcptr b, slong bn, unsigned int sh, const radix_t radix);
0154
0155 /* Multiplication */
0156
0157 ulong radix_mul_1(nn_ptr res, nn_srcptr a, slong n, ulong c, const radix_t radix);
0158
0159 void radix_mulmid_fft_small(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong lo, slong hi, const radix_t radix);
0160 void radix_mulmid_classical(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong lo, slong hi, const radix_t radix);
0161 void radix_mulmid_KS(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong lo, slong hi, const radix_t radix);
0162 void radix_mulmid_naive(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong lo, slong hi, const radix_t radix);
0163
0164 RADIX_INLINE void
0165 radix_mulmid(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong lo, slong hi, const radix_t radix)
0166 {
0167     /* todo: tuning */
0168     if (FLINT_MIN(an, bn) < 80 || hi - lo < 80)
0169         radix_mulmid_classical(res, a, an, b, bn, lo, hi, radix);
0170     else
0171 #if FLINT_HAVE_FFT_SMALL
0172         radix_mulmid_fft_small(res, a, an, b, bn, lo, hi, radix);
0173 #else
0174         radix_mulmid_KS(res, a, an, b, bn, lo, hi, radix);
0175 #endif
0176 }
0177
0178 RADIX_INLINE void
0179 radix_mul(nn_ptr res, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix)
0180 {
0181     radix_mulmid(res, a, an, b, bn, 0, an + bn, radix);
0182 }
0183
0184 /* todo: squaring optimisations in all multiplication algorithms */
0185 RADIX_INLINE void
0186 radix_sqr(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)
0187 {
0188     radix_mul(res, a, an, a, an, radix);
0189 }
0190
0191 RADIX_INLINE ulong
0192 radix_mul_two(nn_ptr res, nn_srcptr a, slong an, const radix_t radix)
0193 {
0194     return radix_add(res, a, an, a, an, radix);
0195 }
0196
0197 void
0198 _radix_mulhigh_known_low(nn_ptr out, nn_srcptr x, slong xn, nn_srcptr y, slong yn,
0199     nn_srcptr kl, slong kl_len, slong klo, slong khi, nn_ptr scratch,
0200     const radix_t radix);
0201
0202 /* Division */
0203
0204 ulong radix_divrem_1(nn_ptr res, nn_srcptr a, slong an, ulong d, const radix_t radix);
0205 void radix_divexact_1(nn_ptr res, nn_srcptr a, slong an, ulong d, const radix_t radix);
0206 ulong radix_divrem_two(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0207
0208 void radix_inv_approx_basecase(nn_ptr q, nn_srcptr a, slong an, slong n, const radix_t radix);
0209 void radix_inv_approx(nn_ptr q, nn_srcptr a, slong an, slong n, const radix_t radix);
0210 void radix_div_approx_invmul(nn_ptr Q, nn_srcptr B, slong Bn, nn_srcptr A, slong An, slong n, const radix_t radix);
0211 void radix_div_approx(nn_ptr Q, nn_srcptr B, slong Bn, nn_srcptr A, slong An, slong n, const radix_t radix);
0212
0213 void radix_divrem_preinv(nn_ptr q, nn_ptr r, nn_srcptr a, slong an, nn_srcptr b, slong bn, nn_srcptr binv, slong binvn, const radix_t radix);
0214 void radix_divrem_via_mpn(nn_ptr q, nn_ptr r, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix);
0215 void radix_divrem_newton(nn_ptr q, nn_ptr r, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix);
0216 void radix_divrem_newton_karp_markstein(nn_ptr Q, nn_ptr R, nn_srcptr A, slong An, nn_srcptr B, slong Bn, const radix_t radix);
0217 void radix_divrem(nn_ptr q, nn_ptr r, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix);
0218
0219 int radix_div(nn_ptr q, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix);
0220 void radix_divexact(nn_ptr q, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix);
0221
0222 /* compare (x, n) with floor(B^n / 2) */
0223 RADIX_INLINE int
0224 radix_cmp_bn_half(nn_srcptr x, slong n, const radix_t radix)
0225 {
0226     FLINT_ASSERT(n >= 1);
0227
0228     ulong B = LIMB_RADIX(radix);
0229     ulong d, B2;
0230     slong i;
0231
0232     d = x[n - 1];
0233     B2 = B / 2;
0234     if (d != B2)
0235         return (d < B2 ? -1 : 1);
0236
0237     B2 = (B % 2) ? B2 : 0;
0238     for (i = n - 2; i >= 0; i--)
0239     {
0240         d = x[i];
0241         if (d != B2)
0242             return (d < B2 ? -1 : 1);
0243     }
0244
0245     return 0;
0246 }
0247
0248 /* Square roots */
0249
0250 void radix_rsqrt_1_approx_basecase(nn_ptr res, ulong a, slong n, const radix_t radix);
0251 void radix_rsqrt_1_approx(nn_ptr res, ulong a, slong n, const radix_t radix);
0252
0253 /* Modular (Hensel) division */
0254
0255 int radix_invmod_bn(nn_ptr res, nn_srcptr x, slong xn, slong n, const radix_t radix);
0256
0257 int radix_divmod_bn_1(nn_ptr q, nn_ptr rem, nn_srcptr a, slong an, ulong b, slong n, const radix_t radix);
0258 int radix_divmod_bn_classical(nn_ptr q, nn_ptr rem, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong n, const radix_t radix);
0259 int radix_divmod_bn_karp_markstein(nn_ptr q, nn_ptr rem, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong n, const radix_t radix);
0260 int radix_divmod_bn(nn_ptr q, nn_ptr rem, nn_srcptr a, slong an, nn_srcptr b, slong bn, slong n, const radix_t radix);
0261
0262 /* Modular square roots */
0263
0264 int radix_rsqrtmod_bn(nn_ptr res, nn_srcptr x, slong xn, slong n, const radix_t radix);
0265 int radix_sqrtmod_bn(nn_ptr res, nn_srcptr x, slong xn, slong n, const radix_t radix);
0266
0267 /* Radix conversion */
0268
0269 typedef struct
0270 {
0271     slong len;
0272     ulong exps[FLINT_BITS];
0273     nn_ptr pows[FLINT_BITS];
0274     slong sizes[FLINT_BITS];
0275     slong val_limbs[FLINT_BITS];
0276     nn_ptr buf;
0277 }
0278 radix_powers_struct;
0279
0280 typedef radix_powers_struct radix_powers_t[1];
0281
0282 void radix_powers_clear(radix_powers_t powers);
0283
0284 slong radix_get_mpn_basecase(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0285 slong radix_get_mpn_divconquer(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0286 slong radix_get_mpn(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0287
0288 slong radix_set_mpn_basecase(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0289 slong radix_set_mpn_divconquer(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0290 slong radix_set_mpn(nn_ptr res, nn_srcptr a, slong an, const radix_t radix);
0291
0292 slong radix_set_mpn_need_alloc(slong n, const radix_t radix);
0293
0294 /* String conversion */
0295
0296 char * radix_get_str_decimal(char * res, nn_srcptr x, slong n, int negative, const radix_t radix);
0297 char * radix_get_str_sum(char * res, nn_srcptr x, slong n, int negative, int ascending, const radix_t radix);
0298
0299 /* Memory-managed integers */
0300
0301 typedef struct
0302 {
0303     nn_ptr d;
0304     slong alloc;
0305     slong size;
0306 }
0307 radix_integer_struct;
0308
0309 typedef radix_integer_struct radix_integer_t[1];
0310
0311 void gr_ctx_init_radix_integer(gr_ctx_t ctx, ulong b, unsigned int exp);
0312
0313 void radix_integer_init(radix_integer_t res, const radix_t radix);
0314 void radix_integer_clear(radix_integer_t res, const radix_t radix);
0315 nn_ptr radix_integer_fit_limbs(radix_integer_t res, slong nlimbs, const radix_t radix);
0316 void radix_integer_zero(radix_integer_t res, const radix_t radix);
0317 void radix_integer_rand_limbs(radix_integer_t res, flint_rand_t state, slong max_limbs, const radix_t radix);
0318 void radix_integer_randtest_limbs(radix_integer_t res, flint_rand_t state, slong max_limbs, const radix_t radix);
0319 void radix_integer_one(radix_integer_t res, const radix_t radix);
0320 void radix_integer_neg_one(radix_integer_t res, const radix_t radix);
0321 int radix_integer_is_zero(const radix_integer_t x, const radix_t radix);
0322 int radix_integer_is_one(const radix_integer_t x, const radix_t radix);
0323 int radix_integer_is_neg_one(const radix_integer_t x, const radix_t radix);
0324 int radix_integer_equal(const radix_integer_t x, const radix_integer_t y, const radix_t radix);
0325 int radix_integer_cmp(const radix_integer_t x, const radix_integer_t y, const radix_t radix);
0326 int radix_integer_cmpabs(const radix_integer_t x, const radix_integer_t y, const radix_t radix);
0327 void radix_integer_set(radix_integer_t res, const radix_integer_t x, const radix_t radix);
0328 void radix_integer_set_ui(radix_integer_t res, ulong x, const radix_t radix);
0329 void radix_integer_set_si(radix_integer_t res, slong x, const radix_t radix);
0330 void radix_integer_set_fmpz(radix_integer_t res, const fmpz_t x, const radix_t radix);
0331 void radix_integer_get_fmpz(fmpz_t res, const radix_integer_t x, const radix_t radix);
0332 void radix_integer_neg(radix_integer_t res, const radix_integer_t x, const radix_t radix);
0333 void radix_integer_abs(radix_integer_t res, const radix_integer_t x, const radix_t radix);
0334 int radix_integer_sgn(const radix_integer_t x, const radix_t radix);
0335 void radix_integer_add(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, const radix_t radix);
0336 void radix_integer_sub(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, const radix_t radix);
0337 void radix_integer_mul(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, const radix_t radix);
0338
0339 int radix_integer_is_normalised(const radix_integer_t x, const radix_t radix);
0340
0341 RADIX_INLINE slong
0342 radix_integer_size_limbs(const radix_integer_t x, const radix_t FLINT_UNUSED(radix))
0343 {
0344     return FLINT_ABS(x->size);
0345 }
0346
0347 RADIX_INLINE slong
0348 radix_integer_ssize_limbs(const radix_integer_t x, const radix_t FLINT_UNUSED(radix))
0349 {
0350     return x->size;
0351 }
0352
0353 RADIX_INLINE slong
0354 radix_integer_size_digits(const radix_integer_t x, const radix_t radix)
0355 {
0356     return (x->size == 0) ? 0 : _radix_size_digits(x->d, FLINT_ABS(x->size), radix);
0357 }
0358
0359 RADIX_INLINE slong
0360 radix_integer_ssize_digits(const radix_integer_t x, const radix_t radix)
0361 {
0362     if (x->size == 0) return 0;
0363     if (x->size > 0) return _radix_size_digits(x->d, x->size, radix);
0364     return -_radix_size_digits(x->d, -x->size, radix);
0365 }
0366
0367 RADIX_INLINE ulong
0368 radix_integer_get_limb(const radix_integer_t x, slong n, const radix_t FLINT_UNUSED(radix))
0369 {
0370     FLINT_ASSERT(n >= 0);
0371     return (n >= FLINT_ABS(x->size)) ? 0 : x->d[n];
0372 }
0373
0374 void radix_integer_set_limb(radix_integer_t res, const radix_integer_t x, slong index, ulong c, const radix_t radix);
0375
0376 ulong radix_integer_get_digit(const radix_integer_t x, slong index, const radix_t radix);
0377 void radix_integer_set_digit(radix_integer_t res, const radix_integer_t x, slong index, ulong c, const radix_t radix);
0378
0379 void radix_integer_lshift_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0380 void radix_integer_rshift_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0381 void radix_integer_lshift_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0382 void radix_integer_rshift_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0383
0384 void radix_integer_addlsh(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, slong v, const radix_t radix);
0385 void radix_integer_sublsh(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, slong v, const radix_t radix);
0386
0387 RADIX_INLINE slong
0388 radix_integer_valuation_limbs(const radix_integer_t x, const radix_t FLINT_UNUSED(radix))
0389 {
0390     slong xn = FLINT_ABS(x->size);
0391     if (xn == 0)
0392         return 0;
0393     slong v = 0;
0394     while (x->d[v] == 0)
0395         v++;
0396     return v;
0397 }
0398
0399 RADIX_INLINE slong
0400 radix_integer_valuation_digits(const radix_integer_t x, const radix_t radix)
0401 {
0402     slong xn = FLINT_ABS(x->size);
0403     if (xn == 0)
0404         return 0;
0405     slong v = 0;
0406     while (x->d[v] == 0)
0407         v++;
0408     return v * radix->exp + _radix_valuation_digits_1(x->d[v], radix);
0409 }
0410
0411 void radix_integer_trunc_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0412 void radix_integer_mod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0413 void radix_integer_smod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0414
0415 void radix_integer_trunc_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0416 void radix_integer_mod_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0417
0418 void radix_integer_mullow_limbs(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, slong n, const radix_t radix);
0419 int radix_integer_invmod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0420 int radix_integer_divmod_limbs(radix_integer_t res, const radix_integer_t a, const radix_integer_t b, slong n, const radix_t radix);
0421 int radix_integer_rsqrtmod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0422 int radix_integer_sqrtmod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix);
0423
0424 int radix_integer_div(radix_integer_t q, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0425 void radix_integer_divexact(radix_integer_t q, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0426
0427 void radix_integer_tdiv_qr(radix_integer_t q, radix_integer_t r, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0428 void radix_integer_fdiv_qr(radix_integer_t q, radix_integer_t r, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0429 void radix_integer_cdiv_qr(radix_integer_t q, radix_integer_t r, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0430 void radix_integer_tdiv_q(radix_integer_t q, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0431 void radix_integer_fdiv_q(radix_integer_t q, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0432 void radix_integer_cdiv_q(radix_integer_t q, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0433 void radix_integer_tdiv_r(radix_integer_t r, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0434 void radix_integer_fdiv_r(radix_integer_t r, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0435 void radix_integer_cdiv_r(radix_integer_t r, const radix_integer_t a, const radix_integer_t b, const radix_t radix);
0436
0437
0438 /* Utilities */
0439
0440 RADIX_INLINE void
0441 n_div_precomp_init(n_div_precomp_t pre, ulong d)
0442 {
0443     pre->e = FLINT_BIT_COUNT(d) - 1;
0444
0445     if ((d & (d - 1)) == 0)
0446     {
0447         pre->m = 0;
0448         pre->c = 0;
0449     }
0450     else
0451     {
0452         ulong q, r;
0453         udiv_qrnnd(q, r, UWORD(1) << pre->e, 0, d);
0454         if (d - r < (UWORD(1) << pre->e))
0455         {
0456             pre->m = q + 1;
0457             pre->c = 0;
0458         }
0459         else
0460         {
0461             pre->m = q;
0462             pre->c = 1;
0463         }
0464     }
0465 }
0466
0467 RADIX_INLINE ulong
0468 n_div_precomp_m0(ulong x, const n_div_precomp_t pre)
0469 {
0470     return x >> pre->e;
0471 }
0472
0473 RADIX_INLINE ulong
0474 n_div_precomp_c0(ulong x, const n_div_precomp_t pre)
0475 {
0476     return n_mulhi(x, pre->m) >> pre->e;
0477 }
0478
0479 RADIX_INLINE ulong n_incsat(ulong x)
0480 {
0481     ulong t = x + 1;
0482     return (t > x) ? t : x;
0483 }
0484
0485 /* Assumes x != UWORD_MAX */
0486 RADIX_INLINE ulong
0487 n_div_precomp_c1_bounded(ulong x, const n_div_precomp_t pre)
0488 {
0489     return n_mulhi(x + 1, pre->m) >> pre->e;
0490 }
0491
0492 RADIX_INLINE ulong
0493 n_div_precomp_c1(ulong x, const n_div_precomp_t pre)
0494 {
0495     return n_mulhi(n_incsat(x), pre->m) >> pre->e;
0496 }
0497
0498 RADIX_INLINE ulong
0499 n_div_precomp(ulong x, const n_div_precomp_t pre)
0500 {
0501     if (pre->m == 0)
0502         return x >> pre->e;
0503     else if (pre->c == 0)
0504         return n_mulhi(x, pre->m) >> pre->e;
0505     else
0506         return n_mulhi(n_incsat(x), pre->m) >> pre->e;
0507 }
0508
0509 RADIX_INLINE ulong
0510 n_div_precomp_bounded(ulong x, const n_div_precomp_t pre)
0511 {
0512     if (pre->m == 0)
0513         return x >> pre->e;
0514     else if (pre->c == 0)
0515         return n_mulhi(x, pre->m) >> pre->e;
0516     else
0517         return n_mulhi(x + 1, pre->m) >> pre->e;
0518 }
0519
0520 RADIX_INLINE ulong
0521 n_rem_precomp_m0(ulong x, ulong FLINT_UNUSED(d), const n_div_precomp_t pre)
0522 {
0523     return x & ((UWORD(1) << pre->e) - 1);
0524 }
0525
0526 RADIX_INLINE ulong
0527 n_rem_precomp_c0(ulong x, ulong d, const n_div_precomp_t pre)
0528 {
0529     return x - d * n_div_precomp_c0(x, pre);
0530 }
0531
0532 RADIX_INLINE ulong
0533 n_rem_precomp_c1(ulong x, ulong d, const n_div_precomp_t pre)
0534 {
0535     return x - d * n_div_precomp_c1(x, pre);
0536 }
0537
0538 RADIX_INLINE ulong
0539 n_rem_precomp_bounded(ulong x, ulong d, const n_div_precomp_t pre)
0540 {
0541     return x - d * n_div_precomp_bounded(x, pre);
0542 }
0543
0544 RADIX_INLINE ulong
0545 n_rem_precomp_c1_bounded(ulong x, ulong d, const n_div_precomp_t pre)
0546 {
0547     return x - d * n_div_precomp_c1_bounded(x, pre);
0548 }
0549
0550 RADIX_INLINE ulong
0551 n_rem_precomp(ulong x, ulong d, const n_div_precomp_t pre)
0552 {
0553     return x - d * n_div_precomp(x, pre);
0554 }
0555
0556 RADIX_INLINE ulong
0557 n_divrem_precomp_m0(ulong * r, ulong x, ulong d, const n_div_precomp_t pre)
0558 {
0559     ulong q = n_div_precomp_m0(x, pre);
0560     *r = x - d * q;
0561     return q;
0562 }
0563
0564 RADIX_INLINE ulong
0565 n_divrem_precomp_c0(ulong * r, ulong x, ulong d, const n_div_precomp_t pre)
0566 {
0567     ulong q = n_div_precomp_c0(x, pre);
0568     *r = x - d * q;
0569     return q;
0570 }
0571
0572 RADIX_INLINE ulong
0573 n_divrem_precomp_c1(ulong * r, ulong x, ulong d, const n_div_precomp_t pre)
0574 {
0575     ulong q = n_div_precomp_c1(x, pre);
0576     *r = x - d * q;
0577     return q;
0578 }
0579
0580 RADIX_INLINE ulong
0581 n_divrem_precomp_c1_bounded(ulong * r, ulong x, ulong d, const n_div_precomp_t pre)
0582 {
0583     ulong q = n_div_precomp_c1_bounded(x, pre);
0584     *r = x - d * q;
0585     return q;
0586 }
0587
0588 RADIX_INLINE ulong
0589 n_divrem_precomp(ulong * r, ulong x, ulong d, const n_div_precomp_t pre)
0590 {
0591     ulong q = n_div_precomp(x, pre);
0592     *r = x - d * q;
0593     return q;
0594 }
0595
0596 RADIX_INLINE ulong
0597 n_divrem_precomp_bounded(ulong * r, ulong x, ulong d, const n_div_precomp_t pre)
0598 {
0599     ulong q = n_div_precomp_bounded(x, pre);
0600     *r = x - d * q;
0601     return q;
0602 }
0603
0604
0605 #ifdef __cplusplus
0606 }
0607 #endif
0608
0609 #endif