include/flint/ca_mat.h

0001 /*
0002     Copyright (C) 2020 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 CA_MAT_H
0013 #define CA_MAT_H
0014
0015 #ifdef CA_MAT_INLINES_C
0016 #define CA_MAT_INLINE
0017 #else
0018 #define CA_MAT_INLINE static inline
0019 #endif
0020
0021 #include "ca.h"
0022
0023 #ifdef __cplusplus
0024 extern "C" {
0025 #endif
0026
0027 /* Matrix object */
0028
0029 #define ca_mat_entry(mat,i,j) ((mat)->entries + (i) * (mat)->stride + (j))
0030 #define ca_mat_nrows(mat) ((mat)->r)
0031 #define ca_mat_ncols(mat) ((mat)->c)
0032
0033 CA_MAT_INLINE ca_ptr
0034 ca_mat_entry_ptr(ca_mat_t mat, slong i, slong j)
0035 {
0036     return ca_mat_entry(mat, i, j);
0037 }
0038
0039 /* Memory management */
0040
0041 void ca_mat_init(ca_mat_t mat, slong r, slong c, ca_ctx_t ctx);
0042
0043 void ca_mat_clear(ca_mat_t mat, ca_ctx_t ctx);
0044
0045 CA_MAT_INLINE void
0046 ca_mat_swap(ca_mat_t mat1, ca_mat_t mat2, ca_ctx_t FLINT_UNUSED(ctx))
0047 {
0048     FLINT_SWAP(ca_mat_struct, *mat1, *mat2);
0049 }
0050
0051 /* Window matrices */
0052
0053 void ca_mat_window_init(ca_mat_t window, const ca_mat_t mat, slong r1, slong c1, slong r2, slong c2, ca_ctx_t ctx);
0054
0055 CA_MAT_INLINE void
0056 ca_mat_window_clear(ca_mat_t FLINT_UNUSED(window), ca_ctx_t FLINT_UNUSED(ctx))
0057 {
0058 }
0059
0060 /* Shape */
0061
0062 CA_MAT_INLINE int
0063 ca_mat_is_empty(const ca_mat_t mat)
0064 {
0065     return (mat->r == 0) || (mat->c == 0);
0066 }
0067
0068 CA_MAT_INLINE int
0069 ca_mat_is_square(const ca_mat_t mat)
0070 {
0071     return (mat->r == mat->c);
0072 }
0073
0074 /* Conversions */
0075
0076 void ca_mat_set(ca_mat_t dest, const ca_mat_t src, ca_ctx_t ctx);
0077 void ca_mat_set_fmpz_mat(ca_mat_t dest, const fmpz_mat_t src, ca_ctx_t ctx);
0078 void ca_mat_set_fmpq_mat(ca_mat_t dest, const fmpq_mat_t src, ca_ctx_t ctx);
0079 void ca_mat_set_ca(ca_mat_t y, const ca_t x, ca_ctx_t ctx);
0080 void ca_mat_get_fexpr(fexpr_t res, const ca_mat_t A, ulong flags, ca_ctx_t ctx);
0081
0082 void ca_mat_transfer(ca_mat_t res, ca_ctx_t res_ctx, const ca_mat_t src, ca_ctx_t src_ctx);
0083
0084 /* Random generation */
0085
0086 void ca_mat_randtest(ca_mat_t mat, flint_rand_t state, slong len, slong bits, ca_ctx_t ctx);
0087 void ca_mat_randtest_rational(ca_mat_t mat, flint_rand_t state, slong bits, ca_ctx_t ctx);
0088 void ca_mat_randops(ca_mat_t mat, flint_rand_t state, slong count, ca_ctx_t ctx);
0089
0090 /* I/O */
0091
0092 void ca_mat_print(const ca_mat_t mat, ca_ctx_t ctx);
0093 void ca_mat_printn(const ca_mat_t mat, slong digits, ca_ctx_t ctx);
0094
0095 /* Special matrices */
0096
0097 void ca_mat_zero(ca_mat_t mat, ca_ctx_t ctx);
0098 void ca_mat_one(ca_mat_t mat, ca_ctx_t ctx);
0099 void ca_mat_ones(ca_mat_t mat, ca_ctx_t ctx);
0100 void ca_mat_pascal(ca_mat_t mat, int triangular, ca_ctx_t ctx);
0101 void ca_mat_stirling(ca_mat_t mat, int kind, ca_ctx_t ctx);
0102 void ca_mat_hilbert(ca_mat_t mat, ca_ctx_t ctx);
0103 void ca_mat_dft(ca_mat_t res, int type, ca_ctx_t ctx);
0104
0105 /* Comparisons and properties */
0106
0107 truth_t ca_mat_check_equal(const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0108 truth_t ca_mat_check_is_zero(const ca_mat_t A, ca_ctx_t ctx);
0109 truth_t ca_mat_check_is_one(const ca_mat_t A, ca_ctx_t ctx);
0110
0111 /* Conjugate and transpose */
0112
0113 void ca_mat_transpose(ca_mat_t B, const ca_mat_t A, ca_ctx_t ctx);
0114 void ca_mat_conj(ca_mat_t B, const ca_mat_t A, ca_ctx_t ctx);
0115 void ca_mat_conj_transpose(ca_mat_t mat1, const ca_mat_t mat2, ca_ctx_t ctx);
0116
0117 /* Arithmetic */
0118
0119 void ca_mat_add_ca(ca_mat_t y, const ca_mat_t a, const ca_t x, ca_ctx_t ctx);
0120 void ca_mat_sub_ca(ca_mat_t y, const ca_mat_t a, const ca_t x, ca_ctx_t ctx);
0121 void ca_mat_addmul_ca(ca_mat_t y, const ca_mat_t a, const ca_t x, ca_ctx_t ctx);
0122 void ca_mat_submul_ca(ca_mat_t y, const ca_mat_t a, const ca_t x, ca_ctx_t ctx);
0123
0124 void ca_mat_neg(ca_mat_t dest, const ca_mat_t src, ca_ctx_t ctx);
0125 void ca_mat_add(ca_mat_t res, const ca_mat_t mat1, const ca_mat_t mat2, ca_ctx_t ctx);
0126 void ca_mat_sub(ca_mat_t res, const ca_mat_t mat1, const ca_mat_t mat2, ca_ctx_t ctx);
0127
0128 void ca_mat_mul(ca_mat_t C, const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0129
0130 void ca_mat_mul_classical(ca_mat_t C, const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0131 void ca_mat_mul_same_nf(ca_mat_t C, const ca_mat_t A, const ca_mat_t B, ca_field_t K, ca_ctx_t ctx);
0132
0133 CA_MAT_INLINE void
0134 ca_mat_mul_si(ca_mat_t B, const ca_mat_t A, slong c, ca_ctx_t ctx)
0135 {
0136     slong i, j;
0137
0138     for (i = 0; i < ca_mat_nrows(A); i++)
0139         for (j = 0; j < ca_mat_ncols(A); j++)
0140             ca_mul_si(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0141 }
0142
0143 CA_MAT_INLINE void
0144 ca_mat_mul_fmpz(ca_mat_t B, const ca_mat_t A, const fmpz_t c, ca_ctx_t ctx)
0145 {
0146     slong i, j;
0147
0148     for (i = 0; i < ca_mat_nrows(A); i++)
0149         for (j = 0; j < ca_mat_ncols(A); j++)
0150             ca_mul_fmpz(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0151 }
0152
0153 CA_MAT_INLINE void
0154 ca_mat_mul_fmpq(ca_mat_t B, const ca_mat_t A, const fmpq_t c, ca_ctx_t ctx)
0155 {
0156     slong i, j;
0157
0158     for (i = 0; i < ca_mat_nrows(A); i++)
0159         for (j = 0; j < ca_mat_ncols(A); j++)
0160             ca_mul_fmpq(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0161 }
0162
0163 CA_MAT_INLINE void
0164 ca_mat_mul_ca(ca_mat_t B, const ca_mat_t A, const ca_t c, ca_ctx_t ctx)
0165 {
0166     slong i, j;
0167
0168     for (i = 0; i < ca_mat_nrows(A); i++)
0169         for (j = 0; j < ca_mat_ncols(A); j++)
0170             ca_mul(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0171 }
0172
0173 CA_MAT_INLINE void
0174 ca_mat_div_si(ca_mat_t B, const ca_mat_t A, slong c, ca_ctx_t ctx)
0175 {
0176     slong i, j;
0177
0178     for (i = 0; i < ca_mat_nrows(A); i++)
0179         for (j = 0; j < ca_mat_ncols(A); j++)
0180             ca_div_si(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0181 }
0182
0183 CA_MAT_INLINE void
0184 ca_mat_div_fmpz(ca_mat_t B, const ca_mat_t A, const fmpz_t c, ca_ctx_t ctx)
0185 {
0186     slong i, j;
0187
0188     for (i = 0; i < ca_mat_nrows(A); i++)
0189         for (j = 0; j < ca_mat_ncols(A); j++)
0190             ca_div_fmpz(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0191 }
0192
0193 CA_MAT_INLINE void
0194 ca_mat_div_fmpq(ca_mat_t B, const ca_mat_t A, const fmpq_t c, ca_ctx_t ctx)
0195 {
0196     slong i, j;
0197
0198     for (i = 0; i < ca_mat_nrows(A); i++)
0199         for (j = 0; j < ca_mat_ncols(A); j++)
0200             ca_div_fmpq(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0201 }
0202
0203 CA_MAT_INLINE void
0204 ca_mat_div_ca(ca_mat_t B, const ca_mat_t A, const ca_t c, ca_ctx_t ctx)
0205 {
0206     slong i, j;
0207
0208     for (i = 0; i < ca_mat_nrows(A); i++)
0209         for (j = 0; j < ca_mat_ncols(A); j++)
0210             ca_div(ca_mat_entry(B, i, j), ca_mat_entry(A, i, j), c, ctx);
0211 }
0212
0213 CA_MAT_INLINE void
0214 ca_mat_sqr(ca_mat_t res, const ca_mat_t A, ca_ctx_t ctx)
0215 {
0216     ca_mat_mul(res, A, A, ctx);
0217 }
0218
0219 void ca_mat_pow_ui_binexp(ca_mat_t B, const ca_mat_t A, ulong exp, ca_ctx_t ctx);
0220
0221 /* Polynomial evaluation */
0222
0223 void _ca_mat_ca_poly_evaluate(ca_mat_t y, ca_srcptr poly, slong len, const ca_mat_t x, ca_ctx_t ctx);
0224 void ca_mat_ca_poly_evaluate(ca_mat_t res, const ca_poly_t f, const ca_mat_t a, ca_ctx_t ctx);
0225
0226 /* Trace */
0227
0228 void ca_mat_trace(ca_t trace, const ca_mat_t mat, ca_ctx_t ctx);
0229
0230 /* Gaussian elimination, solving and inverse */
0231
0232 truth_t ca_mat_find_pivot(slong * pivot_row, ca_mat_t mat, slong start_row, slong end_row, slong column, ca_ctx_t ctx);
0233
0234 CA_MAT_INLINE void
0235 _ca_mat_swap_rows(ca_mat_t mat, slong * perm, slong r, slong s)
0236 {
0237     if (r != s)
0238     {
0239         ca_ptr u, v;
0240         slong i;
0241
0242         if (perm != NULL)
0243             FLINT_SWAP(slong, perm[r], perm[s]);
0244
0245         u = ca_mat_entry(mat, r, 0);
0246         v = ca_mat_entry(mat, s, 0);
0247
0248         for (i = 0; i < mat->c; i++)
0249             FLINT_SWAP(ca_struct, u[i], v[i]);
0250
0251     }
0252 }
0253
0254 int ca_mat_lu_classical(slong * rank, slong * P, ca_mat_t LU, const ca_mat_t A, int rank_check, ca_ctx_t ctx);
0255 int ca_mat_lu_recursive(slong * rank, slong * P, ca_mat_t LU, const ca_mat_t A, int rank_check, ca_ctx_t ctx);
0256 int ca_mat_lu(slong * rank, slong * P, ca_mat_t LU, const ca_mat_t A, int rank_check, ca_ctx_t ctx);
0257
0258 int ca_mat_fflu(slong * rank, slong * P, ca_mat_t LU, ca_t den, const ca_mat_t A, int rank_check, ca_ctx_t ctx);
0259
0260 int ca_mat_rref_fflu(slong * rank, ca_mat_t R, const ca_mat_t A, ca_ctx_t ctx);
0261 int ca_mat_rref_lu(slong * rank, ca_mat_t R, const ca_mat_t A, ca_ctx_t ctx);
0262 int ca_mat_rref(slong * rank, ca_mat_t R, const ca_mat_t A, ca_ctx_t ctx);
0263
0264 truth_t ca_mat_nonsingular_lu(slong * P, ca_mat_t LU, const ca_mat_t A, ca_ctx_t ctx);
0265 truth_t ca_mat_nonsingular_fflu(slong * P, ca_mat_t LU, ca_t den, const ca_mat_t A, ca_ctx_t ctx);
0266
0267 truth_t ca_mat_nonsingular_solve_adjugate(ca_mat_t X, const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0268 truth_t ca_mat_nonsingular_solve_fflu(ca_mat_t X, const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0269 truth_t ca_mat_nonsingular_solve_lu(ca_mat_t X, const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0270 truth_t ca_mat_nonsingular_solve(ca_mat_t X, const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0271
0272 truth_t ca_mat_inv(ca_mat_t X, const ca_mat_t A, ca_ctx_t ctx);
0273
0274 void ca_mat_solve_tril_classical(ca_mat_t X, const ca_mat_t L, const ca_mat_t B, int unit, ca_ctx_t ctx);
0275 void ca_mat_solve_tril_recursive(ca_mat_t X, const ca_mat_t L, const ca_mat_t B, int unit, ca_ctx_t ctx);
0276 void ca_mat_solve_tril(ca_mat_t X, const ca_mat_t L, const ca_mat_t B, int unit, ca_ctx_t ctx);
0277
0278 void ca_mat_solve_triu_classical(ca_mat_t X, const ca_mat_t U, const ca_mat_t B, int unit, ca_ctx_t ctx);
0279 void ca_mat_solve_triu_recursive(ca_mat_t X, const ca_mat_t U, const ca_mat_t B, int unit, ca_ctx_t ctx);
0280 void ca_mat_solve_triu(ca_mat_t X, const ca_mat_t U, const ca_mat_t B, int unit, ca_ctx_t ctx);
0281
0282 void ca_mat_solve_lu_precomp(ca_mat_t X, const slong * perm, const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0283 void ca_mat_solve_fflu_precomp(ca_mat_t X, const slong * perm, const ca_mat_t A, const ca_t den, const ca_mat_t B, ca_ctx_t ctx);
0284
0285 /* Rank and kernel */
0286
0287 int ca_mat_rank(slong * rank, const ca_mat_t A, ca_ctx_t ctx);
0288
0289 int ca_mat_right_kernel(ca_mat_t X, const ca_mat_t A, ca_ctx_t ctx);
0290
0291 /* Determinant */
0292
0293 void ca_mat_det_berkowitz(ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0294 int ca_mat_det_lu(ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0295 int ca_mat_det_bareiss(ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0296 void ca_mat_det_cofactor(ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0297 void ca_mat_det(ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0298
0299 void ca_mat_adjugate_cofactor(ca_mat_t adj, ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0300 void ca_mat_adjugate_charpoly(ca_mat_t adj, ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0301 void ca_mat_adjugate(ca_mat_t adj, ca_t det, const ca_mat_t A, ca_ctx_t ctx);
0302
0303 /* Characteristic polynomial */
0304
0305 void _ca_mat_charpoly_berkowitz(ca_ptr cp, const ca_mat_t mat, ca_ctx_t ctx);
0306 void ca_mat_charpoly_berkowitz(ca_poly_t cp, const ca_mat_t mat, ca_ctx_t ctx);
0307
0308 int _ca_mat_charpoly_danilevsky(ca_ptr p, const ca_mat_t A, ca_ctx_t ctx);
0309 int ca_mat_charpoly_danilevsky(ca_poly_t cp, const ca_mat_t mat, ca_ctx_t ctx);
0310
0311 void _ca_mat_charpoly(ca_ptr cp, const ca_mat_t mat, ca_ctx_t ctx);
0312 void ca_mat_charpoly(ca_poly_t cp, const ca_mat_t mat, ca_ctx_t ctx);
0313
0314
0315
0316 int ca_mat_companion(ca_mat_t A, const ca_poly_t poly, ca_ctx_t ctx);
0317
0318 /* Eigenvalues and eigenvectors */
0319
0320 int ca_mat_eigenvalues(ca_vec_t lambda, ulong * exp, const ca_mat_t mat, ca_ctx_t ctx);
0321
0322 /* Diagonalization */
0323
0324 truth_t ca_mat_diagonalization(ca_mat_t D, ca_mat_t P, const ca_mat_t A, ca_ctx_t ctx);
0325
0326 void ca_mat_set_jordan_blocks(ca_mat_t mat, const ca_vec_t lambda, slong num_blocks, slong * block_lambda, slong * block_size, ca_ctx_t ctx);
0327 int ca_mat_jordan_blocks(ca_vec_t lambda, slong * num_blocks, slong * block_lambda, slong * block_size, const ca_mat_t A, ca_ctx_t ctx);
0328 int ca_mat_jordan_transformation(ca_mat_t mat, const ca_vec_t lambda, slong num_blocks, slong * block_lambda, slong * block_size, const ca_mat_t A, ca_ctx_t ctx);
0329 int ca_mat_jordan_form(ca_mat_t J, ca_mat_t P, const ca_mat_t A, ca_ctx_t ctx);
0330
0331 /* Matrix functions */
0332
0333 int ca_mat_exp(ca_mat_t res, const ca_mat_t A, ca_ctx_t ctx);
0334 truth_t ca_mat_log(ca_mat_t res, const ca_mat_t A, ca_ctx_t ctx);
0335
0336 /* Internal representation */
0337
0338 /* todo: document, make consistent */
0339 ca_field_ptr _ca_mat_same_field(const ca_mat_t A, ca_ctx_t ctx);
0340 ca_field_ptr _ca_mat_same_field2(const ca_mat_t A, const ca_mat_t B, ca_ctx_t ctx);
0341
0342 #ifdef __cplusplus
0343 }
0344 #endif
0345
0346 #endif