Skip to content

Move row and column sorting to decompose_relations_among #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sweeneyde merged 9 commits into from
Mar 16, 2026
Merged

Move row and column sorting to decompose_relations_among #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
a107ad7
decompose_relations_among: get the easy relations out of the way and …
sweeneyde Mar 10, 2026
0b80512
update pyi
sweeneyde Mar 10, 2026
802ed89
more tests
sweeneyde Mar 11, 2026
6022e9f
Refactor. Do a BFS instead of using a union-find.
sweeneyde Mar 12, 2026
69fffda
move sorting to decompose_relations_among
sweeneyde Mar 15, 2026
5a8a563
merge with main
sweeneyde Mar 15, 2026
2bb06a2
tweak assertions
sweeneyde Mar 15, 2026
e6e5bc6
tweak error handling
sweeneyde Mar 16, 2026
c1bd5ec
stacks-->queues
sweeneyde Mar 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
196 changes: 99 additions & 97 deletions _mutable_lattice.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4003,11 +4003,11 @@ static PyTypeObject Lattice_Type = {
};

/*********************************************************************/
/* relations_among, decompose_relations_among, transpose */
/* relations_among_c, decompose_relations_among, transpose */
/*********************************************************************/

static PyObject *
relations_among(PyObject *Py_UNUSED(mod), PyObject *arg)
relations_among_c(PyObject *Py_UNUSED(mod), PyObject *arg)
{
if (!PyList_CheckExact(arg)) {
PyErr_SetString(PyExc_TypeError, "relations_among(vecs) argument must be a list");
Expand Down Expand Up @@ -4038,91 +4038,46 @@ relations_among(PyObject *Py_UNUSED(mod), PyObject *arg)
// zero rows are at the bottom. To make these zero rows, we had to take
// Z-linear combination of the input rows. The right portion keeps
// track of the coefficients needed to do so.

Py_ssize_t *nonzero_in_column = PyMem_Calloc(N, sizeof(Py_ssize_t));
Py_ssize_t *column_order = NULL;
Py_ssize_t *index_of_first_nonzero = PyMem_Malloc(R * sizeof(Py_ssize_t));
Py_ssize_t *rows_by_index_of_first_nonzero = NULL;
TagInt *scratch = PyMem_Malloc((N + R) * sizeof(TagInt *));
Lattice *L = NULL;
Lattice *relations = NULL;
PyObject *result = NULL;
if (nonzero_in_column == NULL || index_of_first_nonzero == NULL || scratch == NULL) {
Lattice *relations = NULL;
if (scratch == NULL) {
PyErr_NoMemory();
goto error;
}

// A first tweak:
// Earlier columns influence what row operations happen on what that follows them,
// so we reshuffle the columns so that less populated columns come first.
for (Py_ssize_t i = 0; i < R; i++) {
TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i));
for (Py_ssize_t j = 0; j < N; j++) {
if (!TagInt_is_zero(vec[j])) {
nonzero_in_column[j]++;
}
}
}
column_order = counting_sort_order(nonzero_in_column, N);
if (column_order == NULL) {
goto error;
}
Py_ssize_t k0 = 0; // the first nonzero column
while (k0 < N && nonzero_in_column[column_order[k0]] == 0) {
k0++;
}

// A second tweak:
// Sort rows by decreasing index of first nonzero entry. This allows
// the bottom right corner to stay HNF while new rows are added above it.
for (Py_ssize_t i = 0; i < R; i++) {
TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i));
Py_ssize_t k = k0;
while (k < N && TagInt_is_zero(vec[column_order[k]])) {
k++;
}
index_of_first_nonzero[i] = k;
}
rows_by_index_of_first_nonzero = counting_sort_order(index_of_first_nonzero, R);
if (rows_by_index_of_first_nonzero == NULL) {
goto error;
return NULL;
}

// Now add everything in in an appropriate order
L = (Lattice *)Lattice_new_impl(&Lattice_Type, N - k0 + R, 1, R);
Lattice *L = (Lattice *)Lattice_new_impl(&Lattice_Type, N + R, 1, R);
if (L == NULL) {
goto error;
PyMem_Free(scratch);
return NULL;
}
for (Py_ssize_t ii = R-1; ii >= 0; ii--) {
Py_ssize_t i = rows_by_index_of_first_nonzero[ii];
for (Py_ssize_t i = R-1; i >= 0; i--) {
TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, i));
for (Py_ssize_t k = k0; k < N; k++) {
scratch[k - k0] = vec[column_order[k]];
}
memset(scratch + N - k0, 0, R*sizeof(TagInt *));
scratch[N - k0 + i] = TagInt_ONE;
memcpy(scratch, vec, N*sizeof(TagInt));
memset(scratch + N, 0, R*sizeof(TagInt *));
scratch[N + i] = TagInt_ONE;
if (Lattice_add_vector_impl(L, scratch)) {
goto error;
}
if (R > 1000 && L->rank % 16 == 0 && PyErr_CheckSignals() < 0) {
goto error;
}
if (R > 1000 && L->rank % 100 == 0) {
printf("relations_among_c: %ld/%ld vectors added\n", (long)L->rank, (long)R);
}
}
relations = (Lattice *)Lattice_new_impl(&Lattice_Type, R, 1, R);
if (relations == NULL) {
goto error;
}
assert(L->rank == R);
for (Py_ssize_t i = R - 1; i >= 0; i--) {
if (L->row_to_pivot[i] >= N - k0) {
if (Lattice_add_vector_impl(relations, L->basis[i] + N - k0)) {
goto error;
}
for (Py_ssize_t i = R - 1; i >= 0 && L->row_to_pivot[i] >= N; i--) {
if (Lattice_add_vector_impl(relations, L->basis[i] + N)) {
goto error;
}
}
result = Py_NewRef((PyObject *)relations);
error:
PyMem_Free(nonzero_in_column);
PyMem_Free(column_order);
PyMem_Free(index_of_first_nonzero);
PyMem_Free(rows_by_index_of_first_nonzero);
PyMem_Free(scratch);
Py_XDECREF(L);
Py_XDECREF(relations);
Expand Down Expand Up @@ -4400,24 +4355,23 @@ partition_active_columns(
}
}
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_queue = PyMem_Malloc(num_active_rows * sizeof(Py_ssize_t));
Py_ssize_t *col_queue = 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 &&
if (!(row_queue && col_queue &&
row_to_component && col_to_component &&
component_to_num_rows && component_to_num_cols &&
rows_by_component && cols_by_component
)) {
PyErr_NoMemory();
PyMem_Free(rowstack); PyMem_Free(colstack);
PyMem_Free(row_queue); PyMem_Free(col_queue);
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);
Expand All @@ -4429,52 +4383,56 @@ partition_active_columns(
for (Py_ssize_t j = 0; j < N; j++) {
col_to_component[j] = -1;
}
// Explore the bipartite graph of rows and columns
// BFS 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 row_queue_length = 1, col_queue_length = 0;
Py_ssize_t k = num_components++;
assert(rowstacktop < num_active_rows);
rowstack[rowstacktop++] = i0;
row_queue[0] = 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];
do {
assert(row_queue_length <= num_active_rows);
for (Py_ssize_t ii = 0; ii < row_queue_length; ii++) {
Py_ssize_t i = row_queue[ii];
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_queue[col_queue_length++] = 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];
row_queue_length = 0;
assert(col_queue_length <= num_active_cols);
for (Py_ssize_t jj = 0; jj < col_queue_length; jj++) {
Py_ssize_t j = col_queue[jj];
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_queue[row_queue_length++] = i;
row_to_component[i] = k;
rows_by_component[num_rows_done++] = i;
component_to_num_rows[k]++;
}
}
}
}
col_queue_length = 0;
} while (row_queue_length);
}
assert(num_rows_done == num_active_rows);
assert(num_cols_done == num_active_cols);
PyMem_Free(rowstack); PyMem_Free(colstack);
assert(num_components <= max_num_components);
PyMem_Free(row_queue); PyMem_Free(col_queue);
PyMem_Free(row_to_component); PyMem_Free(col_to_component);
*p_num_components = num_components;
*p_component_to_num_rows = component_to_num_rows;
Expand All @@ -4493,44 +4451,88 @@ append_subproblem(
PyObject *subproblems_list
)
{
// shuffle the columns to put the most-used columns first
Py_ssize_t *col_to_count = PyMem_Calloc(num_cols, sizeof(Py_ssize_t));
Py_ssize_t *row_to_index = PyMem_Malloc(num_rows * sizeof(Py_ssize_t));
if (col_to_count == NULL || row_to_index == NULL) {
PyErr_NoMemory();
return true;
}
for (Py_ssize_t ii = 0; ii < num_rows; ii++) {
TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, rows[ii]));
for (Py_ssize_t jj = 0; jj < num_cols; jj++) {
if (!TagInt_is_zero(vec[cols[jj]])) {
col_to_count[jj]++;
}
}
}
Py_ssize_t *col_order = counting_sort_order(col_to_count, num_cols);
PyMem_Free(col_to_count);
if (col_order == NULL) {
PyMem_Free(row_to_index);
return true;
}
// sort the rows by increasing index of first nonzero entry
for (Py_ssize_t ii = 0; ii < num_rows; ii++) {
TagInt *vec = Vector_get_vec(PyList_GET_ITEM(arg, rows[ii]));
Py_ssize_t jjj = 0;
while (TagInt_is_zero(vec[cols[col_order[jjj]]])) {
jjj++;
}
assert(jjj < num_cols);
row_to_index[ii] = jjj;
}
Py_ssize_t *row_order = counting_sort_order(row_to_index, num_rows);
PyMem_Free(row_to_index);
if (row_order == NULL) {
PyMem_Free(col_order);
return true;
}
// Here we begin using "goto error"
PyObject *subproblem_rows = Vector_zero_impl(num_rows);
if (subproblem_rows == NULL) {
return true;
goto error;
}
TagInt *rows_vec = Vector_get_vec(subproblem_rows);
for (Py_ssize_t ii = 0; ii < num_rows; ii++) {
for (Py_ssize_t iii = 0; iii < num_rows; iii++) {
static_assert(sizeof(intptr_t) == sizeof(Py_ssize_t));
rows_vec[ii] = pack_integer((intptr_t)rows[ii]);
rows_vec[iii] = pack_integer((intptr_t)rows[row_order[iii]]);
}
if (PyList_Append(subproblem_rows_list, subproblem_rows) < 0) {
Py_DECREF(subproblem_rows);
return true;
goto error;
}
Py_DECREF(subproblem_rows);

PyObject *subproblem = PyList_New(num_rows);
if (subproblem == NULL) {
return true;
goto error;
}
for (Py_ssize_t ii = 0; ii < num_rows; ii++) {
for (Py_ssize_t iii = 0; iii < num_rows; iii++) {
PyObject *v = Vector_zero_impl(num_cols);
if (v == NULL) {
Py_DECREF(subproblem);
return true;
goto error;
}
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]]);
TagInt *source = Vector_get_vec(PyList_GET_ITEM(arg, rows[row_order[iii]]));
for (Py_ssize_t jjj = 0; jjj < num_cols; jjj++) {
dest[jjj] = TagInt_copy(source[cols[col_order[jjj]]]);
}
PyList_SET_ITEM(subproblem, ii, v);
PyList_SET_ITEM(subproblem, iii, v);
}
if (PyList_Append(subproblems_list, subproblem) < 0) {
Py_DECREF(subproblem);
return true;
goto error;
}
Py_DECREF(subproblem);
PyMem_Free(row_order);
PyMem_Free(col_order);
return false;
error:
PyMem_Free(row_order);
PyMem_Free(col_order);
return true;
}

static PyObject *
Expand Down Expand Up @@ -4736,7 +4738,7 @@ 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_c", relations_among, METH_O,
{"relations_among_c", relations_among_c, 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. "
Expand Down
9 changes: 7 additions & 2 deletions mutable_lattice/test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1629,8 +1629,13 @@ def random_cases():
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)
rn = range(len(subproblem[0]))
indexes = [next(filter(v.__getitem__, rn)) for v in subproblem]
self.assertEqual(indexes, sorted(indexes))
cols = transpose(len(subproblem[0]), subproblem)
nonzeros_in_cols = [len(list(filter(None, col))) for col in cols]
self.assertGreaterEqual(min(nonzeros_in_cols), 2)
self.assertEqual(nonzeros_in_cols, sorted(nonzeros_in_cols))
for k in relations_among_c(subproblem).get_basis():
v = Vector.zero(N)
for i, x in zip(rows, k, strict=True):
Expand Down
Loading