+ isl_mat_free(M);
+ isl_mat_free(C);
+ isl_mat_free(U);
+ return NULL;
+}
+
+/* Compute and return the matrix
+ *
+ * U_1^{-1} diag(d_1, 1, ..., 1)
+ *
+ * with U_1 the unimodular completion of the first (and only) row of B.
+ * The columns of this matrix generate the lattice that satisfies
+ * the single (linear) modulo constraint.
+ */
+static struct isl_mat *parameter_compression_1(
+ struct isl_mat *B, struct isl_vec *d)
+{
+ struct isl_mat *U;
+
+ U = isl_mat_alloc(B->ctx, B->n_col - 1, B->n_col - 1);
+ if (!U)
+ return NULL;
+ isl_seq_cpy(U->row[0], B->row[0] + 1, B->n_col - 1);
+ U = isl_mat_unimodular_complete(U, 1);
+ U = isl_mat_right_inverse(U);
+ if (!U)
+ return NULL;
+ isl_mat_col_mul(U, 0, d->block.data[0], 0);
+ U = isl_mat_lin_to_aff(U);
+ return U;
+}
+
+/* Compute a common lattice of solutions to the linear modulo
+ * constraints specified by B and d.
+ * See also the documentation of isl_mat_parameter_compression.
+ * We put the matrix
+ *
+ * A = [ L_1^{-T} L_2^{-T} ... L_k^{-T} ]
+ *
+ * on a common denominator. This denominator D is the lcm of modulos d.
+ * Since L_i = U_i^{-1} diag(d_i, 1, ... 1), we have
+ * L_i^{-T} = U_i^T diag(d_i, 1, ... 1)^{-T} = U_i^T diag(1/d_i, 1, ..., 1).
+ * Putting this on the common denominator, we have
+ * D * L_i^{-T} = U_i^T diag(D/d_i, D, ..., D).
+ */
+static struct isl_mat *parameter_compression_multi(
+ struct isl_mat *B, struct isl_vec *d)
+{
+ int i, j, k;
+ isl_int D;
+ struct isl_mat *A = NULL, *U = NULL;
+ struct isl_mat *T;
+ unsigned size;
+
+ isl_int_init(D);
+
+ isl_vec_lcm(d, &D);
+
+ size = B->n_col - 1;
+ A = isl_mat_alloc(B->ctx, size, B->n_row * size);
+ U = isl_mat_alloc(B->ctx, size, size);
+ if (!U || !A)
+ goto error;
+ for (i = 0; i < B->n_row; ++i) {
+ isl_seq_cpy(U->row[0], B->row[i] + 1, size);
+ U = isl_mat_unimodular_complete(U, 1);
+ if (!U)
+ goto error;
+ isl_int_divexact(D, D, d->block.data[i]);
+ for (k = 0; k < U->n_col; ++k)
+ isl_int_mul(A->row[k][i*size+0], D, U->row[0][k]);
+ isl_int_mul(D, D, d->block.data[i]);
+ for (j = 1; j < U->n_row; ++j)
+ for (k = 0; k < U->n_col; ++k)
+ isl_int_mul(A->row[k][i*size+j],
+ D, U->row[j][k]);
+ }
+ A = isl_mat_left_hermite(A, 0, NULL, NULL);
+ T = isl_mat_sub_alloc(A, 0, A->n_row, 0, A->n_row);
+ T = isl_mat_lin_to_aff(T);
+ if (!T)
+ goto error;
+ isl_int_set(T->row[0][0], D);
+ T = isl_mat_right_inverse(T);
+ if (!T)
+ goto error;
+ isl_assert(T->ctx, isl_int_is_one(T->row[0][0]), goto error);
+ T = isl_mat_transpose(T);
+ isl_mat_free(A);
+ isl_mat_free(U);
+
+ isl_int_clear(D);
+ return T;
+error:
+ isl_mat_free(A);
+ isl_mat_free(U);
+ isl_int_clear(D);