include/flint/nmod_mat.h

0001 /*
0002     Copyright (C) 2010 William Hart
0003     Copyright (C) 2010,2011 Fredrik Johansson
0004     Copyright (C) 2014 Ashish Kedia
0005
0006     This file is part of FLINT.
0007
0008     FLINT is free software: you can redistribute it and/or modify it under
0009     the terms of the GNU Lesser General Public License (LGPL) as published
0010     by the Free Software Foundation; either version 3 of the License, or
0011     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
0012 */
0013
0014 #ifndef NMOD_MAT_H
0015 #define NMOD_MAT_H
0016
0017 #ifdef NMOD_MAT_INLINES_C
0018 #define NMOD_MAT_INLINE
0019 #else
0020 #define NMOD_MAT_INLINE static inline
0021 #endif
0022
0023 #include "nmod_types.h"
0024
0025 #ifdef __cplusplus
0026 extern "C" {
0027 #endif
0028
0029 #define nmod_mat_entry(mat,i,j) ((mat)->entries[(i) * (mat)->stride + (j)])
0030
0031 NMOD_MAT_INLINE
0032 ulong nmod_mat_get_entry(const nmod_mat_t mat, slong i, slong j)
0033 {
0034    return nmod_mat_entry(mat, i, j);
0035 }
0036
0037 NMOD_MAT_INLINE
0038 ulong * nmod_mat_entry_ptr(const nmod_mat_t mat, slong i, slong j)
0039 {
0040    return &nmod_mat_entry(mat, i, j);
0041 }
0042
0043 NMOD_MAT_INLINE
0044 ulong * nmod_mat_row_ptr(const nmod_mat_t mat, slong i)
0045 {
0046    return &nmod_mat_entry(mat, i, 0);
0047 }
0048 /* See inlines.c */
0049
0050 NMOD_MAT_INLINE
0051 slong nmod_mat_nrows(const nmod_mat_t mat)
0052 {
0053    return mat->r;
0054 }
0055
0056 NMOD_MAT_INLINE
0057 slong nmod_mat_ncols(const nmod_mat_t mat)
0058 {
0059    return mat->c;
0060 }
0061
0062 void nmod_mat_set_mod(nmod_mat_t mat, ulong n);
0063
0064 NMOD_MAT_INLINE
0065 nmod_t nmod_mat_mod(const nmod_mat_t mat)
0066 {
0067    return mat->mod;
0068 }
0069
0070 /* Memory management */
0071 void nmod_mat_init(nmod_mat_t mat, slong rows, slong cols, ulong n);
0072 void nmod_mat_init_set(nmod_mat_t mat, const nmod_mat_t src);
0073 void nmod_mat_clear(nmod_mat_t mat);
0074 void nmod_mat_one(nmod_mat_t mat);
0075
0076 void nmod_mat_swap(nmod_mat_t mat1, nmod_mat_t mat2);
0077
0078 NMOD_MAT_INLINE void
0079 nmod_mat_swap_entrywise(nmod_mat_t mat1, nmod_mat_t mat2)
0080 {
0081     slong i, j;
0082     for (i = 0; i < nmod_mat_nrows(mat1); i++)
0083     {
0084        ulong * row1 = nmod_mat_entry_ptr(mat1, i, 0);
0085        ulong * row2 = nmod_mat_entry_ptr(mat2, i, 0);
0086        for (j = 0; j < nmod_mat_ncols(mat1); j++)
0087           FLINT_SWAP(ulong, row1[j], row2[j]);
0088     }
0089 }
0090
0091 /* Windows and concatenation */
0092
0093 NMOD_MAT_INLINE void
0094 nmod_mat_window_init(nmod_mat_t window, const nmod_mat_t mat,
0095     slong r1, slong c1, slong r2, slong c2)
0096 {
0097     FLINT_ASSERT(r1 >= 0 && r1 <= r2 && r2 <= mat->r);
0098     FLINT_ASSERT(c2 >= 0 && c1 <= c2 && c2 <= mat->c);
0099
0100     window->entries = nmod_mat_entry_ptr(mat, r1, c1);
0101     window->r = r2 - r1;
0102     window->c = c2 - c1;
0103     window->stride = mat->stride;
0104     window->mod = mat->mod;
0105 }
0106
0107 NMOD_MAT_INLINE void
0108 nmod_mat_window_clear(nmod_mat_t FLINT_UNUSED(window))
0109 {
0110 }
0111
0112 void nmod_mat_concat_horizontal(nmod_mat_t res,
0113                            const nmod_mat_t mat1,  const nmod_mat_t mat2);
0114 void nmod_mat_concat_vertical(nmod_mat_t res,
0115                            const nmod_mat_t mat1,  const nmod_mat_t mat2);
0116
0117 /* Random matrix generation */
0118 void nmod_mat_randtest(nmod_mat_t mat, flint_rand_t state);
0119 void nmod_mat_randfull(nmod_mat_t mat, flint_rand_t state);
0120 void nmod_mat_rand(nmod_mat_t mat, flint_rand_t state);
0121 int nmod_mat_randpermdiag(nmod_mat_t mat, flint_rand_t state,
0122                  nn_srcptr diag, slong n);
0123 void nmod_mat_randrank(nmod_mat_t mat, flint_rand_t state, slong rank);
0124 void nmod_mat_randops(nmod_mat_t mat, flint_rand_t state, slong count);
0125 void nmod_mat_randtril(nmod_mat_t mat, flint_rand_t state, int unit);
0126 void nmod_mat_randtriu(nmod_mat_t mat, flint_rand_t state, int unit);
0127
0128 #ifdef FLINT_HAVE_FILE
0129 int nmod_mat_fprint_pretty(FILE* file, const nmod_mat_t mat);
0130 int nmod_mat_fprint(FILE* f, const nmod_mat_t mat);
0131 #endif
0132
0133 void nmod_mat_print_pretty(const nmod_mat_t mat);
0134 int nmod_mat_print(const nmod_mat_t mat);
0135
0136 int nmod_mat_equal(const nmod_mat_t mat1, const nmod_mat_t mat2);
0137
0138 void nmod_mat_zero(nmod_mat_t mat);
0139
0140 int nmod_mat_is_zero(const nmod_mat_t mat);
0141 int nmod_mat_is_one(const nmod_mat_t mat);
0142 int nmod_mat_is_zero_row(const nmod_mat_t mat, slong i);
0143
0144 NMOD_MAT_INLINE
0145 int nmod_mat_is_empty(const nmod_mat_t mat)
0146 {
0147     return (mat->r == 0) || (mat->c == 0);
0148 }
0149
0150 NMOD_MAT_INLINE
0151 int nmod_mat_is_square(const nmod_mat_t mat)
0152 {
0153     return (mat->r == mat->c);
0154 }
0155
0156
0157 void nmod_mat_set(nmod_mat_t B, const nmod_mat_t A);
0158 void nmod_mat_transpose(nmod_mat_t B, const nmod_mat_t A);
0159
0160 /* Addition and subtraction */
0161
0162 void nmod_mat_add(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B);
0163 void nmod_mat_sub(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B);
0164 void nmod_mat_neg(nmod_mat_t B, const nmod_mat_t A);
0165
0166 /* Matrix-scalar arithmetic */
0167
0168 void _nmod_mat_scalar_mul_generic(nmod_mat_t B, const nmod_mat_t A, ulong c);
0169 void _nmod_mat_scalar_mul_precomp(nmod_mat_t B, const nmod_mat_t A, ulong c, ulong c_pr);
0170 void nmod_mat_scalar_mul(nmod_mat_t B, const nmod_mat_t A, ulong c);
0171 void nmod_mat_scalar_mul_fmpz(nmod_mat_t B, const nmod_mat_t A, const fmpz_t c);
0172 void _nmod_mat_scalar_addmul_ui_generic(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B, const ulong c);
0173 void _nmod_mat_scalar_addmul_ui_precomp(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B, const ulong c, const ulong c_pr);
0174 void nmod_mat_scalar_addmul_ui(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B, const ulong c);
0175
0176 /* Matrix multiplication */
0177
0178 void nmod_mat_mul(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B);
0179
0180 int nmod_mat_mul_blas(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B);
0181
0182 void nmod_mat_mul_classical(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B);
0183
0184 void
0185 _nmod_mat_mul_classical_threaded_pool_op(nmod_mat_t D, const nmod_mat_t C,
0186                     const nmod_mat_t A, const nmod_mat_t B, int op,
0187                   thread_pool_handle * threads, slong num_threads);
0188
0189 void nmod_mat_mul_classical_threaded(nmod_mat_t C,
0190                                const nmod_mat_t A, const nmod_mat_t B);
0191 void nmod_mat_mul_strassen(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B);
0192
0193 void _nmod_mat_mul_classical_op(nmod_mat_t D, const nmod_mat_t C,
0194                                 const nmod_mat_t A, const nmod_mat_t B, int op);
0195
0196 void nmod_mat_addmul(nmod_mat_t D, const nmod_mat_t C,
0197                                 const nmod_mat_t A, const nmod_mat_t B);
0198
0199 void nmod_mat_submul(nmod_mat_t D, const nmod_mat_t C,
0200                                 const nmod_mat_t A, const nmod_mat_t B);
0201
0202 void nmod_mat_mul_nmod_vec(ulong * c, const nmod_mat_t A,
0203                                               const ulong * b, slong blen);
0204
0205 void nmod_mat_mul_nmod_vec_ptr(ulong * const * c,
0206                   const nmod_mat_t A, const ulong * const * b, slong blen);
0207
0208 void nmod_mat_nmod_vec_mul(ulong * c, const ulong * a,
0209                                                slong alen, const nmod_mat_t B);
0210
0211 void nmod_mat_nmod_vec_mul_ptr(ulong * const * c,
0212                   const ulong * const * a, slong alen, const nmod_mat_t B);
0213
0214 /* Exponent */
0215
0216 void _nmod_mat_pow(nmod_mat_t dest, const nmod_mat_t mat, ulong pow);
0217 void nmod_mat_pow(nmod_mat_t dest, const nmod_mat_t mat, ulong pow);
0218
0219 /* Trace */
0220
0221 ulong nmod_mat_trace(const nmod_mat_t mat);
0222
0223 /* Determinant */
0224
0225 ulong _nmod_mat_det(nmod_mat_t A);
0226 ulong nmod_mat_det(const nmod_mat_t A);
0227
0228 ulong _nmod_mat_det_howell(nmod_mat_t A);
0229 ulong nmod_mat_det_howell(const nmod_mat_t A);
0230
0231 /* Rank */
0232
0233 slong nmod_mat_rank(const nmod_mat_t A);
0234
0235 /* Inverse */
0236
0237 int nmod_mat_inv(nmod_mat_t B, const nmod_mat_t A);
0238
0239 /* Permutations */
0240
0241 NMOD_MAT_INLINE
0242 void nmod_mat_swap_rows(nmod_mat_t mat, slong * perm, slong r, slong s)
0243 {
0244     if (r != s && !nmod_mat_is_empty(mat))
0245     {
0246         slong i;
0247         nn_ptr u, v;
0248
0249         if (perm)
0250             FLINT_SWAP(slong, perm[r], perm[s]);
0251
0252         u = nmod_mat_entry_ptr(mat, r, 0);
0253         v = nmod_mat_entry_ptr(mat, s, 0);
0254
0255         for (i = 0; i < mat->c; i++)
0256             FLINT_SWAP(ulong, u[i], v[i]);
0257     }
0258 }
0259
0260 NMOD_MAT_INLINE
0261 void nmod_mat_invert_rows(nmod_mat_t mat, slong * perm)
0262 {
0263     slong i;
0264
0265     for (i = 0; i < mat->r/2; i++)
0266         nmod_mat_swap_rows(mat, perm, i, mat->r - i - 1);
0267 }
0268
0269 NMOD_MAT_INLINE
0270 void nmod_mat_swap_cols(nmod_mat_t mat, slong * perm, slong r, slong s)
0271 {
0272     if (r != s && !nmod_mat_is_empty(mat))
0273     {
0274         slong i;
0275
0276         if (perm != NULL)
0277             FLINT_SWAP(slong, perm[r], perm[s]);
0278
0279         for (i = 0; i < mat->r; i++)
0280             FLINT_SWAP(ulong, nmod_mat_entry(mat, i, r), nmod_mat_entry(mat, i, s));
0281     }
0282 }
0283
0284 NMOD_MAT_INLINE
0285 void nmod_mat_invert_cols(nmod_mat_t mat, slong * perm)
0286 {
0287     if (!nmod_mat_is_empty(mat))
0288     {
0289         slong t, i;
0290         slong c = mat->c;
0291         slong k = mat->c/2;
0292
0293         if (perm != NULL)
0294             for (i = 0; i < k; i++)
0295                 FLINT_SWAP(slong, perm[i], perm[c - i - 1]);
0296
0297         for (t = 0; t < mat->r; t++)
0298             for (i = 0; i < k; i++)
0299                 FLINT_SWAP(ulong, nmod_mat_entry(mat, t, i), nmod_mat_entry(mat, t, c - i - 1));
0300     }
0301 }
0302
0303 void nmod_mat_permute_rows(nmod_mat_t mat, const slong * perm_act, slong * perm_store);
0304
0305 /* Triangular solving */
0306
0307 void nmod_mat_solve_tril(nmod_mat_t X, const nmod_mat_t L, const nmod_mat_t B, int unit);
0308 void nmod_mat_solve_tril_recursive(nmod_mat_t X, const nmod_mat_t L, const nmod_mat_t B, int unit);
0309 void nmod_mat_solve_tril_classical(nmod_mat_t X, const nmod_mat_t L, const nmod_mat_t B, int unit);
0310
0311 void nmod_mat_solve_triu(nmod_mat_t X, const nmod_mat_t U, const nmod_mat_t B, int unit);
0312 void nmod_mat_solve_triu_recursive(nmod_mat_t X, const nmod_mat_t U, const nmod_mat_t B, int unit);
0313 void nmod_mat_solve_triu_classical(nmod_mat_t X, const nmod_mat_t U, const nmod_mat_t B, int unit);
0314
0315 /* LU decomposition */
0316
0317 slong nmod_mat_lu(slong * P, nmod_mat_t A, int rank_check);
0318 slong nmod_mat_lu_classical(slong * P, nmod_mat_t A, int rank_check);
0319 slong nmod_mat_lu_classical_delayed(slong * P, nmod_mat_t A, int rank_check);
0320 slong nmod_mat_lu_recursive(slong * P, nmod_mat_t A, int rank_check);
0321
0322 slong nmod_mat_lu_with_pivots(slong * P, slong * pivots_nonpivots, nmod_mat_t A);
0323
0324 /* Nonsingular solving */
0325
0326 int nmod_mat_solve(nmod_mat_t X, const nmod_mat_t A, const nmod_mat_t B);
0327 int nmod_mat_solve_vec(nn_ptr x, const nmod_mat_t A, nn_srcptr b);
0328
0329 /* Solving */
0330
0331 int nmod_mat_can_solve_inner(slong * rank, slong * prm, slong * piv,
0332                           nmod_mat_t X, const nmod_mat_t A, const nmod_mat_t B);
0333
0334 int nmod_mat_can_solve(nmod_mat_t X, const nmod_mat_t A,
0335                                                             const nmod_mat_t B);
0336
0337 /* Reduced row echelon form */
0338
0339 slong nmod_mat_rref(nmod_mat_t A);
0340 slong _nmod_mat_rref(nmod_mat_t A, slong * pivots_nonpivots, slong * P);
0341 slong nmod_mat_rref_classical(nmod_mat_t A);
0342 slong _nmod_mat_rref_classical(nmod_mat_t A, slong * pivots_nonpivots);
0343 slong nmod_mat_rref_storjohann(nmod_mat_t A);
0344 slong _nmod_mat_rref_storjohann(nmod_mat_t A, slong * pivots_nonpivots);
0345
0346 slong nmod_mat_reduce_row(nmod_mat_t M, slong * P, slong * L, slong m);
0347
0348 /* Nullspace */
0349
0350 slong nmod_mat_nullspace(nmod_mat_t X, const nmod_mat_t A);
0351 slong nmod_mat_left_nullspace(nmod_mat_t X, const nmod_mat_t A);
0352
0353 /* Howell form */
0354
0355 void nmod_mat_strong_echelon_form(nmod_mat_t A);
0356
0357 slong nmod_mat_howell_form(nmod_mat_t A);
0358
0359 /* Transforms */
0360
0361 void nmod_mat_similarity(nmod_mat_t M, slong r, ulong d);
0362
0363 /* Characteristic polynomial and minimal polynomial */
0364
0365 /* The following prototype actually lives in nmod_poly.h
0366  *
0367  * int nmod_mat_charpoly_danilevsky(nmod_poly_t p, const nmod_mat_t M);
0368  *
0369  * void nmod_mat_minpoly(nmod_poly_t p, const nmod_mat_t M);
0370 */
0371
0372 /* Tuning parameters *********************************************************/
0373
0374 /* Size at which pre-transposing becomes faster in classical multiplication */
0375 #define NMOD_MAT_MUL_TRANSPOSE_CUTOFF 20
0376
0377 /* Cutoff between classical and recursive triangular solving */
0378 #define NMOD_MAT_SOLVE_TRI_ROWS_CUTOFF 64
0379 #define NMOD_MAT_SOLVE_TRI_COLS_CUTOFF 64
0380
0381 /*
0382    Suggested initial modulus size for multimodular algorithms. This should
0383    be chosen so that we get the most number of bits per cycle
0384    in matrix multiplication. On x86-64 it appears to be optimal to use
0385    moduli giving nlimbs = 2. This should hold both in the classical
0386    range and in Strassen blocks.
0387  */
0388 #define NMOD_MAT_OPTIMAL_MODULUS_BITS (FLINT_BITS-5)
0389
0390 /* Inlines *******************************************************************/
0391
0392 void nmod_mat_set_entry(nmod_mat_t mat, slong i, slong j, ulong x);
0393
0394 #ifdef __cplusplus
0395 }
0396 #endif
0397
0398 #endif