Skip to content
Browse files

Implement Dixon's p-adic algorithm and fast matrix rational reconstru…

…ction for nonsingular fmpz_mat solving (initial code)
  • Loading branch information...
1 parent 7692a85 commit efbe2fa3dc3acd9a441b594b7e7f7a2fd40740f3 @fredrik-johansson fredrik-johansson committed Apr 4, 2011
View
5 fmpq_mat.h
@@ -71,9 +71,14 @@ void fmpq_mat_get_fmpz_mat_rowwise_2(fmpz_mat_t num, fmpz_mat_t num2,
void fmpq_mat_mul_direct(fmpq_mat_t C, const fmpq_mat_t A, const fmpq_mat_t B);
void fmpq_mat_mul_cleared(fmpq_mat_t C, const fmpq_mat_t A, const fmpq_mat_t B);
void fmpq_mat_mul(fmpq_mat_t C, const fmpq_mat_t A, const fmpq_mat_t B);
+void fmpq_mat_mul_fmpz_mat(fmpq_mat_t C, const fmpq_mat_t A, const fmpz_mat_t B);
+void fmpq_mat_mul_r_fmpz_mat(fmpq_mat_t C, const fmpz_mat_t A, const fmpq_mat_t B);
void fmpq_mat_det(fmpq_t det, fmpq_mat_t mat);
void fmpq_mat_solve_mat(fmpq_mat_t X, const fmpq_mat_t A, const fmpq_mat_t B);
+
+int fmpq_mat_set_fmpz_mat_mod_fmpz(fmpq_mat_t X, const fmpz_mat_t Xmod, const fmpz_t mod);
+
void fmpq_mat_inv(fmpq_mat_t B, const fmpq_mat_t A);
void fmpq_mat_rref(fmpq_mat_t B, const fmpq_mat_t A);
View
68 fmpq_mat/mul_fmpz_mat.c
@@ -0,0 +1,68 @@
+/*=============================================================================
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2010 William Hart
+ Copyright (C) 2011 Fredrik Johansson
+
+******************************************************************************/
+
+#include <mpir.h>
+#include "flint.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "fmpz_mat.h"
+#include "fmpq.h"
+#include "fmpq_mat.h"
+
+
+void
+fmpq_mat_mul_fmpz_mat(fmpq_mat_t C, const fmpq_mat_t A, const fmpz_mat_t B)
+{
+ long i, j;
+
+ fmpz_mat_t Aclear;
+ fmpz_mat_t Cclear;
+
+ fmpz * Aden;
+
+ fmpz_mat_init(Aclear, A->r, A->c);
+ fmpz_mat_init(Cclear, A->r, B->c);
+
+ Aden = _fmpz_vec_init(A->r);
+
+ fmpq_mat_get_fmpz_mat_rowwise(Aclear, Aden, A);
+ fmpz_mat_mul(Cclear, Aclear, B);
+
+ for (i = 0; i < C->r; i++)
+ {
+ for (j = 0; j < C->c; j++)
+ {
+ fmpz_set(fmpq_mat_entry_num(C, i, j), fmpz_mat_entry(Cclear, i, j));
+ fmpz_set(fmpq_mat_entry_den(C, i, j), Aden + i);
+ fmpq_canonicalise(fmpq_mat_entry(C, i, j));
+ }
+ }
+
+ fmpz_mat_clear(Aclear);
+ fmpz_mat_clear(Cclear);
+
+ _fmpz_vec_clear(Aden, A->r);
+}
View
68 fmpq_mat/mul_r_fmpz_mat.c
@@ -0,0 +1,68 @@
+/*=============================================================================
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2010 William Hart
+ Copyright (C) 2011 Fredrik Johansson
+
+******************************************************************************/
+
+#include <mpir.h>
+#include "flint.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "fmpz_mat.h"
+#include "fmpq.h"
+#include "fmpq_mat.h"
+
+
+void
+fmpq_mat_mul_r_fmpz_mat(fmpq_mat_t C, const fmpz_mat_t A, const fmpq_mat_t B)
+{
+ long i, j;
+
+ fmpz_mat_t Bclear;
+ fmpz_mat_t Cclear;
+
+ fmpz * Bden;
+
+ fmpz_mat_init(Bclear, B->r, B->c);
+ fmpz_mat_init(Cclear, A->r, B->c);
+
+ Bden = _fmpz_vec_init(B->c);
+
+ fmpq_mat_get_fmpz_mat_colwise(Bclear, Bden, B);
+ fmpz_mat_mul(Cclear, A, Bclear);
+
+ for (i = 0; i < C->r; i++)
+ {
+ for (j = 0; j < C->c; j++)
+ {
+ fmpz_set(fmpq_mat_entry_num(C, i, j), fmpz_mat_entry(Cclear, i, j));
+ fmpz_set(fmpq_mat_entry_den(C, i, j), Bden + j);
+ fmpq_canonicalise(fmpq_mat_entry(C, i, j));
+ }
+ }
+
+ fmpz_mat_clear(Bclear);
+ fmpz_mat_clear(Cclear);
+
+ _fmpz_vec_clear(Bden, B->c);
+}
View
82 fmpq_mat/set_fmpz_mat_mod_fmpz.c
@@ -0,0 +1,82 @@
+/*=============================================================================
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2011 Fredrik Johansson
+
+******************************************************************************/
+
+#include <mpir.h>
+#include "flint.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "fmpz_mat.h"
+#include "fmpq.h"
+#include "fmpq_mat.h"
+
+int
+fmpq_mat_set_fmpz_mat_mod_fmpz(fmpq_mat_t X,
+ const fmpz_mat_t Xmod, const fmpz_t mod)
+{
+ fmpz_t num, den, t, u, d;
+ long i, j;
+
+ int success = 1;
+
+ fmpz_init(num);
+ fmpz_init(den);
+ fmpz_init(d);
+ fmpz_init(t);
+ fmpz_init(u);
+
+ fmpz_set_ui(d, 1UL);
+
+ for (i = 0; i < Xmod->r; i++)
+ {
+ for (j = 0; j < Xmod->c; j++)
+ {
+ /* TODO: handle various special cases efficiently; zeros,
+ small integers, etc. */
+ fmpz_mul(t, d, fmpz_mat_entry(Xmod, i, j));
+ fmpz_fdiv_qr(u, t, t, mod);
+
+ success = _fmpq_reconstruct_fmpz(num, den, t, mod);
+
+ fmpz_mul(den, den, d);
+ fmpz_set(d, den);
+
+ if (!success)
+ goto cleanup;
+
+ fmpz_set(fmpq_mat_entry_num(X, i, j), num);
+ fmpz_set(fmpq_mat_entry_den(X, i, j), den);
+ fmpq_canonicalise(fmpq_mat_entry(X, i, j));
+ }
+ }
+
+cleanup:
+ fmpz_clear(num);
+ fmpz_clear(den);
+ fmpz_clear(d);
+ fmpz_clear(t);
+ fmpz_clear(u);
+
+ return success;
+}
View
4 fmpz_mat.h
@@ -140,6 +140,8 @@ void fmpz_mat_scalar_divexact_fmpz(fmpz_mat_t B, const fmpz_mat_t A, const fmpz_
void fmpz_mat_scalar_divexact_si(fmpz_mat_t B, const fmpz_mat_t A, long c);
void fmpz_mat_scalar_divexact_ui(fmpz_mat_t B, const fmpz_mat_t A, ulong c);
+void fmpz_mat_scalar_mod_fmpz(fmpz_mat_t B, const fmpz_mat_t A, const fmpz_t m);
+
/* Multiplication */
void fmpz_mat_mul(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B);
@@ -183,6 +185,8 @@ void fmpz_mat_solve_fraction_free_LU(fmpz * x, fmpz_t den, const fmpz_mat_t A, c
void _fmpz_mat_solve_fraction_free_LU_precomp(fmpz * b, const fmpz_mat_t LU);
void fmpz_mat_solve_mat(fmpz_mat_t X, fmpz_t den, const fmpz_mat_t A, const fmpz_mat_t B);
+void fmpz_mat_solve_dixon(fmpz_mat_t X, fmpz_t mod, const fmpz_mat_t A, const fmpz_mat_t B);
+
/* Kernel */
long fmpz_mat_kernel(fmpz_mat_t res, const fmpz_mat_t mat);
View
44 fmpz_mat/scalar_mod_fmpz.c
@@ -0,0 +1,44 @@
+/*=============================================================================
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2011 Fredrik Johansson
+
+******************************************************************************/
+
+#include <stdlib.h>
+#include "flint.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "fmpz_mat.h"
+
+void
+fmpz_mat_scalar_mod_fmpz(fmpz_mat_t B, const fmpz_mat_t A, const fmpz_t m)
+{
+ long i, j;
+ fmpz_t t;
+ fmpz_init(t);
+
+ for (i = 0; i < A->r; i++)
+ for (j = 0; j < A->c; j++)
+ fmpz_fdiv_qr(t, fmpz_mat_entry(B,i,j), fmpz_mat_entry(A,i,j), m);
+
+ fmpz_init(t);
+}
View
133 fmpz_mat/solve_dixon.c
@@ -0,0 +1,133 @@
+/*=============================================================================
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2011 Fredrik Johansson
+
+******************************************************************************/
+
+#include <stdlib.h>
+#include "flint.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "fmpz_mat.h"
+#include "nmod_mat.h"
+#include "ulong_extras.h"
+
+
+void
+fmpz_mat_solve_dixon(fmpz_mat_t X, fmpz_t mod,
+ const fmpz_mat_t A, const fmpz_mat_t B)
+{
+ long n, cols;
+ mp_limb_t p;
+ fmpz_t N, D, bound, ppow;
+
+ fmpz_mat_t x, y, d, Ay;
+ nmod_mat_t Ainv, dmod, ymod;
+
+ n = A->r;
+ cols = B->c;
+
+ if (!fmpz_mat_is_square(A))
+ {
+ printf("fmpz_mat_solve_dixon: nonsquare system matrix");
+ abort();
+ }
+
+ if (fmpz_mat_is_empty(A) || fmpz_mat_is_empty(B))
+ return;
+
+ p = n_nextprime(1UL << (FLINT_BITS - FLINT_BIT_COUNT(n) - 2), 0);
+
+ fmpz_init(bound);
+ fmpz_init(N);
+ fmpz_init(D);
+
+ fmpz_init(ppow);
+ fmpz_mat_init(x, n, cols);
+ fmpz_mat_init(y, n, cols);
+ fmpz_mat_init(Ay, n, cols);
+ fmpz_mat_init_set(d, B);
+
+ while (1)
+ {
+ nmod_mat_init(Ainv, n, n, p);
+ fmpz_mat_get_nmod_mat(Ainv, A);
+ if (nmod_mat_inv(Ainv, Ainv))
+ break;
+ p = n_nextprime(p, 0);
+ nmod_mat_clear(Ainv);
+ }
+
+ nmod_mat_init(dmod, n, cols, p);
+ nmod_mat_init(ymod, n, cols, p);
+
+ fmpz_mat_solve_bound(N, D, A, B);
+
+ /* TODO: changing the balance here, we would have to pass on N and D */
+ if (fmpz_cmpabs(N, D) < 0)
+ fmpz_mul(bound, D, D);
+ else
+ fmpz_mul(bound, N, N);
+ /* fmpz_mul(bound, N, D); */
+
+ fmpz_mul_ui(bound, bound, 2UL);
+
+ fmpz_set_ui(ppow, 1UL);
+
+ while (fmpz_cmp(ppow, bound) <= 0)
+ {
+ /* y = A^(-1) * d (mod p) */
+ fmpz_mat_get_nmod_mat(dmod, d);
+ nmod_mat_mul(ymod, Ainv, dmod);
+
+ /* x = x + y * p^i [= A^(-1) * b mod p^(i+1)] */
+ fmpz_mat_scalar_addmul_nmod_mat_fmpz(x, ymod, ppow);
+
+ /* ppow = p^(i+1) */
+ fmpz_mul_ui(ppow, ppow, p);
+ if (fmpz_cmp(ppow, bound) > 0)
+ break;
+
+ /* d = (d - Ay) / p */
+ /* TODO: implement fmpz_mat_mul_nmod_mat to avoid y? */
+ fmpz_mat_set_nmod_mat(y, ymod);
+ fmpz_mat_mul(Ay, A, y);
+ fmpz_mat_sub(d, d, Ay);
+ fmpz_mat_scalar_divexact_ui(d, d, p);
+ }
+
+ fmpz_set(mod, ppow);
+ fmpz_mat_set(X, x);
+
+ fmpz_clear(bound);
+ fmpz_clear(N);
+ fmpz_clear(D);
+
+ fmpz_clear(ppow);
+ fmpz_mat_clear(x);
+ fmpz_mat_clear(y);
+ fmpz_mat_clear(d);
+ fmpz_mat_clear(Ay);
+ nmod_mat_clear(Ainv);
+ nmod_mat_clear(ymod);
+ nmod_mat_clear(dmod);
+}
View
97 fmpz_mat/test/t-solve_dixon.c
@@ -0,0 +1,97 @@
+/*=============================================================================
+
+ This file is part of FLINT.
+
+ FLINT 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.
+
+ FLINT 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 FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2010 Fredrik Johansson
+
+******************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <mpir.h>
+#include "flint.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+#include "fmpz_mat.h"
+#include "ulong_extras.h"
+
+int
+main(void)
+{
+ fmpz_mat_t A, X, B, AX, Bm;
+ fmpz_t mod;
+ flint_rand_t state;
+ long i, m, n;
+
+ printf("solve_dixon....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (i = 0; i < 1000; i++)
+ {
+ m = n_randint(state, 30);
+ n = n_randint(state, 30);
+
+ fmpz_mat_init(A, m, m);
+ fmpz_mat_init(B, m, n);
+ fmpz_mat_init(Bm, m, n);
+ fmpz_mat_init(X, m, n);
+ fmpz_mat_init(AX, m, n);
+ fmpz_init(mod);
+
+ fmpz_mat_randrank(A, state, m, 1+n_randint(state, 2)*n_randint(state, 100));
+ fmpz_mat_randtest(B, state, 1+n_randint(state, 2)*n_randint(state, 100));
+
+ /* Dense */
+ if (n_randint(state, 2))
+ fmpz_mat_randops(A, state, 1+n_randint(state, 1 + m*m));
+
+ fmpz_mat_solve_dixon(X, mod, A, B);
+ fmpz_mat_mul(AX, A, X);
+
+ fmpz_mat_scalar_mod_fmpz(AX, AX, mod);
+ fmpz_mat_scalar_mod_fmpz(Bm, B, mod);
+
+ if (!fmpz_mat_equal(AX, Bm))
+ {
+ printf("FAIL:\n");
+ printf("AX != B!\n");
+ printf("A:\n"), fmpz_mat_print_pretty(A), printf("\n");
+ printf("B:\n"), fmpz_mat_print_pretty(B), printf("\n");
+ printf("X:\n"), fmpz_mat_print_pretty(X), printf("\n");
+ printf("mod = "), fmpz_print(mod), printf("\n");
+ printf("AX:\n"), fmpz_mat_print_pretty(AX), printf("\n");
+ abort();
+ }
+
+ fmpz_mat_clear(A);
+ fmpz_mat_clear(B);
+ fmpz_mat_clear(Bm);
+ fmpz_mat_clear(X);
+ fmpz_mat_clear(AX);
+ fmpz_clear(mod);
+ }
+
+ flint_randclear(state);
+ _fmpz_cleanup();
+ printf("PASS\n");
+ return 0;
+}
View
2 fmpz_mat/test/t-solve_mat.c
@@ -67,7 +67,7 @@ main(void)
fmpz_mat_solve_mat(X, den, A, B);
fmpz_mat_mul(AX, A, X);
- _fmpz_vec_scalar_divexact_fmpz(AX->entries, AX->entries, m*n, den);
+ fmpz_mat_scalar_divexact_fmpz(AX, AX, den);
if (!fmpz_mat_equal(AX, B))
{

0 comments on commit efbe2fa

Please sign in to comment.
Something went wrong with that request. Please try again.