Commit 66179232 authored by Quentin HAMMERER's avatar Quentin HAMMERER
Browse files

J'ai du virer une option gnu=c11 dans le makefile pour pouvoir compiler avec...

J'ai du virer une option gnu=c11 dans le makefile pour pouvoir compiler avec strdup, toujours modifier le chemin vers bwc.pl dans le ger_ker.sh, ou simplement utiliser lanczos si il est au point
parent 91261ba5
CC = gcc
CFLAGS = -std=c11 -g -Wall -Wextra -O3 -Wno-unused-parameter -Wno-maybe-uninitialized -Werror -march=native -mavx2
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
all: parser_demo parser_demo_alt moebius_demo monica monica_vector macaulay_gen lanczos 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 macaulay_gen.o
macaulay_gen: parser.o tools.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: 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
macaulay_gen.o: 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`
clean:
rm -f *.o
rm -f parser_demo moebius_demo parser_demo_alt monica monica_vector macaulay_gen lanczos
rm -rf *.o my_mat/ *log.bin
rm -f parser_demo moebius_demo parser_demo_alt monica monica_vector macaulay_gen lanczos reconstruct_monomials
#!/bin/bash
echo "Start script"
rm -rf my_mat/ kernel.bin
mkdir my_mat
mv mat.bin my_mat/
mv mat.cw.bin my_mat/
mv mat.rw.bin my_mat/
~/stage/cado-nfs-master/build/poste-almasty-02/linalg/bwc/./bwc.pl :complete balancing_options="reorder=columns" matrix=`pwd`/my_mat/mat.bin mn=64
cp my_mat/mat.bin-1x1/W kernel.bin
echo "End script"
\ No newline at end of file
......@@ -14,12 +14,8 @@
#include <m4ri/m4ri.h>
#include "parser.h"
typedef __int128 u128;
typedef uint64_t u64;
typedef uint32_t u32;
typedef uint16_t u16;
#include "monomials.h"
#include "tools.h"
/*
* This program generates a degree-D macaulay matrix from a quadratic system,
......@@ -28,382 +24,71 @@ typedef uint16_t u16;
* is the basis of the crossbred algorithm.
*/
#define MAXDEG 8
// #define MAXDEG 8
int n; /* number of variables */
int v; /* number of non-enumerated variables (the first ones) */
int m; /* number of polynomials */
int D;
int Nd[MAXDEG]; /* number of degree-d monomials */
int Nd_p; /* Nd[D+1] */
int Nd_l; /* Nd[D-1] */
// int Nd[MAXDEG]; /* number of degree-d monomials */
int log_level = 99;
#define log(level, fmt, ...) \
do { if (level < log_level) fprintf(stderr, fmt, ##__VA_ARGS__); } while (0)
double wtime()
{
struct timeval ts;
gettimeofday(&ts, NULL);
return (double)ts.tv_sec + ts.tv_usec / 1e6;
}
/* represent n in 4 bytes */
void human_format(u64 n, char *target) {
if (n < 1000) {
sprintf(target, "%" PRId64, n);
return;
}
if (n < 1000000) {
sprintf(target, "%.1fK", n / 1e3);
return;
}
if (n < 1000000000) {
sprintf(target, "%.1fM", n / 1e6);
return;
}
if (n < 1000000000000ll) {
sprintf(target, "%.1fG", n / 1e9);
return;
}
if (n < 1000000000000000ll) {
sprintf(target, "%.1fT", n / 1e12);
return;
}
}
void * alloc_array(u64 n_entries, u64 entry_size, const char *role)
{
void *A = malloc(entry_size * n_entries);
if (A == NULL)
err(1, "cannot allocate %s", role);
return A;
}
void write_array(const char *filename, const int *A, int size)
{
log(1, " - writing %s\n", filename);
FILE *file = fopen(filename, "w");
if (file == NULL)
errx(1, "cannot open %s", filename);
size_t check = fwrite(A, sizeof(*A), size, file);
if (check != (size_t) size)
err(1, "error writing %s", filename);
fclose(file);
}
/*********************** parser interface ************************/
struct parser_monomial_t {
int degree;
int vars[2];
struct parser_monomial_t * next;
};
struct parser_poly_t {
struct parser_monomial_t * terms_head;
struct parser_monomial_t * terms_tail;
struct parser_poly_t * next;
};
char ** variable_name;
struct parser_poly_t * poly_head;
struct parser_poly_t * poly_tail;
struct parser_monomial_t * new_monomial()
{
struct parser_monomial_t *m = malloc(sizeof(*m));
m->degree = 0;
m->next = NULL;
return m;
}
struct parser_poly_t * new_polynomial()
{
struct parser_poly_t * p = malloc(sizeof(*p));
p->next = NULL;
// allocate dummy monomial
struct parser_monomial_t *dummy = new_monomial();
p->terms_head = dummy;
p->terms_tail = dummy;
return p;
}
bool polynomial_started = false;
/* callbacks */
void parser_setup(void *opaque, int nvar, const char **vars)
{
n = nvar;
variable_name = malloc(n * sizeof(char *));
for (int i = 0; i < n; i++)
variable_name[i] = strdup(vars[i]);
// allocate dummy polynomial
struct parser_poly_t *dummy = new_polynomial();
poly_head = dummy;
poly_tail = dummy;
polynomial_started = false;
}
void parser_store_monomial(void *opaque, int line, int column, int degree, const int *variables)
{
if (degree > 2)
errx(1, "This program only supports quadratic polynomials (degree-%d term on line %d)", degree, line);
if (!polynomial_started) {
struct parser_poly_t *n = new_polynomial();
poly_tail->next = n;
poly_tail = n;
polynomial_started = true;
}
struct parser_poly_t *p = poly_tail;
// build the monomial
struct parser_monomial_t * m = new_monomial();
m->degree = degree;
for (int i = 0; i < degree; i++)
m->vars[i] = variables[i];
// store the monomial
p->terms_tail->next = m;
p->terms_tail = m;
}
void parser_store_polynomial(void *opaque, int line)
{
polynomial_started = false;
m += 1;
}
/*********************** monomial numbering *************************/
/*
* Monomials are numbered from 0 to sum(binomial(n, i), i=0..D). The main operation
* we have to provide is multiplication by a single variable (for all monomials except
* those of the greatest degree).
*
* Degree-d monomials are in monomials[Nd[d]:Nd[d + 1]].
* The number of monomials of degree <= d is Nd[d + 1].
* This is NOT a monomial order : a<b but a*b > b*b
*
* "bad" monomial contain more than one special variable.
* bad_renumbering[j] == -1 if j is a good monomial
* bad_renumbering[j] == i (>= 0) if j is the i-th bad monomial
*/
typedef u128 monomial_t;
monomial_t *monomials; /* array to rank/unrank 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 */
static inline monomial_t mkvar(int i)
{
return ((monomial_t) 1) << i;
}
void print_monomial(FILE *stream, monomial_t m)
{
int v[n];
int d = 0;
for (int i = 0; i < n; 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]]);
}
}
/* 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)
{
hi -= 1;
while (lo <= hi) {
int mid = (lo + hi) / 2;
if (target > monomials[mid])
lo = mid + 1;
else if (target < monomials[mid])
hi = mid - 1;
else
return mid;
}
assert(0);
}
int hamming_weight(monomial_t m)
{
return __builtin_popcountl(m) + __builtin_popcountl(m >> 64);
}
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 monomial_rank(monomial_t m)
{
int d = hamming_weight(m);
return binary_search(m, Nd[d], Nd[d + 1]);
}
void test_ranking()
{
for (int i = 0; i < Nd[D + 1]; i++)
assert(monomial_rank(monomials[i]) == i);
}
/* Variants of this code enumerate monomials of degree <= D in lex order */
int monomial_enumerate(int degree, monomial_t *out)
{
if (degree == 0) {
if (out != NULL)
out[0] = 0;
return 1;
}
int count = 0;
monomial_t top = mkvar(n);
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 test_grlex_order()
{
assert(hamming_weight(monomials[0]) == 0);
for (int i = 0; i < n; i++)
assert(hamming_weight(monomials[1 + i]) == 1);
for (int d = 0; d <= D; d++)
for (int k = Nd[d]; k < Nd[d + 1] - 1; k++) {
assert(hamming_weight(monomials[k]) == d);
assert(monomials[k] < monomials[k + 1]);
}
}
/* returns true when m contains more than one special variable */
static inline bool monomial_is_bad(monomial_t m)
{
monomial_t inter = m & vmask;
return (inter & (inter - 1)) != 0;
}
void test_bad_monomials()
{
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++)
assert(monomial_is_bad(mkvar(i) | mkvar(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++)
assert(!monomial_is_bad(mkvar(i) | mkvar(j)));
{
monomial_t inter = (mkvar(i) | mkvar(j)) & vmask;
assert((inter & (inter - 1)) == 0);
}
}
}
void monomial_setup_list()
void renumbering_setup()
{
assert(D >= 2);
double start = wtime();
vmask = (v == 0) ? 0 : mkvar(v) - 1;
log(0, "Monomials setup\n");
log(1, " - Counting monomials\n");
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... */
Nd[d + 1] = Nd[d] + binom_n_d;
log(2, " - %d monomials of degree %d\n", binom_n_d, d);
}
log(1, " - Building monomial list\n");
monomials = alloc_array(Nd[D + 1], sizeof(*monomials), "monomial list");
for (int d = 0; d <= D; d++)
monomial_enumerate(d, &monomials[Nd[d]]);
test_grlex_order();
test_ranking();
log(1, " - Singling out bad monomials\n");
bad_renumbering = alloc_array(Nd[D + 1], sizeof(*bad_renumbering), "bad monomial renumbering");
inv_bad_renumbering = alloc_array(Nd[D + 1], sizeof(*bad_renumbering), "bad monomial renumbering");
for (int j = 0; j < Nd[D + 1]; j++) {
bool bad = monomial_is_bad(monomials[j]);
bad_renumbering[j] = bad ? nbad : -1;
inv_bad_renumbering[nbad] = j;
nbad += bad;
}
log(2, " - %d bad monomials\n", nbad);
test_bad_monomials();
log(1, " - Monomials setup: %.1fs\n", wtime() - start);
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);
bad_renumbering[j] = bad ? nbad : -1;
inv_bad_renumbering[nbad] = j;
nbad += bad;
}
log(2, " - %d bad monomials\n", nbad);
test_bad_monomials();
}
/******************************* input system in nicer form *****************************/
struct input_poly_t {
int nterms;
int terms[8257]; /* ok for 128 variables */
};
struct input_poly_t * input_system;
void convert_input_system()
{
input_system = malloc(m * sizeof(*input_system));
if (input_system == NULL)
errx(1, "cannot allocate input system");
int i = 0;
struct parser_poly_t *head = poly_head->next;
for (struct parser_poly_t * p = head; p != NULL; p = p->next) {
int size = 0;
struct parser_monomial_t *head = p->terms_head->next;
for (struct parser_monomial_t *m = head; m != NULL; m = m->next) {
if (size == 8257)
errx(1, "stupid polynomial is too long (shouldn't happen with <= 128 variables)");
monomial_t term = 0;
for (int i = 0; i < m->degree; i++)
term += mkvar(m->vars[i]);
int mon = monomial_rank(term);
assert(hamming_weight(monomials[mon]) <= 2);
assert(mon < Nd[3]);
input_system[i].terms[size] = mon;
size += 1;
}
input_system[i].nterms = size;
i += 1;
}
}
/********************************** output buffering ********************************/
......@@ -554,20 +239,23 @@ static inline int mzd_ffs(mzd_t *M, rci_t row)
void macaulay_dense(int d)
{
valid_multipliers = Nd[d + 1];
int Nd_p = get_Nd(d+1);
int Nd_l = get_Nd(d-1);
valid_multipliers = Nd_p;
if (d < 2)
return;
int Nd_3 = get_Nd(3);
log(0, "Generating degree-%d (dense) Macaulay matrix \n", d);
log(1, " - matrix size: %d x %d\n", m * Nd[d - 1], Nd[d + 1]);
M = mzd_init(m * Nd[d - 1], Nd[d + 1]);
LM = alloc_array(Nd[d + 1], sizeof(*LM), "leading monomials");
log(1, " - matrix size: %d x %d\n", m * Nd_l, Nd_p);
M = mzd_init(m * Nd_l, Nd_p);
LM = alloc_array(Nd_p, sizeof(*LM), "leading monomials");
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++) {
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[j] | monomials[k];
for (int k = 0; k < Nd_3; k++) {
monomial_t product = monomials_product(j, k);
mmul[k] = monomial_rank(product);
}
for (int i = 0; i < m; i++) {
......@@ -584,7 +272,7 @@ void macaulay_dense(int d)
int rank = mzd_echelonize(M, 0);
if (rank < M->nrows)
errx(1, "Bizarre, rank defect!");
for (int i = 0; i < Nd[d + 1]; i++)
for (int i = 0; i < Nd_p; i++)
LM[i] = 0;
for (int i = 0; i < rank; i++) {
int j = M->ncols - 1 - mzd_ffs(M, i);
......@@ -605,8 +293,11 @@ void macaulay_sparse(int d)
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);
#pragma omp for
for (int j = Nd[d]; j < Nd[d + 1]; j++) {
for (int j = Nd; j < Nd_p; j++) {
if (LM && LM[j]) {
#pragma omp atomic update
skipped += 1;
......@@ -614,8 +305,8 @@ 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[j] | monomials[k];
for (int k = 0; k < Nd_3; k++) {
monomial_t product = monomials_product(j, k);
mmul[k] = monomial_rank(product);
}
......@@ -690,13 +381,19 @@ int main(int argc, char **argv)
}
log(0, "Reading input system\n");
parser(input_file, NULL, &parser_setup, &parser_store_monomial, &parser_store_polynomial, NULL);
n = get_number_of_variables();
log(1, " - Read %d polynomials in %d variables\n", m, n);
log(1, " - Generating degree-%d Macaulay matrix\n", D);
log(1, " - Truncating columns linear in the first %d variables\n", v);
monomial_setup_list();
convert_input_system();
log(1, " - %d total monomials of degree <= %d\n", Nd[D+1], D);
monomial_setup(v, D);
Nd_p = get_Nd(D+1);
Nd_l = get_Nd(D-1);
renumbering_setup();
input_system = convert_input_system(Nd_p, &m);
log(1, " - %d total monomials of degree <= %d\n", Nd_p, D);
/* prepare monomial multiples and setup sparse matrix */
macaulay_dense(D - 2);
......@@ -750,7 +447,7 @@ int main(int argc, char **argv)
ncols = 0;
u64 nnz_check = 0;
/* check empty columns ; enumerate all monomials again */
for (int j = 0; j < Nd[D + 1]; j++) {
for (int j = 0; j < Nd_p; j++) {
int jj = bad_renumbering[j];
assert(jj < 0 || ncols <= jj);
if (jj < 0)
......@@ -813,6 +510,9 @@ int main(int argc, char **argv)
log(0, "Finalization\n");
flush_buffer();
fclose(buffer_stream);
free(inv_bad_renumbering);
free(bad_renumbering);
free_monomial_ctx();
write_array(rw_filename, rw, nrows);
write_array(cw_filename, cw, ncols);
write_array(ilog_filename, ilog, nrows);
......
#include "monomials.h"
#include "tools.h"
struct monomial_ctx_t{
int v;
int D;
int Nd[MAXDEG];
monomial_t *monomials;
};
struct monomial_ctx_t *mono;
bool monomial_is_bad(int index, monomial_t vmask)
{
monomial_t m = mono->monomials[index];
monomial_t inter = m & vmask;
return (inter & (inter - 1)) != 0;
}
void monomial_setup_list()
{
int D = mono->D;