include/flint/nmod.h

0001 /*
0002     Copyright (C) 2010 William Hart
0003     Copyright (C) 2021 Fredrik Johansson
0004
0005     This file is part of FLINT.
0006
0007     FLINT is free software: you can redistribute it and/or modify it under
0008     the terms of the GNU Lesser General Public License (LGPL) as published
0009     by the Free Software Foundation; either version 3 of the License, or
0010     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0011 */
0012
0013 #ifndef NMOD_H
0014 #define NMOD_H
0015
0016 #ifdef NMOD_INLINES_C
0017 #define NMOD_INLINE
0018 #else
0019 #define NMOD_INLINE static inline
0020 #endif
0021
0022 #include "ulong_extras.h"
0023 #include "nmod_types.h"
0024 #include "gr_types.h"
0025
0026 #ifdef __cplusplus
0027 extern "C" {
0028 #endif
0029
0030 #define NMOD_RED2(r, a_hi, a_lo, mod) \
0031   do { \
0032     ulong q0xx, q1xx, r1xx; \
0033     const ulong u1xx = ((a_hi)<<(mod).norm) \
0034      + (((mod).norm == 0) ? UWORD(0) : (a_lo)>>(FLINT_BITS - (mod).norm)); \
0035     const ulong u0xx = (a_lo)<<(mod).norm; \
0036     const ulong nxx = (mod).n<<(mod).norm; \
0037     umul_ppmm(q1xx, q0xx, (mod).ninv, u1xx); \
0038     add_ssaaaa(q1xx, q0xx, q1xx, q0xx, u1xx, u0xx); \
0039     r1xx = (u0xx - (q1xx + 1)*nxx); \
0040     if (r1xx > q0xx) r1xx += nxx; \
0041     if (r1xx < nxx) r = (r1xx>>(mod).norm); \
0042     else r = ((r1xx - nxx)>>(mod).norm); \
0043   } while (0)
0044
0045 #define NMOD_RED(r, a, mod) \
0046   do { \
0047     NMOD_RED2(r, UWORD(0), a, mod); \
0048   } while (0)
0049
0050 #define NMOD2_RED2(r, a_hi, a_lo, mod) \
0051   do { \
0052     ulong v_hi; \
0053     NMOD_RED(v_hi, a_hi, mod); \
0054     NMOD_RED2(r, v_hi, a_lo, mod); \
0055   } while (0)
0056
0057 #define NMOD_RED3(r, a_hi, a_me, a_lo, mod) \
0058   do { \
0059     ulong v_hi; \
0060     NMOD_RED2(v_hi, a_hi, a_me, mod); \
0061     NMOD_RED2(r, v_hi, a_lo, mod); \
0062   } while (0)
0063
0064 #define NMOD_BITS(mod) (FLINT_BITS - ((mod).norm))
0065 #define NMOD_CAN_USE_SHOUP(mod) ((mod).norm > 0)
0066
0067 #define NMOD_RED2_FULLWORD(r, a_hi, a_lo, mod) \
0068   do { \
0069     ulong q0xx, q1xx, r1xx; \
0070     const ulong u1xx = (a_hi); \
0071     const ulong u0xx = (a_lo); \
0072     const ulong nxx = (mod).n; \
0073     umul_ppmm(q1xx, q0xx, (mod).ninv, u1xx); \
0074     add_ssaaaa(q1xx, q0xx, q1xx, q0xx, u1xx, u0xx); \
0075     r1xx = (u0xx - (q1xx + 1)*nxx); \
0076     if (r1xx > q0xx) r1xx += nxx; \
0077     if (r1xx < nxx) r = (r1xx); \
0078     else r = ((r1xx - nxx)); \
0079   } while (0)
0080
0081 #define NMOD_RED2_NONFULLWORD(r, a_hi, a_lo, mod) \
0082   do { \
0083     ulong q0xx, q1xx, r1xx; \
0084     const ulong u1xx = ((a_hi)<<(mod).norm) + ((a_lo)>>(FLINT_BITS - (mod).norm)); \
0085     const ulong u0xx = (a_lo)<<(mod).norm; \
0086     const ulong nxx = (mod).n<<(mod).norm; \
0087     umul_ppmm(q1xx, q0xx, (mod).ninv, u1xx); \
0088     add_ssaaaa(q1xx, q0xx, q1xx, q0xx, u1xx, u0xx); \
0089     r1xx = (u0xx - (q1xx + 1)*nxx); \
0090     if (r1xx > q0xx) r1xx += nxx; \
0091     if (r1xx < nxx) r = (r1xx>>(mod).norm); \
0092     else r = ((r1xx - nxx)>>(mod).norm); \
0093   } while (0)
0094
0095 #define NMOD_MUL_PRENORM(res, a, b, mod) \
0096   do { \
0097     ulong q0xx, q1xx, rxx, p_hixx, p_loxx; \
0098     ulong nxx, ninvxx; \
0099     unsigned int normxx; \
0100     ninvxx = (mod).ninv; \
0101     normxx = (mod).norm; \
0102     nxx = (mod).n << normxx; \
0103     umul_ppmm(p_hixx, p_loxx, (a), (b)); \
0104     umul_ppmm(q1xx, q0xx, ninvxx, p_hixx); \
0105     add_ssaaaa(q1xx, q0xx, q1xx, q0xx, p_hixx, p_loxx); \
0106     rxx = (p_loxx - (q1xx + 1) * nxx); \
0107     if (rxx > q0xx) \
0108       rxx += nxx; \
0109     rxx = (rxx < nxx ? rxx : rxx - nxx) >> normxx; \
0110     (res) = rxx; \
0111   } while (0)
0112
0113 #define NMOD_MUL_FULLWORD(res, a, b, mod) \
0114   do { \
0115     ulong q0xx, q1xx, rxx, p_hixx, p_loxx; \
0116     ulong nxx, ninvxx; \
0117     ninvxx = (mod).ninv; \
0118     nxx = (mod).n; \
0119     umul_ppmm(p_hixx, p_loxx, (a), (b)); \
0120     umul_ppmm(q1xx, q0xx, ninvxx, p_hixx); \
0121     add_ssaaaa(q1xx, q0xx, q1xx, q0xx, p_hixx, p_loxx); \
0122     rxx = (p_loxx - (q1xx + 1) * nxx); \
0123     if (rxx > q0xx) \
0124       rxx += nxx; \
0125     rxx = (rxx < nxx ? rxx : rxx - nxx); \
0126     (res) = rxx; \
0127   } while (0)
0128
0129 NMOD_INLINE ulong nmod_set_ui(ulong x, nmod_t mod)
0130 {
0131     if (x < mod.n)
0132         return x;
0133
0134     NMOD_RED(x, x, mod);
0135     return x;
0136 }
0137
0138 NMOD_INLINE
0139 ulong nmod_set_si(slong x, nmod_t mod)
0140 {
0141     ulong res = (x >= 0) ? (ulong) x : -(ulong) x;
0142     NMOD_RED(res, res, mod);
0143     return (res == 0 || x > 0) ? res : mod.n - res;
0144 }
0145
0146 NMOD_INLINE
0147 ulong _nmod_add(ulong a, ulong b, nmod_t mod)
0148 {
0149    const ulong sum = a + b;
0150    return sum - mod.n + ((((slong)(sum - mod.n))>>(FLINT_BITS - 1)) & mod.n);
0151 }
0152
0153 NMOD_INLINE
0154 ulong _nmod_sub(ulong a, ulong b, nmod_t mod)
0155 {
0156    const ulong diff = a - b;
0157    return  ((((slong)diff)>>(FLINT_BITS - 1)) & mod.n) + diff;
0158 }
0159
0160 NMOD_INLINE
0161 ulong nmod_add(ulong a, ulong b, nmod_t mod)
0162 {
0163     const ulong neg = mod.n - a;
0164     FLINT_ASSERT(a < mod.n);
0165     FLINT_ASSERT(b < mod.n);
0166     if (neg > b)
0167         return a + b;
0168     else
0169         return b - neg;
0170 }
0171
0172 NMOD_INLINE
0173 ulong nmod_ui_add_ui(ulong a, ulong b, nmod_t mod)
0174 {
0175     return nmod_add(nmod_set_ui(a, mod), nmod_set_ui(b, mod), mod);
0176 }
0177
0178 NMOD_INLINE
0179 ulong nmod_sub(ulong a, ulong b, nmod_t mod)
0180 {
0181     const ulong diff = a - b;
0182     FLINT_ASSERT(a < mod.n);
0183     FLINT_ASSERT(b < mod.n);
0184     if (a < b)
0185         return mod.n + diff;
0186     else
0187         return diff;
0188 }
0189
0190 NMOD_INLINE
0191 ulong nmod_neg(ulong a, nmod_t mod)
0192 {
0193    if (a)
0194       return mod.n - a;
0195    else
0196       return 0;
0197 }
0198
0199 NMOD_INLINE
0200 ulong nmod_mul(ulong a, ulong b, nmod_t mod)
0201 {
0202     ulong res;
0203     FLINT_ASSERT(a < mod.n);
0204     FLINT_ASSERT(b < mod.n);
0205     NMOD_MUL_PRENORM(res, a, b << mod.norm, mod);
0206     return res;
0207 }
0208
0209 NMOD_INLINE ulong
0210 nmod_ui_mul_ui(ulong a, ulong b, nmod_t mod)
0211 {
0212     return n_mulmod2_preinv(a, b, mod.n, mod.ninv);
0213 }
0214
0215 NMOD_INLINE
0216 ulong _nmod_mul_fullword(ulong a, ulong b, nmod_t mod)
0217 {
0218     ulong res;
0219     NMOD_MUL_FULLWORD(res, a, b, mod);
0220     return res;
0221 }
0222
0223 NMOD_INLINE
0224 ulong nmod_addmul(ulong s, ulong a, ulong b, nmod_t mod)
0225 {
0226     ulong hi, lo;
0227     FLINT_ASSERT(s < mod.n);
0228     FLINT_ASSERT(a < mod.n);
0229     FLINT_ASSERT(b < mod.n);
0230     umul_ppmm(hi, lo, a, b);
0231     add_ssaaaa(hi, lo, hi, lo, 0, s);
0232     NMOD_RED2(lo, hi, lo, mod);
0233     return lo;
0234 }
0235
0236 #define NMOD_ADDMUL(r, a, b, mod) \
0237     do { \
0238        (r) = nmod_addmul((r), (a), (b), (mod)); \
0239     } while (0)
0240
0241 // TODO doc  a*b + c*d
0242 NMOD_INLINE
0243 ulong nmod_fmma(ulong a, ulong b, ulong c, ulong d, nmod_t mod)
0244 {
0245     a = nmod_mul(a, b, mod);
0246     NMOD_ADDMUL(a, c, d, mod);
0247     return a;
0248 }
0249
0250 NMOD_INLINE
0251 ulong nmod_inv(ulong a, nmod_t mod)
0252 {
0253     return n_invmod(a, mod.n);
0254 }
0255
0256 NMOD_INLINE
0257 ulong nmod_div(ulong a, ulong b, nmod_t mod)
0258 {
0259     return nmod_mul(a, n_invmod(b, mod.n), mod);
0260 }
0261
0262 int nmod_divides(ulong * a, ulong b, ulong c, nmod_t mod);
0263
0264 ulong _nmod_pow_ui_redc(ulong a, ulong exp, nmod_t mod);
0265 ulong _nmod_pow_ui_binexp(ulong a, ulong exp, nmod_t mod);
0266 ulong nmod_pow_ui(ulong a, ulong exp, nmod_t mod);
0267
0268 NMOD_INLINE ulong
0269 nmod_ui_pow_ui(ulong a, ulong exp, nmod_t mod)
0270 {
0271     return nmod_pow_ui(nmod_set_ui(a, mod), exp, mod);
0272 }
0273
0274 ulong _nmod_2_pow_ui_binexp(ulong exp, nmod_t mod);
0275 ulong nmod_2_pow_ui(ulong exp, nmod_t mod);
0276
0277 NMOD_INLINE
0278 ulong nmod_pow_fmpz(ulong a, const fmpz_t exp, nmod_t mod)
0279 {
0280     return n_powmod2_fmpz_preinv(a, exp, mod.n, mod.ninv);
0281 }
0282
0283 NMOD_INLINE
0284 void nmod_init(nmod_t * mod, ulong n)
0285 {
0286     mod->n = n;
0287     mod->norm = flint_clz(n);
0288     mod->ninv = n_preinvert_limb_prenorm(n << (mod->norm));
0289 }
0290
0291 /* Montgomery arithmetic *****************************************************/
0292
0293 /* Internal helpers for Montgomery arithmetic ********************************/
0294 /* Some of these functions could be moved to ulong_extras in the future. */
0295
0296 /* Experimental: see if the compiler generates better code than inline assembly */
0297 #if FLINT_BITS == 64 && defined(__GNUC__)
0298
0299 typedef __uint128_t ull_t;
0300
0301 NMOD_INLINE ulong ull_hi(ull_t x) { return x >> FLINT_BITS; }
0302 NMOD_INLINE ulong ull_lo(ull_t x) { return (ulong) x; }
0303 NMOD_INLINE ull_t ull_add(ull_t x, ull_t y) { return x + y; }
0304 NMOD_INLINE ull_t ull_add_u(ull_t x, ulong y) { return x + (ull_t) y; }
0305 NMOD_INLINE ull_t ull_u_mul_u(ulong x, ulong y) { return (ull_t) x * (ull_t) y; }
0306 NMOD_INLINE ull_t ull(ulong hi, ulong lo) { return ((ull_t) lo) | (((ull_t) hi) << FLINT_BITS); }
0307
0308 #else
0309
0310 typedef struct { ulong lo; ulong hi; } ull_t;
0311
0312 NMOD_INLINE ulong ull_hi(ull_t x) { return x.hi; }
0313 NMOD_INLINE ulong ull_lo(ull_t x) { return x.lo; }
0314 NMOD_INLINE ull_t ull_add(ull_t x, ull_t y) { ull_t t; add_ssaaaa(t.hi, t.lo, x.hi, x.lo, y.hi, y.lo); return t; }
0315 NMOD_INLINE ull_t ull_add_u(ull_t x, ulong y) { ull_t t; add_ssaaaa(t.hi, t.lo, x.hi, x.lo, 0, y); return t; }
0316 NMOD_INLINE ull_t ull_u_mul_u(ulong x, ulong y) { ull_t t; umul_ppmm(t.hi, t.lo, x, y); return t; }
0317 NMOD_INLINE ull_t ull(ulong hi, ulong lo) { ull_t t; t.lo = lo; t.hi = hi; return t; }
0318
0319 #endif
0320
0321 /* Assumes a already reduced mod n. */
0322 FLINT_FORCE_INLINE ulong
0323 n_to_redc_preinv(ulong a, ulong n, ulong ninv, ulong norm)
0324 {
0325     ulong q0, q1, r;
0326     FLINT_ASSERT(a < n);
0327     n <<= norm;
0328     a <<= norm;
0329     ull_t t = ull_u_mul_u(ninv, a);
0330     q1 = ull_hi(t);
0331     q0 = ull_lo(t);
0332     q1 += a;
0333     r = -(q1 + 1) * n;
0334     if (r > q0)
0335         r += n;
0336     return (r < n) ? (r >> norm) : ((r - n) >> norm);
0337 }
0338
0339 /* Computes x/R, reducing to [0,2n). */
0340 /* nred = -1/n mod R, R = 2^FLINT_BITS */
0341 FLINT_FORCE_INLINE
0342 ulong n_redc_fast(ulong x, ulong n, ulong nred)
0343 {
0344     return ull_hi(ull_add_u(ull_u_mul_u(n, x * nred), x));
0345 }
0346
0347 FLINT_FORCE_INLINE
0348 ulong n_ll_redc_fast(ull_t x, ulong n, ulong nred)
0349 {
0350     return ull_hi(ull_add(ull_u_mul_u(n, ull_lo(x) * nred), x));
0351 }
0352
0353 /* Computes x/R mod n, reducing to [0,n). */
0354 FLINT_FORCE_INLINE
0355 ulong n_redc(ulong x, ulong n, ulong nred)
0356 {
0357     ulong y = n_redc_fast(x, n, nred);
0358     if (y >= n)
0359         y -= n;
0360     return y;
0361 }
0362
0363 /* Computes x/R, reducing to [0,n). Assumes x < n * 2^FLINT_BITS. */
0364 FLINT_FORCE_INLINE
0365 ulong n_ll_redc(ull_t x, ulong n, ulong nred)
0366 {
0367     ulong lo = n_mulhi(ull_lo(x) * (-nred), n), hi;
0368     hi = ull_hi(x);
0369
0370     if (hi < lo)
0371         return hi - lo + n;
0372     else
0373         return hi - lo;
0374 }
0375
0376 /* Computes (ab)/R mod n, reducing to [0,n). */
0377 FLINT_FORCE_INLINE
0378 ulong n_mulmod_redc(ulong a, ulong b, ulong n, ulong nred)
0379 {
0380     return n_ll_redc(ull_u_mul_u(a, b), n, nred);
0381 }
0382
0383 /* Computes (ab)/R mod n, reducing to [0,2n). */
0384 FLINT_FORCE_INLINE
0385 ulong n_mulmod_redc_fast(ulong a, ulong b, ulong n, ulong nred)
0386 {
0387     return n_ll_redc_fast(ull_u_mul_u(a, b), n, nred);
0388 }
0389
0390 /* User-friendly functions for Montgomery arithmetic *************************/
0391
0392 typedef struct
0393 {
0394     /* Standard nmod data for odd n. The precomputed inverse is not required
0395        for Montgomery arithmetic itself, but it is useful for conversion
0396        into Montgomery form. */
0397     /* nred = -1/n mod R, R = 2^FLINT_BITS */
0398     nmod_t mod;
0399     ulong nred;
0400     /* Todo: should we store the constant 1? This would slow down
0401        creating the context when we don't need it.
0402     ulong one;
0403     */
0404 }
0405 nmod_redc_ctx_struct;
0406
0407 typedef nmod_redc_ctx_struct nmod_redc_ctx_t[1];
0408
0409 NMOD_INLINE void nmod_redc_ctx_init_nmod(nmod_redc_ctx_t ctx, nmod_t mod)
0410 {
0411     FLINT_ASSERT(mod.n & 1);
0412
0413     ctx->mod = mod;
0414     ctx->nred = n_binvert(-mod.n);
0415 }
0416
0417 NMOD_INLINE void nmod_redc_ctx_init_ui(nmod_redc_ctx_t ctx, ulong n)
0418 {
0419     FLINT_ASSERT(n & 1);
0420
0421     nmod_init(&ctx->mod, n);
0422     ctx->nred = n_binvert(-n);
0423 }
0424
0425 NMOD_INLINE ulong nmod_redc_set_nmod(ulong x, const nmod_redc_ctx_t ctx)
0426 {
0427     return n_to_redc_preinv(x, ctx->mod.n, ctx->mod.ninv, ctx->mod.norm);
0428 }
0429
0430 NMOD_INLINE ulong nmod_redc_set_ui(ulong x, const nmod_redc_ctx_t ctx)
0431 {
0432     return nmod_redc_set_nmod(nmod_set_ui(x, ctx->mod), ctx);
0433 }
0434
0435 NMOD_INLINE ulong nmod_redc_get_nmod(ulong x, const nmod_redc_ctx_t ctx)
0436 {
0437     return n_redc(x, ctx->mod.n, ctx->nred);
0438 }
0439
0440 NMOD_INLINE ulong nmod_redc_neg(ulong x, const nmod_redc_ctx_t ctx)
0441 {
0442     return nmod_neg(x, ctx->mod);
0443 }
0444
0445 NMOD_INLINE ulong nmod_redc_add(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0446 {
0447     return nmod_add(x, y, ctx->mod);
0448 }
0449
0450 NMOD_INLINE ulong nmod_redc_sub(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0451 {
0452     return nmod_sub(x, y, ctx->mod);
0453 }
0454
0455 NMOD_INLINE ulong nmod_redc_mul(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0456 {
0457     return n_mulmod_redc(x, y, ctx->mod.n, ctx->nred);
0458 }
0459
0460 NMOD_INLINE int nmod_redc_can_use_fast(const nmod_redc_ctx_t ctx)
0461 {
0462     return ctx->mod.norm >= 2;
0463 }
0464
0465 NMOD_INLINE ulong nmod_redc_fast_mul(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0466 {
0467     return n_mulmod_redc_fast(x, y, ctx->mod.n, ctx->nred);
0468 }
0469
0470 NMOD_INLINE ulong nmod_redc_fast_neg(ulong x, const nmod_redc_ctx_t ctx)
0471 {
0472     return n_negmod(x, 2 * ctx->mod.n);
0473 }
0474
0475 NMOD_INLINE ulong nmod_redc_fast_add(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0476 {
0477     return n_addmod(x, y, 2 * ctx->mod.n);
0478 }
0479
0480 NMOD_INLINE ulong nmod_redc_fast_sub(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0481 {
0482     return n_submod(x, y, 2 * ctx->mod.n);
0483 }
0484
0485 NMOD_INLINE ulong nmod_redc_fast_normalise(ulong x, const nmod_redc_ctx_t ctx)
0486 {
0487     return x < ctx->mod.n ? x : x - ctx->mod.n;
0488 }
0489
0490 NMOD_INLINE
0491 ulong nmod_redc_fast_mul_two(ulong x, const nmod_redc_ctx_t ctx)
0492 {
0493     ulong n = ctx->mod.n;
0494
0495     if (x >= n)
0496         x -= n;
0497
0498     return 2 * x;
0499 }
0500
0501 ulong _nmod_redc_pow_ui(ulong a, ulong exp, const nmod_redc_ctx_t ctx);
0502 ulong _nmod_redc_fast_pow_ui(ulong a, ulong exp, const nmod_redc_ctx_t ctx);
0503
0504 ulong _nmod_redc_2_pow_ui(ulong exp, const nmod_redc_ctx_t ctx);
0505 ulong _nmod_redc_fast_2_pow_ui(ulong exp, const nmod_redc_ctx_t ctx);
0506
0507 /* Half-limb versions of Montgomery arithmetic */
0508
0509 FLINT_FORCE_INLINE ulong
0510 n_to_redc_half_preinv(ulong a, ulong n, ulong ninv, ulong FLINT_UNUSED(norm))
0511 {
0512     return n_mod2_preinv(((ulong) a) << (FLINT_BITS / 2), n, ninv);
0513 }
0514
0515 FLINT_FORCE_INLINE
0516 ulong n_redc_half_fast(ulong x, ulong n, ulong nred)
0517 {
0518     ulong y = (x * nred) % (UWORD(1) << (FLINT_BITS / 2));
0519     ulong z = x + n * y;
0520     return z >> (FLINT_BITS / 2);
0521 }
0522
0523 FLINT_FORCE_INLINE
0524 ulong n_redc_half(ulong x, ulong n, ulong nred)
0525 {
0526     ulong y = n_redc_half_fast(x, n, nred);
0527     if (y >= n)
0528         y -= n;
0529     return y;
0530 }
0531
0532 FLINT_FORCE_INLINE
0533 ulong n_mulmod_redc_half(ulong a, ulong b, ulong n, ulong nred)
0534 {
0535     return n_redc_half(a * b, n, nred);
0536 }
0537
0538 FLINT_FORCE_INLINE
0539 ulong n_mulmod_redc_half_fast(ulong a, ulong b, ulong n, ulong nred)
0540 {
0541     return n_redc_half_fast(a * b, n, nred);
0542 }
0543
0544 NMOD_INLINE void nmod_redc_half_ctx_init_nmod(nmod_redc_ctx_t ctx, nmod_t mod)
0545 {
0546     FLINT_ASSERT(mod.n & 1);
0547     FLINT_ASSERT(mod.norm >= FLINT_BITS / 2 + 1);
0548     ctx->mod = mod;
0549     /* todo: skip one iteration */
0550     ctx->nred = n_binvert(-mod.n) % (UWORD(1) << (FLINT_BITS / 2));
0551 }
0552
0553 NMOD_INLINE void nmod_redc_half_ctx_init_ui(nmod_redc_ctx_t ctx, ulong n)
0554 {
0555     nmod_init(&ctx->mod, n);
0556     nmod_redc_half_ctx_init_nmod(ctx, ctx->mod);
0557 }
0558
0559 NMOD_INLINE ulong nmod_redc_half_set_nmod(ulong x, const nmod_redc_ctx_t ctx)
0560 {
0561     return n_to_redc_half_preinv(x, ctx->mod.n, ctx->mod.ninv, ctx->mod.norm);
0562 }
0563
0564 NMOD_INLINE ulong nmod_redc_half_set_ui(ulong x, const nmod_redc_ctx_t ctx)
0565 {
0566     return nmod_redc_half_set_nmod(nmod_set_ui(x, ctx->mod), ctx);
0567 }
0568
0569 NMOD_INLINE ulong nmod_redc_half_get_nmod(ulong x, const nmod_redc_ctx_t ctx)
0570 {
0571     return n_redc_half(x, ctx->mod.n, ctx->nred);
0572 }
0573
0574 NMOD_INLINE ulong nmod_redc_half_add(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0575 {
0576     return nmod_add(x, y, ctx->mod);
0577 }
0578
0579 NMOD_INLINE ulong nmod_redc_half_sub(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0580 {
0581     return nmod_sub(x, y, ctx->mod);
0582 }
0583
0584 NMOD_INLINE ulong nmod_redc_half_mul(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0585 {
0586     return n_mulmod_redc_half(x, y, ctx->mod.n, ctx->nred);
0587 }
0588
0589 NMOD_INLINE int nmod_redc_half_can_use_fast(const nmod_redc_ctx_t ctx)
0590 {
0591     return ctx->mod.norm >= (FLINT_BITS / 2 + 2);
0592 }
0593
0594 NMOD_INLINE ulong nmod_redc_half_fast_mul(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0595 {
0596     return n_mulmod_redc_half_fast(x, y, ctx->mod.n, ctx->nred);
0597 }
0598
0599 NMOD_INLINE ulong nmod_redc_half_fast_add(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0600 {
0601     return n_addmod(x, y, 2 * ctx->mod.n);
0602 }
0603
0604 NMOD_INLINE ulong nmod_redc_half_fast_sub(ulong x, ulong y, const nmod_redc_ctx_t ctx)
0605 {
0606     return n_submod(x, y, 2 * ctx->mod.n);
0607 }
0608
0609 /* Rings */
0610
0611 typedef struct
0612 {
0613     nmod_t nmod;
0614     ulong a;   /* when used as finite field with defining polynomial x - a */
0615     truth_t is_prime;
0616 }
0617 _gr_nmod_ctx_struct;
0618
0619 #define NMOD_CTX_REF(ring_ctx) (&((((_gr_nmod_ctx_struct *)(ring_ctx))->nmod)))
0620 #define NMOD_CTX(ring_ctx) (*NMOD_CTX_REF(ring_ctx))
0621 #define NMOD_IS_PRIME(ring_ctx) (((_gr_nmod_ctx_struct *)(ring_ctx))->is_prime)
0622 /* when used as finite field when defining polynomial x - a, allow storing the coefficient a */
0623 #define NMOD_CTX_A(ring_ctx) (&((((_gr_nmod_ctx_struct *)(ring_ctx))->a)))
0624
0625
0626 typedef struct
0627 {
0628     nmod_redc_ctx_struct ctx;
0629     ulong one;
0630     truth_t is_prime;
0631 }
0632 _gr_nmod_redc_ctx_struct;
0633
0634 #define GR_NMOD_REDC_CTX(ring_ctx) (&((((_gr_nmod_redc_ctx_struct *)(ring_ctx))->ctx)))
0635 #define GR_NMOD_REDC_MOD(ring_ctx) (((((_gr_nmod_redc_ctx_struct *)(ring_ctx))->ctx)).mod)
0636 #define GR_NMOD_REDC_N(ring_ctx) ((((((_gr_nmod_redc_ctx_struct *)(ring_ctx))->ctx)).mod).n)
0637 #define GR_NMOD_REDC_NRED(ring_ctx) (GR_NMOD_REDC_CTX(ring_ctx)->nred)
0638 #define GR_NMOD_REDC_ONE(ring_ctx) ((((_gr_nmod_redc_ctx_struct *)(ring_ctx))->one))
0639 #define GR_NMOD_REDC_IS_PRIME(ring_ctx) ((((_gr_nmod_redc_ctx_struct *)(ring_ctx))->is_prime))
0640 #define GR_NMOD_REDC_IS_FAST(ring_ctx) ((ring_ctx)->which_ring == GR_CTX_NMOD_REDC_FAST)
0641
0642 int gr_ctx_init_nmod_redc(gr_ctx_t ctx, ulong n);
0643 int gr_ctx_init_nmod_redc_fast(gr_ctx_t ctx, ulong n);
0644
0645 /* Discrete logs a la Pohlig - Hellman ***************************************/
0646
0647 typedef struct {
0648     ulong gammapow;
0649     ulong cm;
0650 } nmod_discrete_log_pohlig_hellman_table_entry_struct;
0651
0652 typedef struct {
0653     slong exp;
0654     ulong prime;
0655     ulong gamma;
0656     ulong gammainv;
0657     ulong startingbeta;
0658     ulong co;
0659     ulong startinge;
0660     ulong idem;
0661     ulong cbound;
0662     ulong dbound;
0663     nmod_discrete_log_pohlig_hellman_table_entry_struct * table; /* length cbound */
0664 } nmod_discrete_log_pohlig_hellman_entry_struct;
0665
0666 typedef struct {
0667     nmod_t mod;         /* p is mod.n */
0668     ulong alpha;    /* p.r. of p */
0669     ulong alphainv;
0670     slong num_factors;  /* factors of p - 1*/
0671     nmod_discrete_log_pohlig_hellman_entry_struct * entries;
0672 } nmod_discrete_log_pohlig_hellman_struct;
0673
0674 typedef nmod_discrete_log_pohlig_hellman_struct nmod_discrete_log_pohlig_hellman_t[1];
0675
0676 void nmod_discrete_log_pohlig_hellman_init(
0677                 nmod_discrete_log_pohlig_hellman_t L);
0678
0679 void nmod_discrete_log_pohlig_hellman_clear(
0680                 nmod_discrete_log_pohlig_hellman_t L);
0681
0682 double nmod_discrete_log_pohlig_hellman_precompute_prime(
0683                 nmod_discrete_log_pohlig_hellman_t L,
0684                 ulong p);
0685
0686 ulong nmod_discrete_log_pohlig_hellman_run(
0687                 const nmod_discrete_log_pohlig_hellman_t L,
0688                 ulong y);
0689
0690 NMOD_INLINE ulong nmod_discrete_log_pohlig_hellman_primitive_root(
0691                 const nmod_discrete_log_pohlig_hellman_t L)
0692 {
0693     return L->alpha;
0694 }
0695
0696 #ifdef __cplusplus
0697 }
0698 #endif
0699
0700 #endif