Commit 75b300f3 authored by Charles Bouillaguet's avatar Charles Bouillaguet
Browse files

don't precompute multiplication tables anymore

parent 28cc660b
......@@ -196,7 +196,6 @@ void parser_store_polynomial(void *opaque, int line)
typedef u128 monomial_t;
monomial_t *monomials; /* array to rank/unrank monomials */
int *mul_table; /* multiplication table */
monomial_t vmask; /* monomial that contains all special variables */
int nbad; /* total number of bad monomials */
int *bad_renumbering; /* renumber the bad monomials consecutively */
......@@ -312,36 +311,6 @@ void test_grlex_order()
}
}
/* multiplies a degree-(<= D-2) monomial by a degre-(<= 2) monomial */
static inline int mmul(int a, int b)
{
assert(b <= Nd[3]);
assert(a <= Nd[D - 1]);
return mul_table[a * Nd[3] + b];
}
void test_mmul()
{
for (int i = 0; i < n; i++) {
assert(mmul(0, i + 1) == i + 1); /* multiplication by 1 */
if (D <= 2)
continue;
assert(mmul(i + 1, 0) == i + 1); /* multiplication by 1 */
assert(mmul(i + 1, i + 1) == i + 1); /* multiplication by x_i */
}
for (int i = 0; i < Nd[3]; i++)
assert(mmul(0, i) == i);
if (D <= 2)
return;
for (int i = 0; i < n; i++)
for (int j = 0; j < i; j++) {
int u = mmul(i + 1, j + 1);
int v = mmul(j + 1, i + 1);
assert(u == v);
assert(hamming_weight(monomials[u]) == 2);
}
}
/* returns true when m contains more than one special variable */
static inline bool monomial_is_bad(monomial_t m)
{
......@@ -403,24 +372,6 @@ void monomial_setup_list()
log(1, " - Monomials setup: %.1fs\n", wtime() - start);
}
void monomial_setup_multiplication()
{
log(0, "Building monomial multiplication table\n");
u64 size = ((u64) Nd[3]) * ((u64) Nd[D - 1]) * sizeof(*mul_table);
mul_table = malloc(size);
if (mul_table == NULL)
err(1, "cannot allocate monomial multiplication table");
for (int i = 0; i < Nd[3]; i++) {
#pragma omp parallel for
for (int j = 0; j < Nd[D - 1]; j++) {
monomial_t product = monomials[i] | monomials[j];
mul_table[j * Nd[3] + i] = monomial_rank(product);
}
}
test_mmul();
}
/******************************* input system in nicer form *****************************/
struct input_poly_t {
......@@ -602,9 +553,15 @@ void macaulay_dense(int d)
int row = 0;
/* compute m_j * f_i, for all monomial m_j of degree <= d */
for (int j = 0; j < Nd[d - 1]; j++) {
/* generate multiplication table for m_j */
int mmul[9000];
for (int k = 0; k < Nd[3]; k++) {
monomial_t product = monomials[j] | monomials[k];
mmul[k] = monomial_rank(product);
}
for (int i = 0; i < m; i++) {
for (int k = 0; k < input_system[i].nterms; k++) {
int l = mmul(j, input_system[i].terms[k]);
int l = mmul[input_system[i].terms[k]];
mzd_flip_bit(M, row, M->ncols - 1 - l);
}
row += 1;
......@@ -644,11 +601,19 @@ void macaulay_sparse(int d)
skipped += 1;
continue;
}
/* generate multiplication table for m_j */
int mmul[9000];
for (int k = 0; k < Nd[3]; k++) {
monomial_t product = monomials[j] | monomials[k];
mmul[k] = monomial_rank(product);
}
for (int i = 0; i < m; i++) {
/* compute and store m_j * f_i */
int row[9000];
int size = input_system[i].nterms;
for (int k = 0; k < size; k++)
row[k] = mmul(j, input_system[i].terms[k]);
row[k] = mmul[input_system[i].terms[k]];
matrix_store_row(i, j, row, size, scratch);
}
}
......@@ -723,7 +688,6 @@ int main(int argc, char **argv)
log(1, " - %d total monomials of degree <= %d\n", Nd[D+1], D);
/* prepare monomial multiples and setup sparse matrix */
monomial_setup_multiplication();
macaulay_dense(D - 2);
log(1, " - %d monomial multipliers\n", valid_multipliers);
nrows = m * valid_multipliers;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment