Commit 69dd1117 authored by Quentin HAMMERER's avatar Quentin HAMMERER
Browse files

abc

parent bc11158d
CC = gcc
CFLAGS = -g -Wall -Wextra -O3 -Wno-unused-parameter -Wno-maybe-uninitialized -Werror -march=native
CFLAGS = -g -Wall -Wextra -O3 -Wno-unused-parameter -Wno-maybe-uninitialized -Werror -march=native
# OpenMP
CFLAGS += -fopenmp
LDFLAGS += -fopenmp
LDLIBS += -lm
all: parser_demo parser_demo_alt moebius_demo monica monica_vector macaulay_gen lanczos reconstruct_monomials
all: macaulay_gen reconstruct_monomials
monica: parser.o monica.o
monica_vector: parser.o monica_vector.o
parser_demo: parser.o parser_demo.o
parser_demo_alt: parser.o parser_demo_alt.o
moebius_demo: parser.o moebius.o moebius_demo.o
macaulay_gen: parser.o tools.o monomials.o macaulay_gen.o
macaulay_gen: parser.o tools.o tools_parser.o monomials.o macaulay_gen.o
macaulay_gen: LDLIBS += `pkg-config --libs m4ri`
reconstruct_monomials: parser.o tools.o monomials.o reconstruct_monomials.o
reconstruct_monomials: parser.o tools.o tools_parser.o monomials.o reconstruct_monomials.o
reconstruct_monomials: LDLIBS += `pkg-config --libs m4ri`
#test: parser.o tools.o tools_parser.o monomials.o test.o
#test: LDLIBS += `pkg-config --libs m4ri`
lanczos.o: CFLAGS += -flax-vector-conversions
moebius_demo.o: parser.h moebius.h
parser_demo.o: parser.h
parser.o: parser.h
tools.o: tools.h
monomials.o: monomials.h
moebius.o: moebius.h
monica.o: parser.h ffs.h
monica_vector.o: parser.h ffs.h
macaulay_gen.o: parser.h tools.h monomials.h
monomials.o: monomials.h tools.h
tools_parser.o: tools_parser.h monomials.h tools.h
macaulay_gen.o: parser.h tools_parser.h tools.h monomials.h
macaulay_gen.o: CFLAGS += `pkg-config --cflags m4ri`
reconstruct_monomials.o: parser.h tools.h monomials.h
reconstruct_monomials.o: CFLAGS += `pkg-config --cflags m4ri`
reconstruct_monomials.o: parser.h tools_parser.h tools.h monomials.h
reconstruct_monomials.o: CFLAGS += `pkg-config --cflags m4ri`
#test.o: parser.h tools_parser.h tools.h monomials.h
#test.o: CFLAGS += `pkg-config --cflags m4ri`
clean:
rm -rf *.o my_mat/ *log.bin
rm -f parser_demo moebius_demo parser_demo_alt monica monica_vector macaulay_gen lanczos reconstruct_monomials
rm -f parser_demo moebius_demo parser_demo_alt monica monica_vector macaulay_gen lanczos reconstruct_monomials test
......@@ -14,8 +14,9 @@
#include <m4ri/m4ri.h>
#include "parser.h"
#include "monomials.h"
#include "tools.h"
// #include "monomials.h"
// #include "tools.h"
#include "tools_parser.h"
/*
* This program generates a degree-D macaulay matrix from a quadratic system,
......@@ -35,54 +36,33 @@ int Nd_l; /* Nd[D-1] */
// int Nd[MAXDEG]; /* number of degree-d monomials */
// monomial_t *monomials; /* array to rank/unrank monomials */
monomial_t vmask; /* monomial that contains all special variables */
int nbad; /* total number of bad monomials */
int *bad_renumbering; /* renumber the bad monomials consecutively */
int *inv_bad_renumbering; /* renumber the bad monomials consecutively */
void test_bad_monomials()
void test_bad_monomials(struct monomial_ctx_t *mono)
{
assert(bad_renumbering[0] < 0);
for (int i = 0; i < n; i++) {
assert(bad_renumbering[i + 1] < 0);
}
for (int i = 0; i < v; i++)
{
for (int j = 0; j < i; j++)
{
monomial_t inter = (mkvar(i) | mkvar(j)) & vmask;
assert((inter & (inter - 1)) != 0);
}
}
for (int i = v; i < n; i++)
{
for (int j = 0; j < i; j++)
{
monomial_t inter = (mkvar(i) | mkvar(j)) & vmask;
assert((inter & (inter - 1)) == 0);
}
}
test_bad_monomials_2(mono);
}
void renumbering_setup()
void renumbering_setup(struct monomial_ctx_t *mono)
{
vmask = (v == 0) ? 0 : mkvar(v) - 1;
log(1, " - Singling out bad monomials\n");
bad_renumbering = alloc_array(Nd_p, sizeof(*bad_renumbering), "bad monomial renumbering");
inv_bad_renumbering = alloc_array(Nd_p, sizeof(*bad_renumbering), "bad monomial renumbering");
for (int j = 0; j < Nd_p; j++) {
bool bad = monomial_is_bad(j, vmask);
bool bad = monomial_is_bad(j, mono);
bad_renumbering[j] = bad ? nbad : -1;
inv_bad_renumbering[nbad] = j;
nbad += bad;
}
log(2, " - %d bad monomials\n", nbad);
test_bad_monomials();
test_bad_monomials(mono);
}
......@@ -237,14 +217,14 @@ static inline int mzd_ffs(mzd_t *M, rci_t row)
}
void macaulay_dense(int d)
void macaulay_dense(int d, struct monomial_ctx_t *mono)
{
int Nd_p = get_Nd(d+1);
int Nd_l = get_Nd(d-1);
int Nd_p = get_Nd(d+1, mono);
int Nd_l = get_Nd(d-1, mono);
valid_multipliers = Nd_p;
if (d < 2)
return;
int Nd_3 = get_Nd(3);
int Nd_3 = get_Nd(3, mono);
log(0, "Generating degree-%d (dense) Macaulay matrix \n", d);
log(1, " - matrix size: %d x %d\n", m * Nd_l, Nd_p);
M = mzd_init(m * Nd_l, Nd_p);
......@@ -254,10 +234,8 @@ void macaulay_dense(int d)
for (int j = 0; j < Nd_l; j++) {
/* generate multiplication table for m_j */
int mmul[9000];
for (int k = 0; k < Nd_3; k++) {
monomial_t product = monomials_product(j, k);
mmul[k] = monomial_rank(product);
}
for (int k = 0; k < Nd_3; k++)
mmul[k] = monomials_product(j, k, mono);
for (int i = 0; i < m; i++) {
for (int k = 0; k < input_system[i].nterms; k++) {
int l = mmul[input_system[i].terms[k]];
......@@ -287,15 +265,15 @@ void macaulay_dense(int d)
}
/* multiply by degree-d monomials. May run inside a parallel region */
void macaulay_sparse(int d)
void macaulay_sparse(int d, struct monomial_ctx_t *mono)
{
bool * scratch = alloc_array(ncols, sizeof(*scratch), "scratch (per thread)");
for (int j = 0; j < ncols; j++)
scratch[j] = 0;
int Nd = get_Nd(d);
int Nd_p = get_Nd(d+1);
int Nd_3 = get_Nd(3);
int Nd = get_Nd(d, mono);
int Nd_p = get_Nd(d+1, mono);
int Nd_3 = get_Nd(3, mono);
#pragma omp for
for (int j = Nd; j < Nd_p; j++) {
if (LM && LM[j]) {
......@@ -306,8 +284,7 @@ void macaulay_sparse(int d)
/* generate multiplication table for m_j */
int mmul[9000];
for (int k = 0; k < Nd_3; k++) {
monomial_t product = monomials_product(j, k);
mmul[k] = monomial_rank(product);
mmul[k] = monomials_product(j, k, mono);
}
for (int i = 0; i < m; i++) {
......@@ -387,16 +364,16 @@ int main(int argc, char **argv)
log(1, " - Truncating columns linear in the first %d variables\n", v);
monomial_setup(v, D);
Nd_p = get_Nd(D+1);
Nd_l = get_Nd(D-1);
renumbering_setup();
struct monomial_ctx_t *mono = monomial_setup(n, v, D);
Nd_p = get_Nd(D+1, mono);
Nd_l = get_Nd(D-1, mono);
renumbering_setup(mono);
input_system = convert_input_system(Nd_p, &m);
input_system = convert_input_system(Nd_p, &m, mono);
log(1, " - %d total monomials of degree <= %d\n", Nd_p, D);
/* prepare monomial multiples and setup sparse matrix */
macaulay_dense(D - 2);
macaulay_dense(D - 2, mono);
log(1, " - %d monomial multipliers\n", valid_multipliers);
nrows = m * valid_multipliers;
ncols = nbad;
......@@ -420,7 +397,7 @@ int main(int argc, char **argv)
log(1, " - Expanding to degree %d\n", d);
#pragma omp parallel
macaulay_sparse(d - 2);
macaulay_sparse(d - 2, mono);
char hnnz[10];
human_format(rows_written, hnrows);
......@@ -498,7 +475,7 @@ int main(int argc, char **argv)
log(1, " - Expanding to degree %d\n", d);
#pragma omp parallel
macaulay_sparse(d - 2);
macaulay_sparse(d - 2, mono);
char hnnz[10];
human_format(rows_written, hnrows);
......@@ -512,7 +489,7 @@ int main(int argc, char **argv)
fclose(buffer_stream);
free(inv_bad_renumbering);
free(bad_renumbering);
free_monomial_ctx();
free_monomial_ctx(mono);
write_array(rw_filename, rw, nrows);
write_array(cw_filename, cw, ncols);
write_array(ilog_filename, ilog, nrows);
......
#include "monomials.h"
#include "tools.h"
typedef __int128 u128;
typedef u128 monomial_t;
struct monomial_ctx_t{
int v;
int D;
int Nd[MAXDEG];
struct monomial_ctx_t
{
monomial_t *monomials;
int Nd[MAXDEG];
int D;
int n_var;
int v;
};
struct monomial_ctx_t *mono;
static inline monomial_t mkvar(int i)
{
return ((monomial_t) 1) << i;
}
bool monomial_is_bad(int index, monomial_t vmask)
int hamming_weight(monomial_t m)
{
monomial_t m = mono->monomials[index];
monomial_t inter = m & vmask;
return (inter & (inter - 1)) != 0;
return __builtin_popcountl(m) + __builtin_popcountl(m >> 64);
}
void monomial_setup_list()
/* returns the position of target under the assumption that it is in monomials[lo:hi] */
int binary_search(monomial_t target, int lo, int hi, struct monomial_ctx_t *mono)
{
hi -= 1;
while (lo <= hi)
{
int mid = (lo + hi) / 2;
if (target > mono->monomials[mid])
lo = mid + 1;
else if (target < mono->monomials[mid])
hi = mid - 1;
else
return mid;
}
assert(0);
}
int monomial_rank(monomial_t m, struct monomial_ctx_t *mono)
{
int d = hamming_weight(m);
return binary_search(m, mono->Nd[d], mono->Nd[d + 1], mono);
}
void test_grlex_order(struct monomial_ctx_t *mono)
{
assert(hamming_weight(mono->monomials[0]) == 0);
for (int i = 0; i < mono->n_var; i++)
assert(hamming_weight(mono->monomials[1 + i]) == 1);
for (int d = 0; d <= mono->D; d++)
for (int k = mono->Nd[d]; k < mono->Nd[d + 1] - 1; k++) {
assert(hamming_weight(mono->monomials[k]) == d);
assert(mono->monomials[k] < mono->monomials[k + 1]);
}
}
void test_ranking(struct monomial_ctx_t *mono)
{
for (int i = 0; i < mono->Nd[mono->D + 1]; i++)
assert(monomial_rank(mono->monomials[i], mono) == i);
}
int monomial_enumerate(int degree, monomial_t *out, int n_var)
{
if (degree == 0) {
if (out != NULL)
out[0] = 0;
return 1;
}
int count = 0;
monomial_t top = mkvar(n_var);
monomial_t i = mkvar(degree) - 1;
while (i < top) {
if (out != NULL)
out[count] = i;
count++;
monomial_t minbit = i & -i; /* rightmost 1 bit */
monomial_t fillbit = (i + minbit) & ~i; /* rightmost 0 to the left of that bit */
i = (i + minbit) | ((fillbit / (minbit << 1)) - 1);
}
return count;
}
void monomial_setup_list(struct monomial_ctx_t *mono)
{
int D = mono->D;
......@@ -30,7 +97,7 @@ void monomial_setup_list()
log(1, " - Counting monomials\n");
mono->Nd[0] = 0;
for (int d = 0; d <= D; d++) {
int binom_n_d = monomial_enumerate(d, NULL); /* not the smartest way to do the job... */
int binom_n_d = monomial_enumerate(d, NULL, mono->n_var); /* not the smartest way to do the job... */
mono->Nd[d + 1] = mono->Nd[d] + binom_n_d;
log(2, " - %d monomials of degree %d\n", binom_n_d, d);
}
......@@ -38,76 +105,251 @@ void monomial_setup_list()
log(1, " - Building monomial list\n");
mono->monomials = alloc_array(mono->Nd[D + 1], sizeof(*(mono->monomials)), "monomial list");
for (int d = 0; d <= D; d++)
monomial_enumerate(d, &(mono->monomials[mono->Nd[d]]));
// test_grlex_order(32);
test_ranking();
monomial_enumerate(d, &(mono->monomials[mono->Nd[d]]), mono->n_var);
test_grlex_order(mono);
test_ranking(mono);
log(1, " - Monomials setup: %.1fs\n", wtime() - start);
}
void monomial_setup(int v, int D)
struct monomial_ctx_t *monomial_setup(int n_var, int v, int D)
{
mono = (struct monomial_ctx_t *)malloc(sizeof(struct monomial_ctx_t));
struct monomial_ctx_t *mono = (struct monomial_ctx_t *)malloc(sizeof(struct monomial_ctx_t));
mono->v = v;
mono->D = D;
mono->n_var = n_var;
monomial_setup_list(mono);
return mono;
}
void free_monomial_ctx(struct monomial_ctx_t *mono)
{
free(mono->monomials);
free(mono);
}
bool monomial_is_bad(int index, struct monomial_ctx_t *mono)
{
monomial_t vmask = (mono->v == 0) ? 0 : mkvar(mono->v) - 1;
monomial_t m = mono->monomials[index];
monomial_t inter = m & vmask;
return (inter & (inter - 1)) != 0;
}
void test_bad_monomials_2(struct monomial_ctx_t *mono)
{
for (int i = 0; i < mono->v; i++)
for (int j = 0; j < i; j++)
assert(monomial_is_bad(monomial_rank(mkvar(i) | mkvar(j), mono), mono));
for (int i = mono->v; i < mono->n_var; i++)
for (int j = 0; j < i; j++)
assert(!monomial_is_bad(monomial_rank(mkvar(i) | mkvar(j), mono), mono));
}
void test_shift_128()
{
for (int i = 0; i < 128; i++) {
monomial_t x = 1;
x <<= i;
monomial_t y = mkvar(i);
assert(x == y);
}
}
void test_HW()
{
for (int i = 0; i < 128; i++)
assert(hamming_weight(mkvar(i)) == 1);
for (int i = 0; i < 128; i++)
for (int j = 0; j < i; j++)
assert(hamming_weight(mkvar(i) + mkvar(j)) == 2);
}
int get_Nd(int d)
int get_Nd(int d, struct monomial_ctx_t *mono)
{
if (d <= mono->D + 1)
if (d <= mono->D+1)
return mono->Nd[d];
errx(1, "deg in get_Nd too high");
}
monomial_t monomials_product(int j, int k)
void print_monomial(FILE *stream, int ind, char **variable_name, struct monomial_ctx_t *mono)
{
monomial_t m = mono->monomials[ind];
int v[mono->n_var];
int d = 0;
for (int i = 0; i < mono->n_var; i++) {
if (m & mkvar(i)) {
v[d] = i;
d++;
}
}
if (d == 0) {
fprintf(stream, "1");
} else {
for (int i = 0; i < d - 1; i++)
fprintf(stream, "%s*", variable_name[v[i]]);
fprintf(stream, "%s", variable_name[v[d - 1]]);
}
}
int monomials_product(int m1, int m2, struct monomial_ctx_t *mono)
{
return mono->monomials[j] | mono->monomials[k];
return monomial_rank(mono->monomials[m1] | mono->monomials[m2], mono);
}
void free_monomial_ctx()
int create_monomial(int degree, int *vars, struct monomial_ctx_t *mono)
{
free(mono->monomials);
free(mono);
monomial_t term = 0;
for (int i = 0; i < degree; i++)
term += mkvar(vars[i]);
int mon = monomial_rank(term, mono);
assert(hamming_weight(mono->monomials[mon]) <= 2);
assert(mon < mono->Nd[3]);
return mon;
}
int binary_search(monomial_t target, int lo, int hi)
int *get_ker(char *filename, u64 *size_ker, struct monomial_ctx_t *mono)
{
hi -= 1;
while (lo <= hi)
u64 tmp;
u64 *ker_m = load_file_64(filename, &tmp);
int *ker = (int *)malloc(sizeof(int) * tmp);
printf("ker_m[0] = %ld\n", ker_m[17]);
printf("la taille de ker est %lx\n", tmp);
for (u64 i = 0; i < tmp; i++)
ker[i] = monomial_rank(ker_m[i], mono);
free(ker_m);
*size_ker = tmp;
printf("size_ker = %ld\n", *size_ker);
printf("tmp = %ld\n", tmp);
return ker;
}
int monomial_add(int m1, int m2, struct monomial_ctx_t *mono)
{
return monomial_rank(mono->monomials[m1] ^ mono->monomials[m2], mono);
}
bool **matrix_construct(monomial_t *tab, u64 size_tab, struct monomial_ctx_t *mono)
{
bool **ans = (bool **)malloc(sizeof(bool *) * size_tab);
for (u64 i = 0; i < size_tab; i++)
{
int mid = (lo + hi) / 2;
if (target > mono->monomials[mid])
lo = mid + 1;
else if (target < mono->monomials[mid])
hi = mid - 1;
else
return mid;
ans[i] = (bool *)malloc(sizeof(bool) * NB_VEC);
for (u64 j = 0; j < NB_VEC; j++)
ans[i][j] = tab[i] & mkvar(NB_VEC - j - 1);
}
assert(0);
return (ans);
}
int monomial_rank(monomial_t m)
struct input_poly_t mult_poly_monomial(struct input_poly_t p, int index, struct monomial_ctx_t *mono)
{
int d = hamming_weight(m);
return binary_search(m, mono->Nd[d], mono->Nd[d + 1]);
struct input_poly_t res;
res.terms = (int *)malloc(sizeof(int) * p.nterms);
for (int i = 0; i < p.nterms; i++)
res.terms[i] = monomials_product(p.terms[i], index, mono);
res.nterms = p.nterms;
return (res);
}
void test_ranking()
struct input_poly_t *reconstruct(char *filename_ker, char *filename_ilog, char *filename_jlog, int D, struct input_poly_t *f, struct monomial_ctx_t *mono)
{
for (int i = 0; i < mono->Nd[mono->D + 1]; i++)
assert(monomial_rank(mono->monomials[i]) == i);
}
/* n = nb of variables */
void test_grlex_order(int n)
{
assert(hamming_weight(mono->monomials[0]) == 0);
for (int i = 0; i < n; i++)
assert(hamming_weight(mono->monomials[1 + i]) == 1);
for (int d = 0; d <= mono->D; d++)
for (int k = mono->Nd[d]; k < mono->Nd[d + 1] - 1; k++) {
assert(hamming_weight(mono->monomials[k]) == d);
assert(mono->monomials[k] < mono->monomials[k + 1]);
}
}
\ No newline at end of file
// Define variables
u64 size_mono = mono->Nd[D+1];
monomial_t *res = (monomial_t *)malloc(sizeof(monomial_t) * size_mono);
for (u64 i = 0; i < size_mono; i++)
res[i] = 0;
// Get Kernel ans ijlog
u64 size_ker;
u32 size_ilog, size_jlog;
monomial_t *ker = load_file_64(filename_ker, &size_ker);
int *ilog = load_file_32(filename_ilog, &size_ilog);
int *jlog = load_file_32(filename_jlog, &size_jlog);
if ((size_ker != size_ilog) || (size_ker != size_jlog))
err(1, "Error in i/jlog_size");
// Debut de boucle sur le ker
for (u64 i = 0; i < size_ker; i++)
{
// Calculer produit f et monome
struct input_poly_t tmp = mult_poly_monomial(f[ilog[i]], jlog[i], mono);
// on somme
for (int j = 0; j < tmp.nterms; j++)
res[tmp.terms[j]] = res[tmp.terms[j]] ^ ker[i];
free(tmp.terms);
}
// fin de boucle
free(ilog);
free(jlog);
free(ker);
// On a un tableau de monomial_t
// On passe le tableau de monomial en tableau de tableau de bool
bool **tab = matrix_construct(res, size_mono, mono);
free(res);
// On transpose
bool **tab_tr = transpose(tab, size_mono, NB_VEC);
for (u64 i = 0; i < size_mono; i++)
free(tab[i]);
free(tab);
// On get_rank le tab[Nd[D+1]-i-1]
struct input_poly_t *ans = (struct input_poly_t *)malloc(sizeof(struct input_poly_t) * NB_VEC);
for (int i = 0; i < NB_VEC; i++)
{
int cpt = 0;
for (u64 j = 0; j < size_mono; j++)
cpt += tab_tr[i][j];
ans[i].nterms = cpt;
if (cpt == 0)
ans[i].terms = (int *)malloc(sizeof(int));
else
ans[i].terms = (int *)malloc(sizeof(int) * ans[i].nterms);
cpt = 0;
for (u64 j = 0; j < size_mono; j++)
{
if (tab_tr[i][j])
{
ans[i].terms[cpt] = j;
// Maybe :
// ans[i].terms[cpt] = size_mono - 1 - j;
cpt++;
}
}
}