Commit 232dc683 authored by Charles Bouillaguet's avatar Charles Bouillaguet
Browse files

reconstruct poly

parent 69dd1117
......@@ -10,3 +10,4 @@ moebius_demo
/monica_vector
/macaulay
/macaulay_gen
/lanczos
......@@ -6,8 +6,9 @@ CFLAGS += -fopenmp
LDFLAGS += -fopenmp
LDLIBS += -lm
all: macaulay_gen reconstruct_monomials
all: macaulay_gen lanczos reconstruct_monomials
lanczos.o: CFLAGS += -flax-vector-conversions
macaulay_gen: parser.o tools.o tools_parser.o monomials.o macaulay_gen.o
macaulay_gen: LDLIBS += `pkg-config --libs m4ri`
......
#!/usr/bin/env python3
import subprocess
import sys
class Parameters:
outer_hybridation = 0
inner_hybridation = 0
degree = 2
linear_variables = 6
#
verbose = True
dry_run = True
def run(args, parameters):
if parameters.verbose:
print(' '.join(args))
if parameters.dry_run:
return
subprocess.run(args, check=True)
def macaulay_step(input_system, output_base, parameters):
args = ['./macaulay_gen', '--in', input_system, '--degree', str(parameters.degree),
'--inner-hybridation', str(parameters.linear_variables), '--out', output_base]
run(args, parameters)
D = {}
D['matrix'] = output_base + '.bin'
D['ilog'] = output_base + '.ilog.bin'
D['jlog'] = output_base + '.jlog.bin'
return D
def linear_algebra_step(Dma, output, parameters):
args = ['./lanczos', '--matrix', Dma['matrix'], '--output', output]
run(args, parameters)
return {'kernel': output}
def reconstruct_step(input_system, Dma, Dla, output, parameters):
args = ['./reconstruct_monomials', '--in', input_system, '--degree', str(parameters.degree),
'--inner-hybridation', str(parameters.linear_variables), '--out', output,
'--ilog', Dma['ilog'], '--jlog', Dma['jlog']]
run(args, parameters)
return {'polynomials': output}
input_system = "examples/random_32_quad.in"
parameters = Parameters()
Dma = macaulay_step(input_system, "matrix", parameters)
Dla = linear_algebra_step(Dma, 'kernel.bin', parameters)
Dre = reconstruct_step(input_system, Dma, Dla, "new_system.in", parameters)
\ No newline at end of file
......@@ -27,7 +27,7 @@ typedef uint32_t u32;
typedef uint16_t u16;
typedef uint8_t u8;
#define BLOCKING_FACTOR 128
#define BLOCKING_FACTOR 64
static const int n = BLOCKING_FACTOR;
char *matrix_filename;
......
/*
* This program generates a degree-D macaulay matrix from a quadratic system,
* where columns corresponding to monomials having degree at most one in certain
* variables have been removed. Finding vectors in the left-kernel of the matrix
* is the basis of the crossbred algorithm.
*/
// TODO : il faudrait tester ceci
#define _POSIX_C_SOURCE 200809L // for strdup
......@@ -6,50 +13,39 @@
#include <stdlib.h>
#include <err.h>
#include <getopt.h>
#include <inttypes.h>
#include <assert.h>
#include <string.h> // strcmp, strdup
#include <sys/time.h> // gettimeofday
#include <m4ri/m4ri.h>
#include "parser.h"
// #include "monomials.h"
// #include "tools.h"
#include "tools_parser.h"
/*
* This program generates a degree-D macaulay matrix from a quadratic system,
* where columns corresponding to monomials having degree at most one in certain
* variables have been removed. Finding vectors in the left-kernel of the matrix
* is the basis of the crossbred algorithm.
*/
// #define MAXDEG 8
struct input_poly_t * input_system;
int n; /* number of variables */
int v; /* number of non-enumerated variables (the first ones) */
int m; /* number of polynomials */
int D;
int Nd_p; /* Nd[D+1] */
int Nd_l; /* Nd[D-1] */
// int Nd[MAXDEG]; /* number of degree-d monomials */
int Nd_p; /* Nd[D+1] */
int Nd_l; /* Nd[D-1] */
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(struct monomial_ctx_t *mono)
{
assert(bad_renumbering[0] < 0);
/********************************************************************************/
for (int i = 0; i < n; i++) {
assert(bad_renumbering[i + 1] < 0);
}
test_bad_monomials_2(mono);
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);
}
void renumbering_setup(struct monomial_ctx_t *mono)
{
log(1, " - Singling out bad monomials\n");
......@@ -62,17 +58,11 @@ void renumbering_setup(struct monomial_ctx_t *mono)
nbad += bad;
}
log(2, " - %d bad monomials\n", nbad);
test_bad_monomials(mono);
// test_bad_monomials(mono);
}
struct input_poly_t * input_system;
/********************************** output buffering ********************************/
int buffer_capacity;
int buffer_size;
int * buffer;
......@@ -124,7 +114,6 @@ int buffer_write_row(int row_size, int *row)
return i;
}
/********************************** sparse matrix ***********************************/
mzd_t *M; /* dense matrix for degree D-2 multiples */
......@@ -299,15 +288,10 @@ void macaulay_sparse(int d, struct monomial_ctx_t *mono)
free(scratch);
}
/************************************* main ***********************************/
int main(int argc, char **argv)
{
// test code, not useful for the actual computation
test_shift_128();
test_HW();
/* process command-line options */
struct option longopts[5] = {
{"in", required_argument, NULL, 'i'},
......@@ -357,21 +341,21 @@ int main(int argc, char **argv)
err(1, "Cannot open %s", in_filename);
}
log(0, "Reading input system\n");
parser(input_file, NULL, &parser_setup, &parser_store_monomial, &parser_store_polynomial, NULL);
n = get_number_of_variables();
struct parser_system_t * sys = parse_polynomial_system(input_file);
n = sys->n;
m = sys->m;
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);
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, mono);
log(1, " - %d total monomials of degree <= %d\n", Nd_p, D);
renumbering_setup(mono);
input_system = convert_input_system(sys, mono);
/* prepare monomial multiples and setup sparse matrix */
macaulay_dense(D - 2, mono);
log(1, " - %d monomial multipliers\n", valid_multipliers);
......@@ -487,9 +471,6 @@ 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(mono);
write_array(rw_filename, rw, nrows);
write_array(cw_filename, cw, ncols);
write_array(ilog_filename, ilog, nrows);
......
#include <assert.h>
#include <stdlib.h>
#include <err.h>
#include "tools.h"
#include "monomials.h"
typedef __int128 u128;
typedef u128 monomial_t;
......@@ -16,13 +22,13 @@ static inline monomial_t mkvar(int i)
return ((monomial_t) 1) << i;
}
int hamming_weight(monomial_t m)
static int hamming_weight(monomial_t m)
{
return __builtin_popcountl(m) + __builtin_popcountl(m >> 64);
}
/* 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)
static int binary_search(monomial_t target, int lo, int hi, struct monomial_ctx_t *mono)
{
hi -= 1;
while (lo <= hi)
......@@ -38,32 +44,13 @@ int binary_search(monomial_t target, int lo, int hi, struct monomial_ctx_t *mono
assert(0);
}
int monomial_rank(monomial_t m, struct monomial_ctx_t *mono)
static 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)
static int monomial_enumerate(int degree, monomial_t *out, int n_var)
{
if (degree == 0) {
if (out != NULL)
......@@ -84,15 +71,13 @@ int monomial_enumerate(int degree, monomial_t *out, int n_var)
return count;
}
void monomial_setup_list(struct monomial_ctx_t *mono)
static void monomial_setup_list(struct monomial_ctx_t *mono)
{
int D = mono->D;
assert(D >= 2);
double start = wtime();
log(0, "Monomials setup\n");
log(1, " - Counting monomials\n");
mono->Nd[0] = 0;
......@@ -106,11 +91,13 @@ void monomial_setup_list(struct monomial_ctx_t *mono)
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]]), mono->n_var);
test_grlex_order(mono);
test_ranking(mono);
// test_grlex_order(mono);
// test_ranking(mono);
log(1, " - Monomials setup: %.1fs\n", wtime() - start);
}
/* public functions */
struct monomial_ctx_t *monomial_setup(int n_var, int v, int D)
{
struct monomial_ctx_t *mono = (struct monomial_ctx_t *)malloc(sizeof(struct monomial_ctx_t));
......@@ -135,35 +122,6 @@ bool monomial_is_bad(int index, struct monomial_ctx_t *mono)
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, struct monomial_ctx_t *mono)
{
if (d <= mono->D+1)
......@@ -171,9 +129,9 @@ int get_Nd(int d, struct monomial_ctx_t *mono)
errx(1, "deg in get_Nd too high");
}
void print_monomial(FILE *stream, int ind, char **variable_name, struct monomial_ctx_t *mono)
void print_monomial(FILE *stream, int monomial, char **variable_name, struct monomial_ctx_t *mono)
{
monomial_t m = mono->monomials[ind];
monomial_t m = mono->monomials[monomial];
int v[mono->n_var];
int d = 0;
for (int i = 0; i < mono->n_var; i++) {
......@@ -196,8 +154,6 @@ int monomials_product(int m1, int m2, struct monomial_ctx_t *mono)
return monomial_rank(mono->monomials[m1] | mono->monomials[m2], mono);
}
int create_monomial(int degree, int *vars, struct monomial_ctx_t *mono)
{
monomial_t term = 0;
......@@ -207,149 +163,4 @@ int create_monomial(int degree, int *vars, struct monomial_ctx_t *mono)
assert(hamming_weight(mono->monomials[mon]) <= 2);
assert(mon < mono->Nd[3]);
return mon;
}
int *get_ker(char *filename, u64 *size_ker, struct monomial_ctx_t *mono)
{
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++)
{
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);
}
return (ans);
}
struct input_poly_t mult_poly_monomial(struct input_poly_t p, int index, struct monomial_ctx_t *mono)
{
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);
}
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)
{
// 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++;
}
}
}
for (int i = 0; i < NB_VEC; i++)
free(tab_tr[i]);
free(tab_tr);
return (ans);
// On renvoie la liste des rank
}
// int **monomial_transpose(int *mat, u64 size_ker, u64 size_mono, struct monomial_ctx_t *mono)
// {
// int **ans = (int **)malloc(sizeof(int *) * 2 * size_ker);
// // size_ker = 64
// for (u64 i = 0; i < size_ker; i++)
// {
// ans[2*i] = (int *)malloc(sizeof(int));
// ans[2*i][0] = 0;
// // size_mono = Nd[D+1]
// for (u64 j = 0; j < size_mono; j++)
// {
// if (mkvar(i+1) & mono->monomials[mat[j]])
// ans[2*i][0]++;
// }
// ans[2*i+1] = (int *)malloc(sizeof(int) * ans[2*i][0]);
// int cpt = 0;
// for (u64 j = 0; j < size_mono; j++)
// {
// if (mkvar(i+1) & mono->monomials[mat[j]])
// {
// ans[2*i+1][cpt] = mat[j];
// cpt++;
// }
// }
// }
// return (ans);
// }
}
\ No newline at end of file
#include "tools.h"
typedef __int128 u128;
typedef u128 monomial_t;
struct monomial_ctx_t;
#define MAXDEG 8
struct monomial_ctx_t *monomial_setup(int n_var, int v, int D);
struct monomial_ctx_t;
struct monomial_ctx_t * monomial_setup(int n_var, int v, int D);
void free_monomial_ctx(struct monomial_ctx_t *mono);
bool monomial_is_bad(int index, struct monomial_ctx_t *mono);
/* all tests */
void test_bad_monomials_2(struct monomial_ctx_t *mono);
void test_shift_128();
void test_HW();
int get_Nd(int d, struct monomial_ctx_t *mono);
void print_monomial(FILE *stream, int ind, char **variable_name, struct monomial_ctx_t *mono);
/* macaulay part */
int monomials_product(int m1, int m2, struct monomial_ctx_t *mono);
int create_monomial(int degree, int *vars, struct monomial_ctx_t *mono);
bool monomial_is_bad(int index, struct monomial_ctx_t *mono);
/* reconstruct part */
int *get_ker(char *filename, u64 *size_ker, struct monomial_ctx_t *mono);
int monomial_add(int m1, int m2, struct monomial_ctx_t *mono);
// int *get_ker(char *filename, u64 *size_ker, struct monomial_ctx_t *mono);
// int monomial_add(int m1, int m2, struct monomial_ctx_t *mono);
// int **monomial_transpose(int *mat, u64 size_ker, u64 size_mono, struct monomial_ctx_t *mono);
bool **matrix_construct(monomial_t *tab, u64 size_tab, struct monomial_ctx_t *mono);
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);
// bool **matrix_construct(monomial_t *tab, u64 size_tab, struct monomial_ctx_t *mono);
// 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);
//
\ No newline at end of file
CC = gcc
CFLAGS = -std=c11 -g -Wall -Wextra -O3 -Wno-unused-parameter -Wno-maybe-uninitialized -Werror -march=native
# OpenMP
CFLAGS += -fopenmp
LDFLAGS += -fopenmp
LDLIBS += -lm
all: macaulay_gen
macaulay_gen: parser.o tools.o macaulay_gen.o
macaulay_gen: LDLIBS += `pkg-config --libs m4ri`
parser.o: parser.h
tools.o: tools.h
macaulay_gen.o: parser.h tools.h
macaulay_gen.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
#!/bin/bash
echo "Start script"
rm -rf my_mat/ ker.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 ker.bin
echo "End script"
\ No newline at end of file
a, b, c
a
a*b
a*c
b
b*c
c