Merge pull request #95 from argriffing/fmpz-mat-extras

MAINT: move transitive_closure to a new fmpz_mat_extras module

Makefile.in

@@ -14,7 +14,7 @@ AT=@
 
 BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \
    acb_calc acb_hypgeom acb_modular acb_dirichlet fmprb bernoulli hypgeom \
-   fmpz_extras partitions \
+   fmpz_extras fmpz_mat_extras partitions \
    $(EXTRA_BUILD_DIRS)
 
 TEMPLATE_DIRS = 

acb_mat/exp.c

@@ -25,9 +25,9 @@
 
 #include "double_extras.h"
 #include "acb_mat.h"
+#include "fmpz_mat_extras.h"
 
 slong _arb_mat_exp_choose_N(const mag_t norm, slong prec);
-void _fmpz_mat_transitive_closure(fmpz_mat_t A, const fmpz_mat_t B);
 
 int
 _acb_mat_is_diagonal(const acb_mat_t A)
@@ -68,7 +68,7 @@ _acb_mat_exp_get_structure(fmpz_mat_t C, const acb_mat_t A)
             }
         }
     }
-    _fmpz_mat_transitive_closure(C, C);
+    fmpz_mat_transitive_closure(C, C);
 }
 
 void

arb_mat/exp.c

@@ -26,45 +26,10 @@
 #include "fmpz_mat.h"
 #include "double_extras.h"
 #include "arb_mat.h"
+#include "fmpz_mat_extras.h"
 
 #define LOG2_OVER_E 0.25499459743395350926
 
-
-/* Warshall's algorithm */
-void
-_fmpz_mat_transitive_closure(fmpz_mat_t B, fmpz_mat_t A)
-{
-    slong k, i, j, dim;
-    dim = fmpz_mat_nrows(A);
-
-    if (dim != fmpz_mat_ncols(A))
-    {
-        flint_printf("_fmpz_mat_transitive_closure: a square matrix is required!\n");
-        abort();
-    }
-
-    if (A != B)
-    {
-        fmpz_mat_set(B, A);
-    }
-
-    for (k = 0; k < dim; k++)
-    {
-        for (i = 0; i < dim; i++)
-        {
-            for (j = 0; j < dim; j++)
-            {
-                if (fmpz_is_zero(fmpz_mat_entry(B, i, j)) &&
-                    !fmpz_is_zero(fmpz_mat_entry(B, i, k)) &&
-                    !fmpz_is_zero(fmpz_mat_entry(B, k, j)))
-                {
-                    fmpz_one(fmpz_mat_entry(B, i, j));
-                }
-            }
-        }
-    }
-}
-
 int
 _arb_mat_is_diagonal(const arb_mat_t A)
 {
@@ -275,7 +240,7 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
         {
             fmpz_mat_init(S, dim, dim);
             arb_mat_entrywise_not_is_zero(S, A);
-            _fmpz_mat_transitive_closure(S, S);
+            fmpz_mat_transitive_closure(S, S);
         }
 
         q = pow(wp, 0.25);  /* wanted magnitude */

fmpz_mat_extras.h

@@ -0,0 +1,43 @@
+/*=============================================================================
+
+    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) 2016 Arb authors
+
+******************************************************************************/
+
+#ifndef FMPZ_MAT_EXTRAS_H
+#define FMPZ_MAT_EXTRAS_H
+
+#include "flint.h"
+#include "fmpz_mat.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void fmpz_mat_transitive_closure(fmpz_mat_t dest, const fmpz_mat_t src);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
+

fmpz_mat_extras/test/t-transitive_closure.c

@@ -0,0 +1,265 @@
+/*=============================================================================
+
+    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) 2016 Arb authors
+
+******************************************************************************/
+
+#include "fmpz_mat_extras.h"
+#include "perm.h"
+
+/* transitive closure can only turn zeros into ones */
+int
+_transitive_closure_is_ok_entrywise(const fmpz_mat_t X, const fmpz_mat_t Y)
+{
+    slong i, j;
+    if (fmpz_mat_nrows(X) != fmpz_mat_nrows(Y) ||
+        fmpz_mat_ncols(X) != fmpz_mat_ncols(Y))
+    {
+        return 0;
+    }
+    for (i = 0; i < fmpz_mat_nrows(X); i++)
+    {
+        for (j = 0; j < fmpz_mat_ncols(X); j++)
+        {
+            if (!fmpz_equal(
+                        fmpz_mat_entry(X, i, j),
+                        fmpz_mat_entry(Y, i, j)))
+            {
+                if (!fmpz_is_zero(fmpz_mat_entry(X, i, j)))
+                {
+                    return 0;
+                }
+                if (!fmpz_is_one(fmpz_mat_entry(Y, i, j)))
+                {
+                    return 0;
+                }
+            }
+        }
+    }
+    return 1;
+}
+
+/* permute rows and columns of a square matrix */
+void
+_fmpz_mat_permute(fmpz_mat_t B, const fmpz_mat_t A, const slong *perm)
+{
+    slong n, i, j;
+    if (!fmpz_mat_is_square(A)) abort();
+    if (A == B) abort();
+    n = fmpz_mat_nrows(A);
+    for (i = 0; i < n; i++)
+    {
+        for (j = 0; j < n; j++)
+        {
+            fmpz_set(
+                    fmpz_mat_entry(B, perm[i], perm[j]),
+                    fmpz_mat_entry(A, i, j));
+        }
+    }
+}
+
+/* remove all non-sparsity information from an fmpz matrix */
+void
+_fmpz_mat_entrywise_not_is_zero(fmpz_mat_t B, const fmpz_mat_t A)
+{
+    slong i, j;
+    for (i = 0; i < fmpz_mat_nrows(A); i++)
+    {
+        for (j = 0; j < fmpz_mat_ncols(A); j++)
+        {
+            if (fmpz_is_zero(fmpz_mat_entry(A, i, j)))
+            {
+                fmpz_zero(fmpz_mat_entry(B, i, j));
+            }
+            else
+            {
+                fmpz_one(fmpz_mat_entry(B, i, j));
+            }
+        }
+    }
+}
+
+/* this is not efficient */
+void
+_brute_force_transitive_closure(fmpz_mat_t B, const fmpz_mat_t A)
+{
+    slong n, k;
+    fmpz_mat_t S, curr, accum;
+
+    n = fmpz_mat_nrows(A);
+    fmpz_mat_init(S, n, n);
+    fmpz_mat_init(curr, n, n);
+    fmpz_mat_init(accum, n, n);
+    _fmpz_mat_entrywise_not_is_zero(S, A);
+    fmpz_mat_one(curr);
+    fmpz_mat_zero(accum);
+    for (k = 0; k < n; k++)
+    {
+        fmpz_mat_mul(curr, curr, S);
+        fmpz_mat_add(accum, accum, curr);
+    }
+    fmpz_mat_clear(S);
+    fmpz_mat_clear(curr);
+    {
+        slong i, j;
+        for (i = 0; i < n; i++)
+        {
+            for (j = 0; j < n; j++)
+            {
+                if (fmpz_is_zero(fmpz_mat_entry(A, i, j)) &&
+                    !fmpz_is_zero(fmpz_mat_entry(accum, i, j)))
+                {
+                    fmpz_one(fmpz_mat_entry(B, i, j));
+                }
+                else
+                {
+                    fmpz_set(fmpz_mat_entry(B, i, j),
+                             fmpz_mat_entry(A, i, j));
+                }
+            }
+        }
+    }
+    fmpz_mat_clear(accum);
+}
+
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("transitive_closure....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        slong m;
+        fmpz_mat_t A, B, C, D;
+
+        m = n_randint(state, 50);
+
+        fmpz_mat_init(A, m, m);
+        fmpz_mat_init(B, m, m);
+        fmpz_mat_init(C, m, m);
+        fmpz_mat_init(D, m, m);
+
+        fmpz_mat_randtest(A, state, n_randint(state, 20) + 1);
+        fmpz_mat_randtest(B, state, n_randint(state, 20) + 1);
+
+        fmpz_mat_transitive_closure(B, A);
+
+        /* test local properties of the closure */
+        if (!_transitive_closure_is_ok_entrywise(A, B))
+        {
+            flint_printf("FAIL (entrywise)\n");
+            fmpz_mat_print_pretty(A); flint_printf("\n\n");
+            fmpz_mat_print_pretty(B); flint_printf("\n\n");
+            abort();
+        }
+
+        /* test aliasing */
+        {
+            fmpz_mat_set(C, A);
+            fmpz_mat_transitive_closure(C, C);
+            if (!fmpz_mat_equal(B, C))
+            {
+                flint_printf("FAIL (aliasing)\n");
+                fmpz_mat_print_pretty(A); flint_printf("\n\n");
+                fmpz_mat_print_pretty(B); flint_printf("\n\n");
+                fmpz_mat_print_pretty(C); flint_printf("\n\n");
+                abort();
+            }
+        }
+
+        /* test commutativity of permutation with transitive closure */
+        {
+            slong *perm;
+            perm = flint_malloc(m * sizeof(slong));
+            _perm_randtest(perm, m, state);
+
+            /* C is the transitive closure of the permutation of A */
+            fmpz_mat_randtest(C, state, n_randint(state, 20) + 1);
+            _fmpz_mat_permute(C, A, perm);
+            fmpz_mat_transitive_closure(C, C);
+
+            /* D is the permutation of the transitive closure of A */
+            fmpz_mat_randtest(D, state, n_randint(state, 20) + 1);
+            _fmpz_mat_permute(D, B, perm);
+
+            if (!fmpz_mat_equal(C, D))
+            {
+                flint_printf("FAIL (commutativity with permutation)\n");
+                fmpz_mat_print_pretty(A); flint_printf("\n\n");
+                fmpz_mat_print_pretty(B); flint_printf("\n\n");
+                fmpz_mat_print_pretty(C); flint_printf("\n\n");
+                fmpz_mat_print_pretty(D); flint_printf("\n\n");
+                abort();
+            }
+            flint_free(perm);
+        }
+
+        fmpz_mat_clear(A);
+        fmpz_mat_clear(B);
+        fmpz_mat_clear(C);
+        fmpz_mat_clear(D);
+    }
+
+    /* check transitive closure using brute force with smallish matrices */
+    for (iter = 0; iter < 1000; iter++)
+    {
+        slong m;
+        fmpz_mat_t A, B, C;
+
+        m = n_randint(state, 10);
+
+        fmpz_mat_init(A, m, m);
+        fmpz_mat_init(B, m, m);
+        fmpz_mat_init(C, m, m);
+
+        fmpz_mat_randtest(A, state, n_randint(state, 20) + 1);
+        fmpz_mat_randtest(B, state, n_randint(state, 20) + 1);
+        fmpz_mat_randtest(C, state, n_randint(state, 20) + 1);
+
+        fmpz_mat_transitive_closure(B, A);
+        _brute_force_transitive_closure(C, A);
+
+        if (!fmpz_mat_equal(B, C))
+        {
+            flint_printf("FAIL (brute force)\n");
+            fmpz_mat_print_pretty(A); flint_printf("\n\n");
+            fmpz_mat_print_pretty(B); flint_printf("\n\n");
+            fmpz_mat_print_pretty(C); flint_printf("\n\n");
+            abort();
+        }
+
+        fmpz_mat_clear(A);
+        fmpz_mat_clear(B);
+        fmpz_mat_clear(C);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmpz_mat_extras/transitive_closure.c

@@ -0,0 +1,65 @@
+/*=============================================================================
+
+    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) 2016 Arb authors
+
+******************************************************************************/
+
+#include "fmpz_mat_extras.h"
+
+/* Warshall's algorithm */
+void
+fmpz_mat_transitive_closure(fmpz_mat_t dest, const fmpz_mat_t src)
+{
+    slong k, i, j, dim;
+
+    if (fmpz_mat_nrows(dest) != fmpz_mat_nrows(src) ||
+        fmpz_mat_ncols(dest) != fmpz_mat_ncols(src))
+    {
+        flint_printf("fmpz_mat_transitive_closure: incompatible dimensions\n");
+        abort();
+    }
+
+    dim = fmpz_mat_nrows(src);
+    if (dim != fmpz_mat_ncols(src))
+    {
+        flint_printf("fmpz_mat_transitive_closure: a square matrix is required!\n");
+        abort();
+    }
+
+    fmpz_mat_set(dest, src);
+
+    for (k = 0; k < dim; k++)
+    {
+        for (i = 0; i < dim; i++)
+        {
+            for (j = 0; j < dim; j++)
+            {
+                if (fmpz_is_zero(fmpz_mat_entry(dest, i, j)) &&
+                    !fmpz_is_zero(fmpz_mat_entry(dest, i, k)) &&
+                    !fmpz_is_zero(fmpz_mat_entry(dest, k, j)))
+                {
+                    fmpz_one(fmpz_mat_entry(dest, i, j));
+                }
+            }
+        }
+    }
+}