diff --git a/doc/source/fmpcb_mat.rst b/doc/source/fmpcb_mat.rst index 49a64baa..4e9ceb3a 100644 --- a/doc/source/fmpcb_mat.rst +++ b/doc/source/fmpcb_mat.rst @@ -243,3 +243,23 @@ Gaussian elimination and solving determinant of the remaining submatrix is bounded using Hadamard's inequality. + +Special functions +------------------------------------------------------------------------------- + +.. function:: void fmpcb_mat_exp(fmpcb_mat_t B, const fmpcb_mat_t A, long prec) + + Sets *B* to the exponential of the matrix *A*, defined by the Taylor series + + .. math :: + + \exp(A) = \sum_{k=0}^{\infty} \frac{A^k}{k!}. + + The exponential function is evaluated using scaling followed by + rectangular splitting evaluation of the Taylor series. + Scaling amounts to picking a nonnegative integer *r* such that + the Taylor series converges quickly, and then evaluating + `\exp(A/2^r)^{2^r}`. + If `\|A/2^r\| \le c` and `N \ge 2c`, we bound the entrywise error + when truncating the Taylor series before term `N` by `2 c^N / N!`. + diff --git a/doc/source/fmprb_mat.rst b/doc/source/fmprb_mat.rst index 377cd766..448cd0b7 100644 --- a/doc/source/fmprb_mat.rst +++ b/doc/source/fmprb_mat.rst @@ -236,3 +236,24 @@ Gaussian elimination and solving determinant of the remaining submatrix is bounded using Hadamard's inequality. + +Special functions +------------------------------------------------------------------------------- + +.. function:: void fmprb_mat_exp(fmprb_mat_t B, const fmprb_mat_t A, long prec) + + Sets *B* to the exponential of the matrix *A*, defined by the Taylor series + + .. math :: + + \exp(A) = \sum_{k=0}^{\infty} \frac{A^k}{k!}. + + The exponential function is evaluated using scaling followed by + rectangular splitting evaluation of the Taylor series. + Scaling amounts to picking a nonnegative integer *r* such that + the Taylor series converges quickly, and then evaluating + `\exp(A/2^r)^{2^r}`. + If `\|A/2^r\| \le c` and `N \ge 2c`, we bound the entrywise error + when truncating the Taylor series before term `N` by `2 c^N / N!`. + + diff --git a/fmpcb.h b/fmpcb.h index 98d35b72..1a360aec 100644 --- a/fmpcb.h +++ b/fmpcb.h @@ -249,6 +249,13 @@ fmpcb_trim(fmpcb_t z, const fmpcb_t x) fmprb_trim(fmpcb_imagref(z), fmpcb_imagref(x)); } +static __inline__ void +fmpcb_add_error_fmpr(fmpcb_t x, const fmpr_t err) +{ + fmprb_add_error_fmpr(fmpcb_realref(x), err); + fmprb_add_error_fmpr(fmpcb_imagref(x), err); +} + static __inline__ void fmpcb_get_abs_ubound_fmpr(fmpr_t u, const fmpcb_t z, long prec) { @@ -826,10 +833,7 @@ _fmpcb_vec_add_error_fmpr_vec(fmpcb_ptr res, fmpr_srcptr err, long len) { long i; for (i = 0; i < len; i++) - { - fmprb_add_error_fmpr(fmpcb_realref(res + i), err + i); - fmprb_add_error_fmpr(fmpcb_imagref(res + i), err + i); - } + fmpcb_add_error_fmpr(res + i, err + i); } static __inline__ void diff --git a/fmpcb_mat.h b/fmpcb_mat.h index 65dc46d7..3422eb0e 100644 --- a/fmpcb_mat.h +++ b/fmpcb_mat.h @@ -282,6 +282,9 @@ int fmpcb_mat_inv(fmpcb_mat_t X, const fmpcb_mat_t A, long prec); void fmpcb_mat_det(fmpcb_t det, const fmpcb_mat_t A, long prec); +/* Special functions */ + +void fmpcb_mat_exp(fmpcb_mat_t B, const fmpcb_mat_t A, long prec); #ifdef __cplusplus } diff --git a/fmpcb_mat/exp.c b/fmpcb_mat/exp.c new file mode 100644 index 00000000..517e0c7f --- /dev/null +++ b/fmpcb_mat/exp.c @@ -0,0 +1,199 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "double_extras.h" +#include "fmpcb_mat.h" + +long _fmprb_mat_exp_choose_N(const fmpr_t norm, long prec); + +void _fmprb_mat_exp_bound(fmpr_t err, const fmpr_t norm, long N); + +/* evaluates the truncated Taylor series (assumes no aliasing) */ +void +_fmpcb_mat_exp_taylor(fmpcb_mat_t S, const fmpcb_mat_t A, long N, long prec) +{ + if (N == 1) + { + fmpcb_mat_one(S); + } + else if (N == 2) + { + fmpcb_mat_one(S); + fmpcb_mat_add(S, S, A, prec); + } + else if (N == 3) + { + fmpcb_mat_t T; + fmpcb_mat_init(T, fmpcb_mat_nrows(A), fmpcb_mat_nrows(A)); + fmpcb_mat_mul(T, A, A, prec); + fmpcb_mat_scalar_mul_2exp_si(T, T, -1); + fmpcb_mat_add(S, A, T, prec); + fmpcb_mat_one(T); + fmpcb_mat_add(S, S, T, prec); + fmpcb_mat_clear(T); + } + else + { + long i, lo, hi, m, w, dim; + fmpcb_mat_struct * pows; + fmpcb_mat_t T, U; + fmpz_t c, f; + + dim = fmpcb_mat_nrows(A); + m = n_sqrt(N); + w = (N + m - 1) / m; + + fmpz_init(c); + fmpz_init(f); + pows = flint_malloc(sizeof(fmpcb_mat_t) * (m + 1)); + fmpcb_mat_init(T, dim, dim); + fmpcb_mat_init(U, dim, dim); + + for (i = 0; i <= m; i++) + { + fmpcb_mat_init(pows + i, dim, dim); + if (i == 0) + fmpcb_mat_one(pows + i); + else if (i == 1) + fmpcb_mat_set(pows + i, A); + else + fmpcb_mat_mul(pows + i, pows + i - 1, A, prec); + } + + fmpcb_mat_zero(S); + fmpz_one(f); + + for (i = w - 1; i >= 0; i--) + { + lo = i * m; + hi = FLINT_MIN(N - 1, lo + m - 1); + + fmpcb_mat_zero(T); + fmpz_one(c); + + while (hi >= lo) + { + fmpcb_mat_scalar_addmul_fmpz(T, pows + hi - lo, c, prec); + if (hi != 0) + fmpz_mul_ui(c, c, hi); + hi--; + } + + fmpcb_mat_mul(U, pows + m, S, prec); + fmpcb_mat_scalar_mul_fmpz(S, T, f, prec); + fmpcb_mat_add(S, S, U, prec); + fmpz_mul(f, f, c); + } + + fmpcb_mat_scalar_div_fmpz(S, S, f, prec); + + fmpz_clear(c); + fmpz_clear(f); + for (i = 0; i <= m; i++) + fmpcb_mat_clear(pows + i); + flint_free(pows); + fmpcb_mat_clear(T); + fmpcb_mat_clear(U); + } +} + +void +fmpcb_mat_exp(fmpcb_mat_t B, const fmpcb_mat_t A, long prec) +{ + long i, j, dim, wp, N, q, r; + fmpr_t norm, err; + fmpcb_mat_t T; + + dim = fmpcb_mat_nrows(A); + + if (dim != fmpcb_mat_ncols(A)) + { + printf("fmpcb_mat_exp: a square matrix is required!\n"); + abort(); + } + + if (dim == 0) + { + return; + } + else if (dim == 1) + { + fmpcb_exp(fmpcb_mat_entry(B, 0, 0), fmpcb_mat_entry(A, 0, 0), prec); + return; + } + + wp = prec + 3 * FLINT_BIT_COUNT(prec); + + fmpr_init(norm); + fmpr_init(err); + fmpcb_mat_init(T, dim, dim); + + fmpcb_mat_bound_inf_norm(norm, A, FMPRB_RAD_PREC); + + if (fmpr_is_zero(norm)) + { + fmpcb_mat_one(B); + } + else + { + r = fmpr_abs_bound_lt_2exp_si(norm); + q = pow(wp, 0.25); /* wanted magnitude */ + + if (r > 2 * wp) /* too big */ + r = 2 * wp; + else if (r < -q) /* tiny, no need to reduce */ + r = 0; + else + r += q; /* reduce to magnitude 2^(-r) */ + + fmpcb_mat_scalar_mul_2exp_si(T, A, -r); + fmpr_mul_2exp_si(norm, norm, -r); + + N = _fmprb_mat_exp_choose_N(norm, wp); + _fmprb_mat_exp_bound(err, norm, N); + + _fmpcb_mat_exp_taylor(B, T, N, wp); + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + fmpcb_add_error_fmpr(fmpcb_mat_entry(B, i, j), err); + + for (i = 0; i < r; i++) + { + fmpcb_mat_mul(T, B, B, wp); + fmpcb_mat_swap(T, B); + } + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + fmpcb_set_round(fmpcb_mat_entry(B, i, j), + fmpcb_mat_entry(B, i, j), prec); + } + + fmpr_clear(norm); + fmpr_clear(err); + fmpcb_mat_clear(T); +} + diff --git a/fmpcb_mat/test/t-exp.c b/fmpcb_mat/test/t-exp.c new file mode 100644 index 00000000..dc4c592c --- /dev/null +++ b/fmpcb_mat/test/t-exp.c @@ -0,0 +1,105 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_mat.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("exp...."); + fflush(stdout); + + flint_randinit(state); + + /* check exp(A)*exp(c*A) = exp((1+c)*A) */ + for (iter = 0; iter < 500; iter++) + { + fmpcb_mat_t A, E, F, EF, G; + fmpq_mat_t Q; + fmpcb_t c, d; + long n, qbits, prec; + + n = n_randint(state, 5); + qbits = 2 + n_randint(state, 300); + prec = 2 + n_randint(state, 300); + + fmpcb_init(c); + fmpcb_init(d); + fmpq_mat_init(Q, n, n); + fmpcb_mat_init(A, n, n); + fmpcb_mat_init(E, n, n); + fmpcb_mat_init(F, n, n); + fmpcb_mat_init(EF, n, n); + fmpcb_mat_init(G, n, n); + + fmpq_mat_randtest(Q, state, qbits); + fmpcb_mat_set_fmpq_mat(A, Q, prec); + + fmpcb_mat_exp(E, A, prec); + + fmpcb_randtest(c, state, prec, 10); + fmpcb_mat_scalar_mul_fmpcb(F, A, c, prec); + fmpcb_mat_exp(F, F, prec); + + fmpcb_add_ui(d, c, 1, prec); + fmpcb_mat_scalar_mul_fmpcb(G, A, d, prec); + fmpcb_mat_exp(G, G, prec); + + fmpcb_mat_mul(EF, E, F, prec); + + if (!fmpcb_mat_overlaps(EF, G)) + { + printf("FAIL\n\n"); + printf("n = %ld, prec = %ld\n", n, prec); + + printf("c = \n"); fmpcb_printd(c, 15); printf("\n\n"); + + printf("A = \n"); fmpcb_mat_printd(A, 15); printf("\n\n"); + printf("E = \n"); fmpcb_mat_printd(E, 15); printf("\n\n"); + printf("F = \n"); fmpcb_mat_printd(F, 15); printf("\n\n"); + printf("E*F = \n"); fmpcb_mat_printd(EF, 15); printf("\n\n"); + printf("G = \n"); fmpcb_mat_printd(G, 15); printf("\n\n"); + + abort(); + } + + fmpcb_clear(c); + fmpcb_clear(d); + fmpq_mat_clear(Q); + fmpcb_mat_clear(A); + fmpcb_mat_clear(E); + fmpcb_mat_clear(F); + fmpcb_mat_clear(EF); + fmpcb_mat_clear(G); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/fmprb_mat.h b/fmprb_mat.h index 6eddb217..9b61d23d 100644 --- a/fmprb_mat.h +++ b/fmprb_mat.h @@ -254,6 +254,10 @@ int fmprb_mat_inv(fmprb_mat_t X, const fmprb_mat_t A, long prec); void fmprb_mat_det(fmprb_t det, const fmprb_mat_t A, long prec); +/* Special functions */ + +void fmprb_mat_exp(fmprb_mat_t B, const fmprb_mat_t A, long prec); + #ifdef __cplusplus } #endif diff --git a/fmprb_mat/exp.c b/fmprb_mat/exp.c new file mode 100644 index 00000000..d69ac54f --- /dev/null +++ b/fmprb_mat/exp.c @@ -0,0 +1,254 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "double_extras.h" +#include "fmprb_mat.h" + +#define LOG2_OVER_E 0.25499459743395350926 + +void fmpr_gamma_ui_lbound(fmpr_t x, ulong n, long prec); + +long +_fmprb_mat_exp_choose_N(const fmpr_t norm, long prec) +{ + if (fmpr_is_special(norm) || fmpr_cmp_2exp_si(norm, 30) > 0) + { + return 1; + } + else if (fmpr_cmp_2exp_si(norm, -300) < 0) + { + long N = -fmpr_abs_bound_lt_2exp_si(norm); + + return (prec + N - 1) / N; + } + else + { + double c, t; + + c = fmpr_get_d(norm, FMPR_RND_UP); + + t = d_lambertw(prec * LOG2_OVER_E / c); + t = c * exp(t + 1.0); + + return FLINT_MIN((long) (t + 1.0), 2 * prec); + } +} + +void +_fmprb_mat_exp_bound(fmpr_t err, const fmpr_t norm, long N) +{ + fmpr_t t, u; + + fmpr_init(t); + fmpr_init(u); + + fmpr_set_si_2exp_si(t, N, -1); + + /* bound by geometric series when N >= 2*c <=> N/2 >= c */ + if (N > 0 && fmpr_cmp(t, norm) >= 0) + { + /* 2 c^N / N! */ + fmpr_pow_sloppy_ui(t, norm, N, FMPRB_RAD_PREC, FMPR_RND_UP); + fmpr_gamma_ui_lbound(u, N + 1, FMPRB_RAD_PREC); + fmpr_div(err, t, u, FMPRB_RAD_PREC, FMPR_RND_UP); + fmpr_mul_2exp_si(err, err, 1); + } + else + { + fmpr_pos_inf(err); + } + + fmpr_clear(t); + fmpr_clear(u); +} + + +/* evaluates the truncated Taylor series (assumes no aliasing) */ +void +_fmprb_mat_exp_taylor(fmprb_mat_t S, const fmprb_mat_t A, long N, long prec) +{ + if (N == 1) + { + fmprb_mat_one(S); + } + else if (N == 2) + { + fmprb_mat_one(S); + fmprb_mat_add(S, S, A, prec); + } + else if (N == 3) + { + fmprb_mat_t T; + fmprb_mat_init(T, fmprb_mat_nrows(A), fmprb_mat_nrows(A)); + fmprb_mat_mul(T, A, A, prec); + fmprb_mat_scalar_mul_2exp_si(T, T, -1); + fmprb_mat_add(S, A, T, prec); + fmprb_mat_one(T); + fmprb_mat_add(S, S, T, prec); + fmprb_mat_clear(T); + } + else + { + long i, lo, hi, m, w, dim; + fmprb_mat_struct * pows; + fmprb_mat_t T, U; + fmpz_t c, f; + + dim = fmprb_mat_nrows(A); + m = n_sqrt(N); + w = (N + m - 1) / m; + + fmpz_init(c); + fmpz_init(f); + pows = flint_malloc(sizeof(fmprb_mat_t) * (m + 1)); + fmprb_mat_init(T, dim, dim); + fmprb_mat_init(U, dim, dim); + + for (i = 0; i <= m; i++) + { + fmprb_mat_init(pows + i, dim, dim); + if (i == 0) + fmprb_mat_one(pows + i); + else if (i == 1) + fmprb_mat_set(pows + i, A); + else + fmprb_mat_mul(pows + i, pows + i - 1, A, prec); + } + + fmprb_mat_zero(S); + fmpz_one(f); + + for (i = w - 1; i >= 0; i--) + { + lo = i * m; + hi = FLINT_MIN(N - 1, lo + m - 1); + + fmprb_mat_zero(T); + fmpz_one(c); + + while (hi >= lo) + { + fmprb_mat_scalar_addmul_fmpz(T, pows + hi - lo, c, prec); + if (hi != 0) + fmpz_mul_ui(c, c, hi); + hi--; + } + + fmprb_mat_mul(U, pows + m, S, prec); + fmprb_mat_scalar_mul_fmpz(S, T, f, prec); + fmprb_mat_add(S, S, U, prec); + fmpz_mul(f, f, c); + } + + fmprb_mat_scalar_div_fmpz(S, S, f, prec); + + fmpz_clear(c); + fmpz_clear(f); + for (i = 0; i <= m; i++) + fmprb_mat_clear(pows + i); + flint_free(pows); + fmprb_mat_clear(T); + fmprb_mat_clear(U); + } +} + +void +fmprb_mat_exp(fmprb_mat_t B, const fmprb_mat_t A, long prec) +{ + long i, j, dim, wp, N, q, r; + fmpr_t norm, err; + fmprb_mat_t T; + + dim = fmprb_mat_nrows(A); + + if (dim != fmprb_mat_ncols(A)) + { + printf("fmprb_mat_exp: a square matrix is required!\n"); + abort(); + } + + if (dim == 0) + { + return; + } + else if (dim == 1) + { + fmprb_exp(fmprb_mat_entry(B, 0, 0), fmprb_mat_entry(A, 0, 0), prec); + return; + } + + wp = prec + 3 * FLINT_BIT_COUNT(prec); + + fmpr_init(norm); + fmpr_init(err); + fmprb_mat_init(T, dim, dim); + + fmprb_mat_bound_inf_norm(norm, A, FMPRB_RAD_PREC); + + if (fmpr_is_zero(norm)) + { + fmprb_mat_one(B); + } + else + { + r = fmpr_abs_bound_lt_2exp_si(norm); + q = pow(wp, 0.25); /* wanted magnitude */ + + if (r > 2 * wp) /* too big */ + r = 2 * wp; + else if (r < -q) /* tiny, no need to reduce */ + r = 0; + else + r += q; /* reduce to magnitude 2^(-r) */ + + fmprb_mat_scalar_mul_2exp_si(T, A, -r); + fmpr_mul_2exp_si(norm, norm, -r); + + N = _fmprb_mat_exp_choose_N(norm, wp); + _fmprb_mat_exp_bound(err, norm, N); + + _fmprb_mat_exp_taylor(B, T, N, wp); + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + fmprb_add_error_fmpr(fmprb_mat_entry(B, i, j), err); + + for (i = 0; i < r; i++) + { + fmprb_mat_mul(T, B, B, wp); + fmprb_mat_swap(T, B); + } + + for (i = 0; i < dim; i++) + for (j = 0; j < dim; j++) + fmprb_set_round(fmprb_mat_entry(B, i, j), + fmprb_mat_entry(B, i, j), prec); + } + + fmpr_clear(norm); + fmpr_clear(err); + fmprb_mat_clear(T); +} + diff --git a/fmprb_mat/test/t-exp.c b/fmprb_mat/test/t-exp.c new file mode 100644 index 00000000..51ebc7f1 --- /dev/null +++ b/fmprb_mat/test/t-exp.c @@ -0,0 +1,105 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "fmprb_mat.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("exp...."); + fflush(stdout); + + flint_randinit(state); + + /* check exp(A)*exp(c*A) = exp((1+c)*A) */ + for (iter = 0; iter < 1000; iter++) + { + fmprb_mat_t A, E, F, EF, G; + fmpq_mat_t Q; + fmprb_t c, d; + long n, qbits, prec; + + n = n_randint(state, 5); + qbits = 2 + n_randint(state, 300); + prec = 2 + n_randint(state, 300); + + fmprb_init(c); + fmprb_init(d); + fmpq_mat_init(Q, n, n); + fmprb_mat_init(A, n, n); + fmprb_mat_init(E, n, n); + fmprb_mat_init(F, n, n); + fmprb_mat_init(EF, n, n); + fmprb_mat_init(G, n, n); + + fmpq_mat_randtest(Q, state, qbits); + fmprb_mat_set_fmpq_mat(A, Q, prec); + + fmprb_mat_exp(E, A, prec); + + fmprb_randtest(c, state, prec, 10); + fmprb_mat_scalar_mul_fmprb(F, A, c, prec); + fmprb_mat_exp(F, F, prec); + + fmprb_add_ui(d, c, 1, prec); + fmprb_mat_scalar_mul_fmprb(G, A, d, prec); + fmprb_mat_exp(G, G, prec); + + fmprb_mat_mul(EF, E, F, prec); + + if (!fmprb_mat_overlaps(EF, G)) + { + printf("FAIL\n\n"); + printf("n = %ld, prec = %ld\n", n, prec); + + printf("c = \n"); fmprb_printd(c, 15); printf("\n\n"); + + printf("A = \n"); fmprb_mat_printd(A, 15); printf("\n\n"); + printf("E = \n"); fmprb_mat_printd(E, 15); printf("\n\n"); + printf("F = \n"); fmprb_mat_printd(F, 15); printf("\n\n"); + printf("E*F = \n"); fmprb_mat_printd(EF, 15); printf("\n\n"); + printf("G = \n"); fmprb_mat_printd(G, 15); printf("\n\n"); + + abort(); + } + + fmprb_clear(c); + fmprb_clear(d); + fmpq_mat_clear(Q); + fmprb_mat_clear(A); + fmprb_mat_clear(E); + fmprb_mat_clear(F); + fmprb_mat_clear(EF); + fmprb_mat_clear(G); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +