include/flint/mpn_mod.h

0001 /*
0002     Copyright (C) 2024 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 MPN_MOD_H
0013 #define MPN_MOD_H
0014
0015 #ifdef MPN_MOD_INLINES_C
0016 #define MPN_MOD_INLINE
0017 #else
0018 #define MPN_MOD_INLINE static inline
0019 #endif
0020
0021 #include "mpn_extras.h"
0022 #include "gr_types.h"
0023
0024 #ifdef __cplusplus
0025 extern "C" {
0026 #endif
0027
0028 /* Single limbs are already dealt with well by nmod, and excluding them
0029    allows avoiding various special cases. */
0030 #define MPN_MOD_MIN_LIMBS 2
0031
0032 /* This should be small enough that we can stack-allocate mpn_mods
0033    and temporary buffers a small multiple of the size.
0034    For bigger sizes, use fmpz_mod.
0035    Before resizing this, make sure that any lookup tables that go up
0036    to MPN_MOD_MAX_LIMBS (such as tuning tables) are large enough. */
0037 #define MPN_MOD_MAX_LIMBS 16
0038
0039 typedef struct
0040 {
0041     slong nlimbs;
0042     ulong d[MPN_MOD_MAX_LIMBS];
0043     ulong dinv[MPN_MOD_MAX_LIMBS];
0044     ulong dnormed[MPN_MOD_MAX_LIMBS];
0045     flint_bitcnt_t norm;
0046     truth_t is_prime;
0047 }
0048 _mpn_mod_ctx_struct;
0049
0050 #define MPN_MOD_CTX(ctx) ((_mpn_mod_ctx_struct *)(GR_CTX_DATA_AS_PTR(ctx)))
0051 #define MPN_MOD_CTX_NLIMBS(ctx) (MPN_MOD_CTX(ctx)->nlimbs)
0052 #define MPN_MOD_CTX_MODULUS(ctx) (MPN_MOD_CTX(ctx)->d)
0053 #define MPN_MOD_CTX_MODULUS_NORMED(ctx) (MPN_MOD_CTX(ctx)->dnormed)
0054 #define MPN_MOD_CTX_MODULUS_PREINV(ctx) (MPN_MOD_CTX(ctx)->dinv)
0055 #define MPN_MOD_CTX_NORM(ctx) (MPN_MOD_CTX(ctx)->norm)
0056 #define MPN_MOD_CTX_IS_PRIME(ctx) (MPN_MOD_CTX(ctx)->is_prime)
0057 #define MPN_MOD_CTX_MODULUS_BITS(ctx) ((MPN_MOD_CTX_NLIMBS(ctx) - 1) * FLINT_BITS + (FLINT_BITS - MPN_MOD_CTX_NORM(ctx)))
0058
0059 MPN_MOD_INLINE int
0060 mpn_mod_ctx_set_is_field(gr_ctx_t ctx, truth_t is_field)
0061 {
0062     MPN_MOD_CTX_IS_PRIME(ctx) = is_field;
0063     return GR_SUCCESS;
0064 }
0065
0066 /* Basic operations and arithmetic */
0067
0068 int gr_ctx_init_mpn_mod(gr_ctx_t ctx, const fmpz_t n);
0069 int _gr_ctx_init_mpn_mod(gr_ctx_t ctx, nn_srcptr n, slong nlimbs);
0070 void gr_ctx_init_mpn_mod_randtest(gr_ctx_t ctx, flint_rand_t state);
0071
0072 int mpn_mod_ctx_write(gr_stream_t out, gr_ctx_t ctx);
0073 void mpn_mod_ctx_clear(gr_ctx_t ctx);
0074
0075 MPN_MOD_INLINE truth_t
0076 mpn_mod_ctx_is_field(gr_ctx_t ctx)
0077 {
0078     return MPN_MOD_CTX_IS_PRIME(ctx);
0079 }
0080
0081 MPN_MOD_INLINE void
0082 mpn_mod_init(nn_ptr x, gr_ctx_t ctx)
0083 {
0084     flint_mpn_zero(x, MPN_MOD_CTX_NLIMBS(ctx));
0085 }
0086
0087 MPN_MOD_INLINE void
0088 mpn_mod_clear(nn_ptr FLINT_UNUSED(x), gr_ctx_t FLINT_UNUSED(ctx))
0089 {
0090 }
0091
0092 MPN_MOD_INLINE void
0093 mpn_mod_swap(nn_ptr x, nn_ptr y, gr_ctx_t ctx)
0094 {
0095     slong i = 0, n = MPN_MOD_CTX_NLIMBS(ctx);
0096     for (i = 0; i < n; i++)
0097         FLINT_SWAP(ulong, x[i], y[i]);
0098 }
0099
0100 MPN_MOD_INLINE int
0101 mpn_mod_set(nn_ptr res, nn_srcptr x, gr_ctx_t ctx)
0102 {
0103     flint_mpn_copyi(res, x, MPN_MOD_CTX_NLIMBS(ctx));
0104     return GR_SUCCESS;
0105 }
0106
0107 MPN_MOD_INLINE int
0108 mpn_mod_zero(nn_ptr res, gr_ctx_t ctx)
0109 {
0110     flint_mpn_zero(res, MPN_MOD_CTX_NLIMBS(ctx));
0111     return GR_SUCCESS;
0112 }
0113
0114 MPN_MOD_INLINE
0115 int mpn_mod_one(nn_ptr res, gr_ctx_t ctx)
0116 {
0117     res[0] = 1;
0118     flint_mpn_zero(res + 1, MPN_MOD_CTX_NLIMBS(ctx) - 1);
0119     return GR_SUCCESS;
0120 }
0121
0122 int mpn_mod_set_ui(nn_ptr res, ulong x, gr_ctx_t ctx);
0123 int mpn_mod_set_si(nn_ptr res, slong x, gr_ctx_t ctx);
0124 int mpn_mod_neg_one(nn_ptr res, gr_ctx_t ctx);
0125
0126 int mpn_mod_set_mpn(nn_ptr res, nn_srcptr x, slong xn, gr_ctx_t ctx);
0127 int mpn_mod_set_fmpz(nn_ptr res, const fmpz_t x, gr_ctx_t ctx);
0128 int mpn_mod_set_other(nn_ptr res, gr_ptr v, gr_ctx_t v_ctx, gr_ctx_t ctx);
0129 int mpn_mod_randtest(nn_ptr res, flint_rand_t state, gr_ctx_t ctx);
0130 int mpn_mod_write(gr_stream_t out, nn_srcptr x, gr_ctx_t ctx);
0131
0132 int mpn_mod_get_fmpz(fmpz_t res, nn_srcptr x, gr_ctx_t ctx);
0133
0134 MPN_MOD_INLINE truth_t
0135 mpn_mod_is_zero(nn_srcptr x, gr_ctx_t ctx)
0136 {
0137     return flint_mpn_zero_p(x, MPN_MOD_CTX_NLIMBS(ctx)) ? T_TRUE : T_FALSE;
0138 }
0139
0140 MPN_MOD_INLINE truth_t
0141 mpn_mod_is_one(nn_srcptr x, gr_ctx_t ctx)
0142 {
0143     return (x[0] == 1 && flint_mpn_zero_p(x + 1, MPN_MOD_CTX_NLIMBS(ctx) - 1)) ? T_TRUE : T_FALSE;
0144 }
0145
0146 truth_t mpn_mod_is_neg_one(gr_srcptr x, gr_ctx_t ctx);
0147
0148 MPN_MOD_INLINE truth_t
0149 mpn_mod_equal(nn_srcptr x, nn_srcptr y, gr_ctx_t ctx)
0150 {
0151     return flint_mpn_equal_p(x, y, MPN_MOD_CTX_NLIMBS(ctx)) ? T_TRUE : T_FALSE;
0152 }
0153
0154 int mpn_mod_neg(nn_ptr res, nn_srcptr x, gr_ctx_t ctx);
0155 int mpn_mod_add(nn_ptr res, nn_srcptr x, nn_srcptr y, gr_ctx_t ctx);
0156 int mpn_mod_sub(nn_ptr res, nn_srcptr x, nn_srcptr y, gr_ctx_t ctx);
0157 int mpn_mod_add_ui(nn_ptr res, nn_srcptr x, ulong y, gr_ctx_t ctx);
0158 int mpn_mod_sub_ui(nn_ptr res, nn_srcptr x, ulong y, gr_ctx_t ctx);
0159 int mpn_mod_add_si(nn_ptr res, nn_srcptr x, slong y, gr_ctx_t ctx);
0160 int mpn_mod_sub_si(nn_ptr res, nn_srcptr x, slong y, gr_ctx_t ctx);
0161 int mpn_mod_add_fmpz(nn_ptr res, nn_srcptr x, const fmpz_t y, gr_ctx_t ctx);
0162 int mpn_mod_sub_fmpz(nn_ptr res, nn_srcptr x, const fmpz_t y, gr_ctx_t ctx);
0163
0164 int mpn_mod_mul(nn_ptr res, nn_srcptr x, nn_srcptr y, gr_ctx_t ctx);
0165
0166 int mpn_mod_mul_ui(nn_ptr res, nn_srcptr x, ulong y, gr_ctx_t ctx);
0167 int mpn_mod_mul_si(nn_ptr res, nn_srcptr x, slong y, gr_ctx_t ctx);
0168 int mpn_mod_mul_fmpz(nn_ptr res, nn_srcptr x, const fmpz_t y, gr_ctx_t ctx);
0169 int mpn_mod_addmul(nn_ptr res, nn_srcptr x, nn_srcptr y, gr_ctx_t ctx);
0170 int mpn_mod_addmul_ui(nn_ptr res, nn_srcptr x, ulong y, gr_ctx_t ctx);
0171 int mpn_mod_addmul_si(nn_ptr res, nn_srcptr x, slong y, gr_ctx_t ctx);
0172 int mpn_mod_addmul_fmpz(nn_ptr res, nn_srcptr x, const fmpz_t y, gr_ctx_t ctx);
0173 int mpn_mod_submul(nn_ptr res, nn_srcptr x, nn_srcptr y, gr_ctx_t ctx);
0174 int mpn_mod_submul_ui(nn_ptr res, nn_srcptr x, ulong y, gr_ctx_t ctx);
0175 int mpn_mod_submul_si(nn_ptr res, nn_srcptr x, slong y, gr_ctx_t ctx);
0176 int mpn_mod_submul_fmpz(nn_ptr res, nn_srcptr x, const fmpz_t y, gr_ctx_t ctx);
0177
0178 MPN_MOD_INLINE int
0179 mpn_mod_sqr(nn_ptr res, nn_srcptr x, gr_ctx_t ctx)
0180 {
0181     return mpn_mod_mul(res, x, x, ctx);
0182 }
0183
0184 int mpn_mod_fmma(nn_ptr res, nn_srcptr x1, nn_srcptr y1, nn_srcptr x2, nn_srcptr y2, gr_ctx_t ctx);
0185
0186 int mpn_mod_inv(nn_ptr res, nn_srcptr x, gr_ctx_t ctx);
0187 int mpn_mod_div(nn_ptr res, nn_srcptr x, nn_srcptr y, gr_ctx_t ctx);
0188
0189 /* Vector functions */
0190
0191 int _mpn_mod_vec_zero(nn_ptr res, slong len, gr_ctx_t ctx);
0192 int _mpn_mod_vec_clear(nn_ptr FLINT_UNUSED(res), slong FLINT_UNUSED(len), gr_ctx_t FLINT_UNUSED(ctx));
0193 int _mpn_mod_vec_set(nn_ptr res, nn_srcptr x, slong len, gr_ctx_t ctx);
0194 void _mpn_mod_vec_swap(nn_ptr vec1, nn_ptr vec2, slong len, gr_ctx_t ctx);
0195 int _mpn_mod_vec_neg(nn_ptr res, nn_srcptr x, slong len, gr_ctx_t ctx);
0196 int _mpn_mod_vec_add(nn_ptr res, nn_srcptr x, nn_srcptr y, slong len, gr_ctx_t ctx);
0197 int _mpn_mod_vec_sub(nn_ptr res, nn_srcptr x, nn_srcptr y, slong len, gr_ctx_t ctx);
0198 int _mpn_mod_vec_mul(nn_ptr res, nn_srcptr x, nn_srcptr y, slong len, gr_ctx_t ctx);
0199 int _mpn_mod_vec_mul_scalar(nn_ptr res, nn_srcptr x, slong len, nn_srcptr y, gr_ctx_t ctx);
0200 int _mpn_mod_scalar_mul_vec(nn_ptr res, nn_srcptr y, nn_srcptr x, slong len, gr_ctx_t ctx);
0201 int _mpn_mod_vec_addmul_scalar(nn_ptr res, nn_srcptr x, slong len, nn_srcptr y, gr_ctx_t ctx);
0202 int _mpn_mod_vec_submul_scalar(nn_ptr res, nn_srcptr x, slong len, nn_srcptr y, gr_ctx_t ctx);
0203 int _mpn_mod_vec_dot(nn_ptr res, nn_srcptr initial, int subtract, nn_srcptr vec1, nn_srcptr vec2, slong len, gr_ctx_t ctx);
0204 int _mpn_mod_vec_dot_rev(nn_ptr res, nn_srcptr initial, int subtract, nn_srcptr vec1, nn_srcptr vec2, slong len, gr_ctx_t ctx);
0205 int _mpn_mod_vec_dot_strided(nn_ptr res, nn_srcptr initial, int subtract, nn_srcptr vec1, slong stride1, nn_srcptr vec2, slong stride2, slong len, gr_ctx_t ctx);
0206
0207 /* Matrix algorithms */
0208
0209 int mpn_mod_mat_mul_waksman(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
0210 int mpn_mod_mat_mul_multi_mod(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
0211 int mpn_mod_mat_mul(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
0212
0213 int mpn_mod_mat_nonsingular_solve_tril(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx);
0214 int mpn_mod_mat_nonsingular_solve_triu(gr_mat_t X, const gr_mat_t U, const gr_mat_t B, int unit, gr_ctx_t ctx);
0215 int mpn_mod_mat_lu_classical_delayed(slong * res_rank, slong * P, gr_mat_t A, const gr_mat_t A_in, int rank_check, gr_ctx_t ctx);
0216 int mpn_mod_mat_lu(slong * rank, slong * P, gr_mat_t LU, const gr_mat_t A, int rank_check, gr_ctx_t ctx);
0217 int mpn_mod_mat_det(nn_ptr res, const gr_mat_t A, gr_ctx_t ctx);
0218 int _mpn_mod_mat_charpoly(nn_ptr res, const gr_mat_t mat, gr_ctx_t ctx);
0219 int mpn_mod_mat_reduce_row(slong * column, gr_mat_t A, slong * P, slong * L, slong m, gr_ctx_t ctx);
0220
0221 /* Polynomial algorithms */
0222
0223 int _mpn_mod_poly_mulmid_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
0224 int _mpn_mod_poly_mulmid_karatsuba(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, slong cutoff, gr_ctx_t ctx);
0225 int _mpn_mod_poly_mulmid_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
0226 int _mpn_mod_poly_mulmid_fft_small(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
0227 int _mpn_mod_poly_mulmid(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong nlo, slong nhi, gr_ctx_t ctx);
0228
0229 int _mpn_mod_poly_mullow_classical(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx);
0230 int _mpn_mod_poly_mullow_karatsuba(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, slong cutoff, gr_ctx_t ctx);
0231 int _mpn_mod_poly_mullow_KS(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx);
0232 int _mpn_mod_poly_mullow_fft_small(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx);
0233 int _mpn_mod_poly_mullow(nn_ptr res, nn_srcptr poly1, slong len1, nn_srcptr poly2, slong len2, slong len, gr_ctx_t ctx);
0234
0235 int _mpn_mod_poly_inv_series(nn_ptr Q, nn_srcptr B, slong lenB, slong len, gr_ctx_t ctx);
0236 int _mpn_mod_poly_div_series(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, slong len, gr_ctx_t ctx);
0237
0238 int _mpn_mod_poly_divrem_q1_preinv1_fmma(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nn_srcptr invL, gr_ctx_t ctx);
0239 int _mpn_mod_poly_divrem_q1_preinv1_fmma_precond(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nn_srcptr invL, gr_ctx_t ctx);
0240 int _mpn_mod_poly_divrem_q1_preinv1_karatsuba_precond(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nn_srcptr invL, gr_ctx_t ctx);
0241 int _mpn_mod_poly_divrem_q1_preinv1(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nn_srcptr invL, gr_ctx_t ctx);
0242
0243 int _mpn_mod_poly_divrem_basecase_preinv1(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, nn_srcptr invL, gr_ctx_t ctx);
0244 int _mpn_mod_poly_divrem_basecase(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx);
0245 int _mpn_mod_poly_divrem(nn_ptr Q, nn_ptr R, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx);
0246 int _mpn_mod_poly_div(nn_ptr Q, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx);
0247
0248 int _mpn_mod_poly_gcd(nn_ptr G, slong * lenG, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx);
0249 int _mpn_mod_poly_xgcd(slong * lenG, nn_ptr G, nn_ptr S, nn_ptr T, nn_srcptr A, slong lenA, nn_srcptr B, slong lenB, gr_ctx_t ctx);
0250
0251 #ifdef __cplusplus
0252 }
0253 #endif
0254
0255 #endif