From a107ad7a6e4f1caa6d1dba789eb29e61fbea407a Mon Sep 17 00:00:00 2001 From: Dennis Sweeney Date: Tue, 10 Mar 2026 16:15:20 -0400 Subject: [PATCH 1/4] decompose_relations_among: get the easy relations out of the way and isolate the hard subproblems --- _mutable_lattice.c | 548 +++++++++++++++++++++++++++++++++++- mutable_lattice/__init__.py | 15 +- mutable_lattice/test.py | 198 +++++++++++-- pyproject.toml | 2 +- 4 files changed, 723 insertions(+), 40 deletions(-) diff --git a/_mutable_lattice.c b/_mutable_lattice.c index 10d9a77..c0b0159 100644 --- a/_mutable_lattice.c +++ b/_mutable_lattice.c @@ -258,6 +258,7 @@ unpack_integer(TagInt t) #define TagInt_ONE ((TagInt) {.bits = 2}) #define TagInt_ZERO ((TagInt) {.bits = 0}) +#define TagInt_NEGATIVE_ONE ((TagInt) {.bits = -2}) static inline bool TagInt_is_zero(TagInt t) { @@ -3539,19 +3540,16 @@ Lattice_intersection(PyObject *self, PyObject *other) } static inline Py_ssize_t -uf_find(Py_ssize_t *parent, Py_ssize_t x0) +uf_find(Py_ssize_t *parent, Py_ssize_t x) { - Py_ssize_t x = x0, px; + // path halving + Py_ssize_t px; while ((px = parent[x]) != x) { - x = px; + Py_ssize_t ppx = parent[px]; + parent[x] = ppx; + x = ppx; } - Py_ssize_t root = x; - x = x0; - while ((px = parent[x]) != x) { - parent[x] = root; - x = px; - } - return root; + return x; } static inline void @@ -3682,8 +3680,7 @@ Lattice_decompose(PyObject *self, PyObject *const *args, Py_ssize_t nargs) } // flatten for (Py_ssize_t j = 0; j < N; j++) { - (void)uf_find(uf_parent, j); - assert(uf_parent[uf_parent[j]] == uf_parent[j]); + uf_parent[j] = uf_find(uf_parent, j); } // find num_components and root_to_component_index @@ -4007,7 +4004,7 @@ static PyTypeObject Lattice_Type = { }; /*********************************************************************/ -/* relations_among, transpose */ +/* relations_among, decompose_relations_among, transpose */ /*********************************************************************/ static PyObject * @@ -4193,6 +4190,528 @@ transpose(PyObject *Py_UNUSED(mod), PyObject *const *args, Py_ssize_t nargs) return result; } +static bool +Vector_cmp(TagInt *a, TagInt *b, Py_ssize_t N, int *res) +{ + // Compare two vectors of equal length. + // This isn't a useful order--it just has to be a total order + // so that we can efficiently deduplicate + for (Py_ssize_t j = 0; j < N; j++) { + if (a[j].bits == b[j].bits) { + continue; + } + if (TagInt_is_pointer(a[j])) { + if (TagInt_is_pointer(b[j])) { + PyObject *aj = untag_pointer(a[j]); + PyObject *bj = untag_pointer(b[j]); + int eq = PyObject_RichCompareBool(aj, bj, Py_EQ); + if (eq == -1) { + return true; + } + if (eq) { + continue; + } + int lt = PyObject_RichCompareBool(aj, bj, Py_LT); + if (lt == -1) { + return true; + } + *res = lt ? -1 : 1; + return false; + } + else { + *res = 1; + return false; + } + } + else { + if (TagInt_is_pointer(b[j])) { + *res = -1; + return false; + } + else { + *res = (a[j].bits > b[j].bits) ? 1 : -1; + return false; + } + } + } + *res = 0; + return false; +} + +static Py_ssize_t * +mergesort_vectors(PyObject *arg) +{ + Py_ssize_t R = PyList_GET_SIZE(arg); + Py_ssize_t N = Py_SIZE(PyList_GET_ITEM(arg, 0)); + Py_ssize_t *arr1 = PyMem_Malloc(R*sizeof(Py_ssize_t)); + Py_ssize_t *arr2 = PyMem_Malloc(R*sizeof(Py_ssize_t)); + if (arr1 == NULL || arr2 == NULL) { + PyErr_NoMemory(); + return NULL; + } + for (Py_ssize_t i = 0; i < R; i++) { + arr1[i] = i; + } + // bottom-up mergesort + for (Py_ssize_t width = 1; width < R; width *= 2) { + for (Py_ssize_t i0 = 0; i0 < R; i0 += 2*width) { + Py_ssize_t i1 = i0 + width; + Py_ssize_t i2 = i1 + width; + if (i1 > R) { + i1 = R; + } + if (i2 > R) { + i2 = R; + } + // arr2[i0:i2] = merge(arr1[i0:i1], arr1[i1:i2]) + Py_ssize_t left = i0; + Py_ssize_t right = i1; + for (Py_ssize_t k = i0; k < i2; k++) { + assert(left >= i0); + assert(left <= i1); + assert(right >= i1); + assert(right <= i2); + if (left == i1) { + arr2[k] = arr1[right++]; + continue; + } + if (right == i2) { + assert(left < i1); + arr2[k] = arr1[left++]; + continue; + } + int cmp; + if (Vector_cmp( + Vector_get_vec(PyList_GET_ITEM(arg, arr1[left])), + Vector_get_vec(PyList_GET_ITEM(arg, arr1[right])), + N, &cmp + )) { + goto error; + } + if (cmp <= 0) { + arr2[k] = arr1[left++]; + } + else { + arr2[k] = arr1[right++]; + } + } + } + { + Py_ssize_t *tmp = arr2; + arr2 = arr1; + arr1 = tmp; + } + } + PyMem_Free(arr2); + return arr1; +error: + PyMem_Free(arr1); + PyMem_Free(arr2); + return NULL; +} + +static bool +cross_off_rows_and_columns(PyObject *arg, char *col_active, char *row_active) +{ + Py_ssize_t R = PyList_GET_SIZE(arg); + Py_ssize_t N = Py_SIZE(PyList_GET_ITEM(arg, 0)); + Py_ssize_t *stack = PyMem_Malloc(2*N*sizeof(Py_ssize_t)); + Py_ssize_t *column_to_count = PyMem_Calloc(N, sizeof(Py_ssize_t)); + if (stack == NULL || column_to_count == NULL) { + PyMem_Free(stack); + PyMem_Free(column_to_count); + return true; + } + for (Py_ssize_t i = 0; i < R; i++) { + if (!row_active[i]) { + continue; + } + TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); + for (Py_ssize_t j = 0; j < N; j++) { + if (col_active[j] && !TagInt_is_zero(vec[j])) { + column_to_count[j]++; + } + } + } + Py_ssize_t stacktop = 0; + for (Py_ssize_t j0 = 0; j0 < N; j0++) { + // Maintain a stack of column indexes that have reached + // either 0 or 1 nonzero entries. + stack[stacktop++] = j0; + while (stacktop) { + Py_ssize_t j = stack[--stacktop]; + if (!col_active[j] || column_to_count[j] >= 2) { + continue; + } + if (column_to_count[j] == 0) { + col_active[j] = 0; + continue; + } + assert(column_to_count[j] == 1); + Py_ssize_t i = 0; + while (!row_active[i] || TagInt_is_zero(Vector_get_vec(PyList_GET_ITEM(arg, i))[j])) { + i++; + assert(i < R); + } + row_active[i] = 0; + col_active[j] = 0; + TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); + for (Py_ssize_t j1 = 0; j1 < N; j1++) { + if (!col_active[j1] || TagInt_is_zero(vec[j1])) { + continue; + } + column_to_count[j1]--; + assert(column_to_count[j1] >= 0); + if (column_to_count[j1] <= 1) { + stack[stacktop++] = j1; + // The stack can't overflow because each column is added + // at most twice: once when column_to_count[j] becomes 1, + // then perhaps again when it becomes 0. + assert(stacktop <= 2*N); + } + } + } + } + PyMem_Free(stack); + PyMem_Free(column_to_count); + return false; +} + +static bool +partition_active_columns( + PyObject *arg, char *col_active, char *row_active, + Py_ssize_t *p_num_components, + Py_ssize_t **p_component_to_size, + Py_ssize_t **p_cols_by_component, + Py_ssize_t ***p_components +) +{ + Py_ssize_t R = PyList_GET_SIZE(arg); + Py_ssize_t N = Py_SIZE(PyList_GET_ITEM(arg, 0)); + Py_ssize_t *uf_parent = PyMem_Malloc(N*sizeof(Py_ssize_t)); + Py_ssize_t *uf_size = PyMem_Malloc(N*sizeof(Py_ssize_t)); + Py_ssize_t *root_to_component_index = PyMem_Malloc(N*sizeof(Py_ssize_t)); + Py_ssize_t *component_to_size = NULL; + Py_ssize_t *component_to_index = NULL; + Py_ssize_t *cols_by_component = NULL; + Py_ssize_t **components = NULL; + if (!(uf_parent && uf_size && root_to_component_index)) { + PyErr_NoMemory(); + goto error; + } + Py_ssize_t num_cols_active = 0; + for (Py_ssize_t j = 0; j < N; j++) { + if (col_active[j]) { + uf_parent[j] = j; + uf_size[j] = 1; + num_cols_active++; + } + } + for (Py_ssize_t i = 0; i < R; i++) { + if (!row_active[i]) { + continue; + } + TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); + Py_ssize_t j0 = 0; + while (!col_active[j0] || TagInt_is_zero(vec[j0])) { + j0++; + } + assert(j0 < N && col_active[j0] && !TagInt_is_zero(vec[j0])); + for (Py_ssize_t j1 = j0 + 1; j1 < N; j1++) { + if (col_active[j1] && !TagInt_is_zero(vec[j1])) { + uf_union(uf_parent, uf_size, j0, j1); + } + } + } + Py_ssize_t num_components = 0; + for (Py_ssize_t j = 0; j < N; j++) { + if (col_active[j]) { + Py_ssize_t fj = uf_find(uf_parent, j); + uf_parent[j] = fj; + if (j == fj) { + root_to_component_index[j] = num_components++; + } + } + } + component_to_size = PyMem_Calloc(num_components, sizeof(Py_ssize_t)); + component_to_index = PyMem_Malloc(num_components * sizeof(Py_ssize_t)); + cols_by_component = PyMem_Malloc(num_cols_active * sizeof(Py_ssize_t)); + components = PyMem_Malloc(num_components * sizeof(Py_ssize_t *)); + if (!(component_to_size && cols_by_component && components)) { + PyErr_NoMemory(); + goto error; + } + for (Py_ssize_t j = 0; j < N; j++) { + if (col_active[j] && uf_parent[j] == j) { + component_to_size[root_to_component_index[j]] = uf_size[j]; + } + } + // Sort columns into the components, similar to a counting sort. + Py_ssize_t sum = 0; + for (Py_ssize_t k = 0; k < num_components; k++) { + components[k] = &cols_by_component[sum]; // points to the start of each slice + sum += component_to_size[k]; + component_to_index[k] = sum; // indexes of the end of each slice + } + for (Py_ssize_t j = N - 1; j >= 0; j--) { + if (col_active[j]) { + cols_by_component[--component_to_index[root_to_component_index[uf_parent[j]]]] = j; + } + } + for (Py_ssize_t k = 0; k < num_components; k++) { + assert(components[k] == &cols_by_component[component_to_index[k]]); + } + PyMem_Free(uf_parent); + PyMem_Free(uf_size); + PyMem_Free(root_to_component_index); + PyMem_Free(component_to_index); + *p_num_components = num_components; + *p_component_to_size = component_to_size; + *p_cols_by_component = cols_by_component; + *p_components = components; + return false; +error: + PyMem_Free(uf_parent); + PyMem_Free(uf_size); + PyMem_Free(root_to_component_index); + PyMem_Free(component_to_size); + PyMem_Free(component_to_index); + PyMem_Free(cols_by_component); + PyMem_Free(components); + return true; +} + +static bool +append_subproblem( + PyObject *arg, + Py_ssize_t *cols, Py_ssize_t num_cols, + Py_ssize_t *rows, Py_ssize_t num_rows, + PyObject *subproblem_rows_list, + PyObject *subproblems_list +) +{ + PyObject *subproblem_rows = Vector_zero_impl(num_rows); + if (subproblem_rows == NULL) { + return true; + } + TagInt *rows_vec = Vector_get_vec(subproblem_rows); + for (Py_ssize_t ii = 0; ii < num_rows; ii++) { + static_assert(sizeof(intptr_t) == sizeof(Py_ssize_t)); + rows_vec[ii] = pack_integer((intptr_t)rows[ii]); + } + if (PyList_Append(subproblem_rows_list, subproblem_rows) < 0) { + Py_DECREF(subproblem_rows); + return true; + } + Py_DECREF(subproblem_rows); + + PyObject *subproblem = PyList_New(num_rows); + if (subproblem == NULL) { + return true; + } + for (Py_ssize_t ii = 0; ii < num_rows; ii++) { + PyObject *v = Vector_zero_impl(num_cols); + if (v == NULL) { + Py_DECREF(subproblem); + return true; + } + TagInt *dest = Vector_get_vec(v); + TagInt *source = Vector_get_vec(PyList_GET_ITEM(arg, rows[ii])); + for (Py_ssize_t jj = 0; jj < num_cols; jj++) { + dest[jj] = TagInt_copy(source[cols[jj]]); + } + PyList_SET_ITEM(subproblem, ii, v); + } + if (PyList_Append(subproblems_list, subproblem) < 0) { + Py_DECREF(subproblem); + return true; + } + Py_DECREF(subproblem); + return false; +} + +static PyObject * +decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) +{ + if (!PyList_CheckExact(arg)) { + PyErr_SetString(PyExc_TypeError, "decompose_relations_among(vecs) argument must be a list"); + return NULL; + } + Py_ssize_t R = PyList_GET_SIZE(arg); + if (R == 0) { + return Py_BuildValue("([][][])"); + } + Py_ssize_t N = -1; + for (Py_ssize_t i = 0; i < R; i++) { + PyObject *v = PyList_GET_ITEM(arg, i); + if (Py_TYPE(v) != &Vector_Type) { + PyErr_SetString(PyExc_TypeError, "decompose_relations_among(vecs) argument must be a list of Vectors"); + return NULL; + } + if (i == 0) { + N = Py_SIZE(v); + } else { + if (Py_SIZE(v) != N) { + PyErr_SetString(PyExc_ValueError, "length mismatch in decompose_relations_among"); + return NULL; + } + } + } + if (!is_packable_int(R - 1)) { + PyErr_SetString(PyExc_OverflowError, "too many rows to partition"); + return NULL; + } + PyObject *easy_relations_list = NULL; + PyObject *subproblem_rows_list = NULL; + PyObject *subproblems_list = NULL; + PyObject *result = NULL; + char *row_active = NULL; + char *col_active = NULL; + Py_ssize_t *sorted_rows = NULL; + Py_ssize_t *component_to_size = NULL; + Py_ssize_t *cols_by_component = NULL; + Py_ssize_t **components = NULL; + Py_ssize_t *rows = NULL; + + if (!(easy_relations_list = PyList_New(0))) { + goto error; + } + if (!(subproblem_rows_list = PyList_New(0))) { + goto error; + } + if (!(subproblems_list = PyList_New(0))) { + goto error; + } + if (!(row_active = PyMem_Malloc(R)) || + !(col_active = PyMem_Malloc(N)) || + !(rows = PyMem_Malloc(R * sizeof(Py_ssize_t))) + ) { + PyErr_NoMemory(); + return NULL; + } + memset(row_active, 1, R); + memset(col_active, 1, N); + // Step 1: look for zero rows + for (Py_ssize_t i = 0; i < R; i++) { + TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); + Py_ssize_t j = 0; + while (j < N && TagInt_is_zero(vec[j])) { + j++; + } + if (j < N) { + continue; + } + row_active[i] = 0; + PyObject *v = Vector_zero_impl(R); + if (v == NULL) { + goto error; + } + Vector_get_vec(v)[i] = TagInt_ONE; + if (PyList_Append(easy_relations_list, v) < 0) { + Py_DECREF(v); + goto error; + } + Py_DECREF(v); + } + // Step 2: Look for duplicate vectors + if (!(sorted_rows = mergesort_vectors(arg))) { + goto error; + } + for (Py_ssize_t ii0 = 0; ii0 < R;) { + Py_ssize_t ii1 = ii0 + 1; + for (; ii1 < R; ii1++) { + // extend the streak sorted_rows[ii0:ii1] + assert(ii0 < R && ii1 < R); + int cmp; + if (Vector_cmp( + Vector_get_vec(PyList_GET_ITEM(arg, sorted_rows[ii0])), + Vector_get_vec(PyList_GET_ITEM(arg, sorted_rows[ii1])), + N, &cmp + )) { + goto error; + } + assert(cmp <= 0); + if (cmp != 0) { + break; + } + } + // we already handled the zero rows + if (row_active[sorted_rows[ii0]]) { + for (Py_ssize_t ii = ii0; ii < ii1 - 1; ii++) { + row_active[sorted_rows[ii]] = 0; + PyObject *v = Vector_zero_impl(R); + if (v == NULL) { + goto error; + } + Vector_get_vec(v)[sorted_rows[ii]] = TagInt_ONE; + Vector_get_vec(v)[sorted_rows[ii1 - 1]] = TagInt_NEGATIVE_ONE; + if (PyList_Append(easy_relations_list, v) < 0) { + Py_DECREF(v); + goto error; + } + Py_DECREF(v); + } + } + ii0 = ii1; + } + // Step 3: cross off rows and columns with <= 1 nonzero entry. + if (cross_off_rows_and_columns(arg, col_active, row_active)) { + goto error; + } + // Step 4: Find which columns have nonzero entries in the same row. + Py_ssize_t num_components; + if (partition_active_columns( + arg, col_active, row_active, + &num_components, &component_to_size, &cols_by_component, &components + )) { + goto error; + } + // Step 5: Make each subproblem of the decomposition + for (Py_ssize_t k = 0; k < num_components; k++) { + Py_ssize_t num_cols = component_to_size[k]; + Py_ssize_t *cols = components[k]; + Py_ssize_t num_rows = 0; + for (Py_ssize_t i = 0; i < R; i++) { + // find the rows associated to this component. + if (!row_active[i]) { + continue; + } + TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); + Py_ssize_t jj = 0; + for (; jj < num_cols; jj++) { + if (!TagInt_is_zero(vec[cols[jj]])) { + break; + } + } + if (jj < num_cols) { + // broke out early: found a nonzero entry + rows[num_rows++] = i; + } + } + if (append_subproblem( + arg, cols, num_cols, rows, num_rows, + subproblem_rows_list, subproblems_list + )) { + goto error; + } + } + result = PyTuple_Pack(3, easy_relations_list, + subproblem_rows_list, + subproblems_list); +error: + Py_XDECREF(easy_relations_list); + Py_XDECREF(subproblem_rows_list); + Py_XDECREF(subproblems_list); + PyMem_Free(row_active); + PyMem_Free(col_active); + PyMem_Free(sorted_rows); + PyMem_Free(component_to_size); + PyMem_Free(cols_by_component); + PyMem_Free(components); + PyMem_Free(rows); + return result; +} + /*********************************************************************/ /* Module stuff */ /*********************************************************************/ @@ -4232,6 +4751,9 @@ static PyMethodDef mutable_lattice_methods[] = { "xgcd(a, b) returns a triple (x, y, g) of integers with x*a + y*b == g == gcd(a, b)"}, {"relations_among", relations_among, METH_O, "relations_among([v0, ..., vk]) returns the Lattice of coefficient vectors for linear dependencies among the given Vectors"}, + {"decompose_relations_among", decompose_relations_among, METH_O, + "decompose_relations_among([v0, ..., vk]) reduces a relations_among calculation to sub-problems. " + "It returns a 3-tuple (easy_relations: list[Vector], subproblem_rows : list[Vector], subproblems: list[list[Vector]])"}, {"transpose", (PyCFunction)(void(*)(void))transpose, METH_FASTCALL, "transpose(N, [v0, ..., vk]) transposes a length-k list of legth-N vectors into a length-N list of length-k vectors"}, {NULL, NULL, 0, NULL} /* sentinel */ diff --git a/mutable_lattice/__init__.py b/mutable_lattice/__init__.py index 39ee170..d68f10e 100644 --- a/mutable_lattice/__init__.py +++ b/mutable_lattice/__init__.py @@ -4,6 +4,17 @@ generalized_row_op, Lattice, xgcd, - relations_among, + relations_among as relations_among_c, + decompose_relations_among, transpose, -) \ No newline at end of file +) + +def relations_among_with_decomposition(vecs, /): + R = len(vecs) + relations, subproblem_rows, subproblems = decompose_relations_among(vecs) + for rows, subproblem in zip(subproblem_rows, subproblems, strict=True): + for rel in relations_among_c(subproblem).get_basis(): + relations.append(rel.shuffled_by_action(rows, R)) + return Lattice(R, relations, maxrank=len(relations)) + +relations_among = relations_among_with_decomposition diff --git a/mutable_lattice/test.py b/mutable_lattice/test.py index 4f1af42..83b7c2f 100644 --- a/mutable_lattice/test.py +++ b/mutable_lattice/test.py @@ -10,7 +10,9 @@ generalized_row_op, Lattice, xgcd, - relations_among, + relations_among_c, + relations_among_with_decomposition, + decompose_relations_among, transpose, ) from .pylattice import PyLattice @@ -1316,7 +1318,10 @@ def test_lattice_compare_error(self): with self.assertRaisesRegex(ValueError, "different ambient dimensions"): L1 < L2 -class TestKernels(unittest.TestCase): +class TestKernels: + def ker(self, vecs): + raise NotImplementedError + def setUp(self): self.VALUES = make_some_values() @@ -1324,38 +1329,69 @@ def tearDown(self): self.VALUES.clear() del self.VALUES - def test_relations_among(self): + def test_relations_among_basic(self): self.assertEqual( - relations_among([Vector([10]*100), Vector([20]*100)]), - Lattice(2, [[2, -1]]) - ) - self.assertEqual( - relations_among([Vector([60]), Vector([100]), Vector([150])]), - Lattice(3, [[5, -3, 0], [0, 3, -2]]) + self.ker([]), + Lattice(0), ) self.assertEqual( - relations_among([]), - Lattice(0), + self.ker([Vector([17])]), + Lattice(1), ) self.assertEqual( - relations_among([Vector([]), Vector([]), Vector([])]), + self.ker([Vector([]), Vector([]), Vector([])]), Lattice.full(3), ), self.assertEqual( - relations_among([Vector([17])]), - Lattice(1), + self.ker([Vector([17]), Vector([0])]), + Lattice(2, [[0, 1]]), ) self.assertEqual( - relations_among([Vector([17]), Vector([0])]), - Lattice(2, [[0, 1]]), + self.ker([Vector([10]), Vector([20])]), + Lattice(2, [[2, -1]]) + ) + self.assertEqual( + self.ker([Vector([10, 100]), Vector([20, 200])]), + Lattice(2, [[2, -1]]) + ) + self.assertEqual( + self.ker([Vector([10]*100), Vector([20]*100)]), + Lattice(2, [[2, -1]]) ) self.assertEqual( - relations_among([Vector([0, 0, 0])]*100), + self.ker([Vector([60]), Vector([100]), Vector([150])]), + Lattice(3, [[5, -3, 0], [0, 3, -2]]) + ) + self.assertEqual( + self.ker([Vector([0, 0, 0])]*100), Lattice.full(100) ) + self.assertEqual( + self.ker([Vector([-1, -1]), Vector([-2, 0])]), + Lattice(2) + ) + self.assertEqual( + self.ker([Vector([1, 1]), Vector([0, 1])]), + Lattice(2) + ) + self.assertEqual( + self.ker([ + Vector([0,0,0,5,5,5,0]), + Vector([0,1,2,0,0,0,0]), + Vector([0,2,4,0,0,0,0]), + Vector([0,0,0,0,0,0,0]), + Vector([0,0,0,0,3,3,0]), + Vector([0,0,0,0,1,1,0]), + Vector([0,0,0,0,0,0,0]), + ]), + Lattice(7, [[0, 2,-1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1,-3, 0], + [0, 0, 0, 0, 0, 0, 1]]) + ) # From https://github.com/sagemath/sage/blob/develop/src/sage/matrix/matrix_integer_dense.pyx self.assertEqual( - relations_among([ + self.ker([ Vector([4, 1, 0, 4]), Vector([7, 0, 1, 7]), Vector([9, 5, 0, 6]), @@ -1366,6 +1402,7 @@ def test_relations_among(self): Lattice(6, [[26, -31, 30, -21, -2, 10], [-47, -13, 48, -14, -11, 18]]), ) + def test_relations_among_random(self): for values in [ [-2, -1, 0, 0, 0, 0, 1, 2], @@ -1374,26 +1411,139 @@ def test_relations_among_random(self): for N in range(5): zero_N = Vector.zero(N) for R in range(5): - for _ in range(10): + for _ in range(100): vecs = [Vector(random.choices(values, k=N)) for _ in range(R)] rank = Lattice(N, vecs, maxrank=R).rank - kernel = relations_among(vecs) + kernel = self.ker(vecs) # Ensure the kernel has the correct rank self.assertEqual(kernel.rank + rank, R) # Ensure the computed kernel vectors are actually relations for row in kernel.tolist(): lin_combo = sum((c*vec for (c, vec) in zip(row, vecs)), start=zero_N) - self.assertEqual(lin_combo, zero_N) + self.assertEqual(lin_combo, zero_N, vecs) # Ensure the kernel is saturated self.assertEqual(kernel.invariants(), [1]*kernel.rank + [0] * rank) def test_relations_among_errors(self): with self.assertRaisesRegex(TypeError, "must be a list"): - relations_among((1, 2, 3)) + self.ker((1, 2, 3)) with self.assertRaisesRegex(TypeError, "must be a list of Vectors"): - relations_among([1, 2, 3]) + self.ker([1, 2, 3]) with self.assertRaisesRegex(ValueError, "length mismatch"): - relations_among([Vector([1,2,3]), Vector([4,5,6,7])]) + self.ker([Vector([1,2,3]), Vector([4,5,6,7])]) + +class TestKernelsC(TestKernels, unittest.TestCase): + def ker(self, vecs): + return relations_among_c(vecs) + +class TestKernelsWithDecomposition(TestKernels, unittest.TestCase): + def ker(self, vecs): + return relations_among_with_decomposition(vecs) + +class TestDecomposeRelationsAmong(unittest.TestCase): + def test_decompose_relations_among_basic(self): + self.assertEqual( + # find zero rows + decompose_relations_among([ + Vector([]), + Vector([]), + Vector([]), + ]), + ([Vector([1,0,0]), Vector([0,1,0]), Vector([0, 0, 1])], + [], + []) + ) + self.assertEqual( + decompose_relations_among([ + Vector([0,0]), + Vector([0,0]), + Vector([0,0]), + ]), + ([Vector([1,0,0]), Vector([0,1,0]), Vector([0, 0, 1])], + [], + []) + ) + # deduplicate + self.assertEqual( + decompose_relations_among([ + Vector([17]), + Vector([17]), + Vector([17]), + ]), + ([Vector([1,0,-1]), Vector([0,1,-1])], + [], + []) + ) + # split into smaller problems + self.assertEqual( + decompose_relations_among([ + Vector([0, 0, 0]), + Vector([0, 5, 0]), + Vector([0, 6, 0]), + Vector([0, 0, 0]), + ]), + ([Vector([1,0,0,0]), Vector([0,0,0,1])], + [Vector([1, 2])], [[Vector([5]), Vector([6])]]) + ) + self.assertEqual( + decompose_relations_among([ + Vector([0,0,1,1,0,0,0]), + Vector([0,0,2,2,0,0,0]), + Vector([0,0,0,0,0,3,3]), + Vector([0,0,0,0,0,4,4]), + ]), + ([], + [Vector([0, 1]), Vector([2, 3])], + [[Vector([1,1]),Vector([2,2])], [Vector([3,3]),Vector([4,4])]]) + ) + # All together now + self.assertEqual( + decompose_relations_among([ + Vector([0,0,0,5,5,5,0]), + Vector([0,1,2,0,0,0,0]), + Vector([0,2,4,0,0,0,0]), + Vector([0,0,0,0,0,0,0]), + Vector([0,0,0,0,3,3,0]), + Vector([0,0,0,0,1,1,0]), + Vector([0,0,0,0,0,0,0]), + Vector([0,0,0,5,5,5,0]), + ]), + ([Vector([0, 0, 0, 1, 0, 0, 0, 0]), + Vector([0, 0, 0, 0, 0, 0, 1, 0]), + Vector([1, 0, 0, 0, 0, 0, 0, -1])], + [Vector([1, 2]), Vector([4, 5])], + [[Vector([1,2]),Vector([2,4])], [Vector([3,3]),Vector([1,1])]]) + ) + # bigints are okay + vecs = [Vector([-36893488147419103234]), Vector([9223372036854775808])] + self.assertEqual( + decompose_relations_among(vecs), + ([], [Vector([0, 1])], [vecs]) + ) + # Be sure that crossing out rows/columns can let other rows/columns be + # crossed out. + self.assertEqual( + decompose_relations_among([ + Vector([1,2,0,0,0,0]), + Vector([0,3,4,0,0,0]), + Vector([0,0,5,6,0,0]), + ]), + ([], [], []) + ) + self.assertEqual( + decompose_relations_among([Vector([0]*i + [1]*(100-i)) for i in range(100)]), + ([], [], []) + ) + + + +class TestTranspose(unittest.TestCase): + def setUp(self): + self.VALUES = make_some_values() + + def tearDown(self): + self.VALUES.clear() + del self.VALUES def test_transpose(self): self.assertEqual(transpose(2, [Vector([1, 2]), Vector([3, 4])]), [Vector([1, 3]), Vector([2, 4])]) diff --git a/pyproject.toml b/pyproject.toml index c1a52fd..d8bf654 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ test-command = "python -m mutable_lattice.test -v" [project] name = "mutable_lattice" -version = "0.3.0" +version = "0.4.0" authors = [ { name="Dennis Sweeney", email="sweeney.427@osu.edu" }, ] From 0b805125b7561f2afb47084d777abf3ca2eaca8a Mon Sep 17 00:00:00 2001 From: Dennis Sweeney Date: Tue, 10 Mar 2026 16:18:33 -0400 Subject: [PATCH 2/4] update pyi --- mutable_lattice/__init__.pyi | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/mutable_lattice/__init__.pyi b/mutable_lattice/__init__.pyi index 27c6734..d42fbd8 100644 --- a/mutable_lattice/__init__.pyi +++ b/mutable_lattice/__init__.pyi @@ -59,3 +59,11 @@ def transpose(N: int, vectors: list[Vector], /) -> list[Vector]: def relations_among(vectors: list[Vector], /) -> Lattice: ... +def relations_among_c(vectors: list[Vector], /) -> Lattice: + ... + +def relations_among_with_decomposition(vectors: list[Vector], /) -> Lattice: + ... + +def decompose_relations_among(vectors: list[Vector], /) -> tuple[list[Vector], list[Vector], list[list[Vector]]]: + ... From 802ed8948478f14fa5df8dc9047cc08826c12f58 Mon Sep 17 00:00:00 2001 From: Dennis Sweeney Date: Tue, 10 Mar 2026 21:34:15 -0400 Subject: [PATCH 3/4] more tests --- mutable_lattice/test.py | 64 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/mutable_lattice/test.py b/mutable_lattice/test.py index 83b7c2f..2aa073a 100644 --- a/mutable_lattice/test.py +++ b/mutable_lattice/test.py @@ -1442,6 +1442,17 @@ def ker(self, vecs): class TestDecomposeRelationsAmong(unittest.TestCase): def test_decompose_relations_among_basic(self): + self.assertEqual(decompose_relations_among([]), ([], [], [])) + self.assertEqual(decompose_relations_among([Vector([0])]), + ([Vector([1])], [], [])) + self.assertEqual(decompose_relations_among([Vector([17])]), + ([], [], [])) + self.assertEqual(decompose_relations_among([Vector([0, 17, 0])]), + ([], [], [])) + self.assertEqual(decompose_relations_among([Vector([0, 0, 0])]), + ([Vector([1])], [], [])) + self.assertEqual(decompose_relations_among([Vector([1,2,3])]), + ([], [], [])) self.assertEqual( # find zero rows decompose_relations_among([ @@ -1463,6 +1474,15 @@ def test_decompose_relations_among_basic(self): [], []) ) + self.assertEqual( + decompose_relations_among([ + Vector([1,1]), + Vector([0,0]), + Vector([2,2]), + ]), + ([Vector([0,1,0])], [Vector([0,2])], [[Vector([1,1]), + Vector([2,2])]]) + ) # deduplicate self.assertEqual( decompose_relations_among([ @@ -1474,6 +1494,29 @@ def test_decompose_relations_among_basic(self): [], []) ) + self.assertEqual( + decompose_relations_among([ + Vector([10**k for k in range(100)]), + Vector([10**k for k in range(100)]), + Vector([10**k for k in range(100)]), + ]), + ([Vector([1,0,-1]), Vector([0,1,-1])], + [], + []) + ) + self.assertEqual( + decompose_relations_among([ + Vector([1,1]), + Vector([2,2]), + Vector([2,2]), + Vector([3,3]), + Vector([3,3]), + Vector([3,3]), + ]), + ([Vector([0,1,-1,0,0,0]),Vector([0,0,0,1,0,-1]),Vector([0,0,0,0,1,-1])], + [Vector([0, 2, 5])], + [[Vector([1,1]), Vector([2,2]), Vector([3,3])]]) + ) # split into smaller problems self.assertEqual( decompose_relations_among([ @@ -1496,6 +1539,17 @@ def test_decompose_relations_among_basic(self): [Vector([0, 1]), Vector([2, 3])], [[Vector([1,1]),Vector([2,2])], [Vector([3,3]),Vector([4,4])]]) ) + self.assertEqual( + decompose_relations_among([ + Vector([1,0,1,0]), + Vector([2,0,2,0]), + Vector([0,3,0,3]), + Vector([0,4,0,4]), + ]), + ([], + [Vector([0,1]), Vector([2,3])], + [[Vector([1,1]),Vector([2,2])], [Vector([3,3]),Vector([4,4])]]) + ) # All together now self.assertEqual( decompose_relations_among([ @@ -1534,7 +1588,15 @@ def test_decompose_relations_among_basic(self): decompose_relations_among([Vector([0]*i + [1]*(100-i)) for i in range(100)]), ([], [], []) ) - + self.assertEqual( + decompose_relations_among([ + Vector([1,2,0,0]), + Vector([1,2,0,0]), + Vector([3,4,5,6]), + Vector([3,4,5,6]), + ]), + ([Vector([1, -1, 0, 0]), Vector([0, 0, 1, -1])], [], []) + ) class TestTranspose(unittest.TestCase): From 6022e9f8afff5b94b7849054128c2fd25d4d42da Mon Sep 17 00:00:00 2001 From: Dennis Sweeney Date: Thu, 12 Mar 2026 02:52:17 -0400 Subject: [PATCH 4/4] Refactor. Do a BFS instead of using a union-find. --- _mutable_lattice.c | 251 +++++++++++++++++------------------- mutable_lattice/__init__.py | 2 +- mutable_lattice/test.py | 48 ++++++- 3 files changed, 164 insertions(+), 137 deletions(-) diff --git a/_mutable_lattice.c b/_mutable_lattice.c index c0b0159..ea31f0f 100644 --- a/_mutable_lattice.c +++ b/_mutable_lattice.c @@ -1763,7 +1763,6 @@ bisect_left(Py_ssize_t *a, Py_ssize_t lo, Py_ssize_t hi, Py_ssize_t x) static PyObject * Lattice__assert_consistent(PyObject *self, PyObject *Py_UNUSED(other)) { - // printf("asserting consistent\n"); Lattice *L = (Lattice *)self; Py_ssize_t N = L->N; Py_ssize_t R = L->rank; @@ -4191,11 +4190,13 @@ transpose(PyObject *Py_UNUSED(mod), PyObject *const *args, Py_ssize_t nargs) } static bool -Vector_cmp(TagInt *a, TagInt *b, Py_ssize_t N, int *res) +Vector_cmp(PyObject *arg, Py_ssize_t i1, Py_ssize_t i2, Py_ssize_t N, int *res) { + TagInt *a = Vector_get_vec(PyList_GET_ITEM(arg, i1)); + TagInt *b = Vector_get_vec(PyList_GET_ITEM(arg, i2)); // Compare two vectors of equal length. // This isn't a useful order--it just has to be a total order - // so that we can efficiently deduplicate + // so that we can efficiently deduplicate. for (Py_ssize_t j = 0; j < N; j++) { if (a[j].bits == b[j].bits) { continue; @@ -4281,11 +4282,7 @@ mergesort_vectors(PyObject *arg) continue; } int cmp; - if (Vector_cmp( - Vector_get_vec(PyList_GET_ITEM(arg, arr1[left])), - Vector_get_vec(PyList_GET_ITEM(arg, arr1[right])), - N, &cmp - )) { + if (Vector_cmp(arg, arr1[left], arr1[right], N, &cmp)) { goto error; } if (cmp <= 0) { @@ -4353,6 +4350,7 @@ cross_off_rows_and_columns(PyObject *arg, char *col_active, char *row_active) i++; assert(i < R); } + assert(row_active[i] && col_active[j]); row_active[i] = 0; col_active[j] = 0; TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); @@ -4381,104 +4379,109 @@ static bool partition_active_columns( PyObject *arg, char *col_active, char *row_active, Py_ssize_t *p_num_components, - Py_ssize_t **p_component_to_size, - Py_ssize_t **p_cols_by_component, - Py_ssize_t ***p_components + Py_ssize_t **p_component_to_num_rows, + Py_ssize_t **p_component_to_num_cols, + Py_ssize_t **p_rows_by_component, + Py_ssize_t **p_cols_by_component ) { Py_ssize_t R = PyList_GET_SIZE(arg); Py_ssize_t N = Py_SIZE(PyList_GET_ITEM(arg, 0)); - Py_ssize_t *uf_parent = PyMem_Malloc(N*sizeof(Py_ssize_t)); - Py_ssize_t *uf_size = PyMem_Malloc(N*sizeof(Py_ssize_t)); - Py_ssize_t *root_to_component_index = PyMem_Malloc(N*sizeof(Py_ssize_t)); - Py_ssize_t *component_to_size = NULL; - Py_ssize_t *component_to_index = NULL; - Py_ssize_t *cols_by_component = NULL; - Py_ssize_t **components = NULL; - if (!(uf_parent && uf_size && root_to_component_index)) { - PyErr_NoMemory(); - goto error; - } - Py_ssize_t num_cols_active = 0; - for (Py_ssize_t j = 0; j < N; j++) { - if (col_active[j]) { - uf_parent[j] = j; - uf_size[j] = 1; - num_cols_active++; - } - } + Py_ssize_t num_active_rows = 0; for (Py_ssize_t i = 0; i < R; i++) { - if (!row_active[i]) { - continue; - } - TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); - Py_ssize_t j0 = 0; - while (!col_active[j0] || TagInt_is_zero(vec[j0])) { - j0++; - } - assert(j0 < N && col_active[j0] && !TagInt_is_zero(vec[j0])); - for (Py_ssize_t j1 = j0 + 1; j1 < N; j1++) { - if (col_active[j1] && !TagInt_is_zero(vec[j1])) { - uf_union(uf_parent, uf_size, j0, j1); - } + if (row_active[i]) { + num_active_rows++; } } - Py_ssize_t num_components = 0; + Py_ssize_t num_active_cols = 0; for (Py_ssize_t j = 0; j < N; j++) { if (col_active[j]) { - Py_ssize_t fj = uf_find(uf_parent, j); - uf_parent[j] = fj; - if (j == fj) { - root_to_component_index[j] = num_components++; - } - } - } - component_to_size = PyMem_Calloc(num_components, sizeof(Py_ssize_t)); - component_to_index = PyMem_Malloc(num_components * sizeof(Py_ssize_t)); - cols_by_component = PyMem_Malloc(num_cols_active * sizeof(Py_ssize_t)); - components = PyMem_Malloc(num_components * sizeof(Py_ssize_t *)); - if (!(component_to_size && cols_by_component && components)) { + num_active_cols++; + } + } + Py_ssize_t max_num_components = Py_MIN(num_active_rows, num_active_cols); + Py_ssize_t *rowstack = PyMem_Malloc(num_active_rows * sizeof(Py_ssize_t)); + Py_ssize_t *colstack = PyMem_Malloc(num_active_cols * sizeof(Py_ssize_t)); + Py_ssize_t *row_to_component = PyMem_Malloc(R * sizeof(Py_ssize_t)); + Py_ssize_t *col_to_component = PyMem_Malloc(N * sizeof(Py_ssize_t)); + Py_ssize_t *component_to_num_rows = PyMem_Calloc(max_num_components, sizeof(Py_ssize_t)); + Py_ssize_t *component_to_num_cols = PyMem_Calloc(max_num_components, sizeof(Py_ssize_t)); + Py_ssize_t *rows_by_component = PyMem_Malloc(num_active_rows * sizeof(Py_ssize_t)); + Py_ssize_t *cols_by_component = PyMem_Malloc(num_active_cols * sizeof(Py_ssize_t)); + Py_ssize_t num_components = 0; + Py_ssize_t rowstacktop = 0, colstacktop = 0; + Py_ssize_t num_rows_done = 0, num_cols_done = 0; + if (!(rowstack && colstack && + row_to_component && col_to_component && + component_to_num_rows && component_to_num_cols && + rows_by_component && cols_by_component + )) { PyErr_NoMemory(); - goto error; + PyMem_Free(rowstack); PyMem_Free(colstack); + PyMem_Free(row_to_component); PyMem_Free(col_to_component); + PyMem_Free(component_to_num_rows); PyMem_Free(component_to_num_cols); + PyMem_Free(rows_by_component); PyMem_Free(cols_by_component); + return true; } - for (Py_ssize_t j = 0; j < N; j++) { - if (col_active[j] && uf_parent[j] == j) { - component_to_size[root_to_component_index[j]] = uf_size[j]; - } + for (Py_ssize_t i = 0; i < R; i++) { + row_to_component[i] = -1; } - // Sort columns into the components, similar to a counting sort. - Py_ssize_t sum = 0; - for (Py_ssize_t k = 0; k < num_components; k++) { - components[k] = &cols_by_component[sum]; // points to the start of each slice - sum += component_to_size[k]; - component_to_index[k] = sum; // indexes of the end of each slice + for (Py_ssize_t j = 0; j < N; j++) { + col_to_component[j] = -1; } - for (Py_ssize_t j = N - 1; j >= 0; j--) { - if (col_active[j]) { - cols_by_component[--component_to_index[root_to_component_index[uf_parent[j]]]] = j; + // Explore the bipartite graph of rows and columns + // connected whenever they intersect at a nonzero entry. + for (Py_ssize_t i0 = 0; i0 < R; i0++) { + if (!row_active[i0] || row_to_component[i0] >= 0) { + continue; + } + Py_ssize_t k = num_components++; + assert(rowstacktop < num_active_rows); + rowstack[rowstacktop++] = i0; + row_to_component[i0] = k; + assert(rowstacktop < num_active_rows); + rows_by_component[num_rows_done++] = i0; + component_to_num_rows[k]++; + while (rowstacktop || colstacktop) { + while (rowstacktop) { + Py_ssize_t i = rowstack[--rowstacktop]; + TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); + for (Py_ssize_t j = 0; j < N; j++) { + if (col_active[j] && col_to_component[j] < 0 + && !TagInt_is_zero(vec[j])) + { + colstack[colstacktop++] = j; + col_to_component[j] = k; + cols_by_component[num_cols_done++] = j; + component_to_num_cols[k]++; + } + } + } + while (colstacktop) { + Py_ssize_t j = colstack[--colstacktop]; + for (Py_ssize_t i = 0; i < R; i++) { + if (row_active[i] && row_to_component[i] < 0 + && !TagInt_is_zero(Vector_get_vec(PyList_GET_ITEM(arg, i))[j])) + { + rowstack[rowstacktop++] = i; + row_to_component[i] = k; + rows_by_component[num_rows_done++] = i; + component_to_num_rows[k]++; + } + } + } } } - for (Py_ssize_t k = 0; k < num_components; k++) { - assert(components[k] == &cols_by_component[component_to_index[k]]); - } - PyMem_Free(uf_parent); - PyMem_Free(uf_size); - PyMem_Free(root_to_component_index); - PyMem_Free(component_to_index); + assert(num_rows_done == num_active_rows); + assert(num_cols_done == num_active_cols); + PyMem_Free(rowstack); PyMem_Free(colstack); + PyMem_Free(row_to_component); PyMem_Free(col_to_component); *p_num_components = num_components; - *p_component_to_size = component_to_size; + *p_component_to_num_rows = component_to_num_rows; + *p_component_to_num_cols = component_to_num_cols; + *p_rows_by_component = rows_by_component; *p_cols_by_component = cols_by_component; - *p_components = components; return false; -error: - PyMem_Free(uf_parent); - PyMem_Free(uf_size); - PyMem_Free(root_to_component_index); - PyMem_Free(component_to_size); - PyMem_Free(component_to_index); - PyMem_Free(cols_by_component); - PyMem_Free(components); - return true; } static bool @@ -4534,7 +4537,8 @@ static PyObject * decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) { if (!PyList_CheckExact(arg)) { - PyErr_SetString(PyExc_TypeError, "decompose_relations_among(vecs) argument must be a list"); + PyErr_SetString(PyExc_TypeError, + "decompose_relations_among(vecs) argument must be a list"); return NULL; } Py_ssize_t R = PyList_GET_SIZE(arg); @@ -4545,14 +4549,16 @@ decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) for (Py_ssize_t i = 0; i < R; i++) { PyObject *v = PyList_GET_ITEM(arg, i); if (Py_TYPE(v) != &Vector_Type) { - PyErr_SetString(PyExc_TypeError, "decompose_relations_among(vecs) argument must be a list of Vectors"); + PyErr_SetString(PyExc_TypeError, + "decompose_relations_among(vecs) argument must be a list of Vectors"); return NULL; } if (i == 0) { N = Py_SIZE(v); } else { if (Py_SIZE(v) != N) { - PyErr_SetString(PyExc_ValueError, "length mismatch in decompose_relations_among"); + PyErr_SetString(PyExc_ValueError, + "length mismatch in decompose_relations_among"); return NULL; } } @@ -4568,11 +4574,10 @@ decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) char *row_active = NULL; char *col_active = NULL; Py_ssize_t *sorted_rows = NULL; - Py_ssize_t *component_to_size = NULL; + Py_ssize_t *component_to_num_rows = NULL; + Py_ssize_t *component_to_num_cols = NULL; + Py_ssize_t *rows_by_component = NULL; Py_ssize_t *cols_by_component = NULL; - Py_ssize_t **components = NULL; - Py_ssize_t *rows = NULL; - if (!(easy_relations_list = PyList_New(0))) { goto error; } @@ -4582,10 +4587,7 @@ decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) if (!(subproblems_list = PyList_New(0))) { goto error; } - if (!(row_active = PyMem_Malloc(R)) || - !(col_active = PyMem_Malloc(N)) || - !(rows = PyMem_Malloc(R * sizeof(Py_ssize_t))) - ) { + if (!(row_active = PyMem_Malloc(R)) || !(col_active = PyMem_Malloc(N))) { PyErr_NoMemory(); return NULL; } @@ -4617,17 +4619,13 @@ decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) if (!(sorted_rows = mergesort_vectors(arg))) { goto error; } - for (Py_ssize_t ii0 = 0; ii0 < R;) { + for (Py_ssize_t ii0 = 0; ii0 < R; ) { Py_ssize_t ii1 = ii0 + 1; for (; ii1 < R; ii1++) { // extend the streak sorted_rows[ii0:ii1] assert(ii0 < R && ii1 < R); int cmp; - if (Vector_cmp( - Vector_get_vec(PyList_GET_ITEM(arg, sorted_rows[ii0])), - Vector_get_vec(PyList_GET_ITEM(arg, sorted_rows[ii1])), - N, &cmp - )) { + if (Vector_cmp(arg, sorted_rows[ii0], sorted_rows[ii1], N, &cmp)) { goto error; } assert(cmp <= 0); @@ -4661,39 +4659,28 @@ decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) // Step 4: Find which columns have nonzero entries in the same row. Py_ssize_t num_components; if (partition_active_columns( - arg, col_active, row_active, - &num_components, &component_to_size, &cols_by_component, &components + arg, col_active, row_active, &num_components, + &component_to_num_rows, &component_to_num_cols, + &rows_by_component, &cols_by_component )) { goto error; } // Step 5: Make each subproblem of the decomposition + Py_ssize_t num_rows_done = 0; + Py_ssize_t num_cols_done = 0; for (Py_ssize_t k = 0; k < num_components; k++) { - Py_ssize_t num_cols = component_to_size[k]; - Py_ssize_t *cols = components[k]; - Py_ssize_t num_rows = 0; - for (Py_ssize_t i = 0; i < R; i++) { - // find the rows associated to this component. - if (!row_active[i]) { - continue; - } - TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i)); - Py_ssize_t jj = 0; - for (; jj < num_cols; jj++) { - if (!TagInt_is_zero(vec[cols[jj]])) { - break; - } - } - if (jj < num_cols) { - // broke out early: found a nonzero entry - rows[num_rows++] = i; - } - } + Py_ssize_t num_cols = component_to_num_cols[k]; + Py_ssize_t num_rows = component_to_num_rows[k]; + Py_ssize_t *cols = &cols_by_component[num_cols_done]; + Py_ssize_t *rows = &rows_by_component[num_rows_done]; if (append_subproblem( arg, cols, num_cols, rows, num_rows, subproblem_rows_list, subproblems_list )) { goto error; } + num_cols_done += num_cols; + num_rows_done += num_rows; } result = PyTuple_Pack(3, easy_relations_list, subproblem_rows_list, @@ -4705,10 +4692,10 @@ decompose_relations_among(PyObject *Py_UNUSED(mod), PyObject *arg) PyMem_Free(row_active); PyMem_Free(col_active); PyMem_Free(sorted_rows); - PyMem_Free(component_to_size); + PyMem_Free(component_to_num_rows); + PyMem_Free(component_to_num_cols); + PyMem_Free(rows_by_component); PyMem_Free(cols_by_component); - PyMem_Free(components); - PyMem_Free(rows); return result; } @@ -4749,8 +4736,8 @@ static PyMethodDef mutable_lattice_methods[] = { "generalized_row_op(v, w, a, b, c, d) does (v, w) = (a*v+b*w, c*v+d*w) when v and w are Vectors"}, {"xgcd", (PyCFunction)(void(*)(void))xgcd, METH_FASTCALL, "xgcd(a, b) returns a triple (x, y, g) of integers with x*a + y*b == g == gcd(a, b)"}, - {"relations_among", relations_among, METH_O, - "relations_among([v0, ..., vk]) returns the Lattice of coefficient vectors for linear dependencies among the given Vectors"}, + {"relations_among_c", relations_among, METH_O, + "relations_among_c([v0, ..., vk]) returns the Lattice of coefficient vectors for linear dependencies among the given Vectors"}, {"decompose_relations_among", decompose_relations_among, METH_O, "decompose_relations_among([v0, ..., vk]) reduces a relations_among calculation to sub-problems. " "It returns a 3-tuple (easy_relations: list[Vector], subproblem_rows : list[Vector], subproblems: list[list[Vector]])"}, diff --git a/mutable_lattice/__init__.py b/mutable_lattice/__init__.py index d68f10e..899ff4e 100644 --- a/mutable_lattice/__init__.py +++ b/mutable_lattice/__init__.py @@ -4,7 +4,7 @@ generalized_row_op, Lattice, xgcd, - relations_among as relations_among_c, + relations_among_c, decompose_relations_among, transpose, ) diff --git a/mutable_lattice/test.py b/mutable_lattice/test.py index 2aa073a..39480c1 100644 --- a/mutable_lattice/test.py +++ b/mutable_lattice/test.py @@ -369,11 +369,10 @@ def test_generalized_row_op_single(self): def test_generalized_row_op_many(self): VALUES = list(self.VALUES) - random.Random(0).shuffle(VALUES) vx_data = [x for x in VALUES for _ in VALUES] vy_data = [y for _ in VALUES for y in VALUES] - for i in range(len(VALUES) - 4): - a, b, c, d = VALUES[i:i+4] + for _ in range(20): + a, b, c, d = random.choices(VALUES, k=4) vx = Vector(vx_data) vy = Vector(vy_data) new_vx = a*vx + b*vy @@ -1453,8 +1452,8 @@ def test_decompose_relations_among_basic(self): ([Vector([1])], [], [])) self.assertEqual(decompose_relations_among([Vector([1,2,3])]), ([], [], [])) + # find zero rows self.assertEqual( - # find zero rows decompose_relations_among([ Vector([]), Vector([]), @@ -1598,6 +1597,47 @@ def test_decompose_relations_among_basic(self): ([Vector([1, -1, 0, 0]), Vector([0, 0, 1, -1])], [], []) ) + def test_decompose_relations_among_random(self): + def random_cases(): + for R in range(5): + for N in range(5): + for num_nonzero in range(N*R): + for _ in range(3): + xs = random.choices((-1,0,1,2), k=num_nonzero) + places = random.sample(range(N*R), k=num_nonzero) + data = [Vector.zero(N) for _ in range(R)] + for place, x in zip(places, xs): + i, j = divmod(place, N) + data[i][j] = x + yield R, N, data + for R, N, data in random_cases(): + easy_relations, subproblem_rows, subproblems = decompose_relations_among(data) + for rel in easy_relations: + nz = list(filter(rel.__getitem__, range(R))) + if len(nz) == 1: + [j] = nz + self.assertEqual(rel[j], 1) + self.assertEqual(data[j], Vector.zero(N)) + else: + [j1, j2] = nz + self.assertEqual([rel[j1], rel[j2]], [1, -1]) + self.assertEqual(data[j1], data[j2]) + total_rank = Lattice(N, data).rank + total_nullity = len(easy_relations) + for rows, subproblem in zip(subproblem_rows, subproblems, strict=True): + for i, v in zip(rows, subproblem): + self.assertEqual(sorted(filter(None, data[i])), + sorted(filter(None, v))) + self.assertTrue(any(v)) + for col in zip(*subproblem): + self.assertGreaterEqual(len(list(filter(None, col))), 2) + for k in relations_among_c(subproblem).get_basis(): + v = Vector.zero(N) + for i, x in zip(rows, k, strict=True): + v += x * data[i] + self.assertEqual(v, Vector.zero(N)) + total_nullity += 1 + self.assertEqual(total_rank + total_nullity, R) class TestTranspose(unittest.TestCase): def setUp(self):