Commit 6ccc0672 authored by QuentinH's avatar QuentinH
Browse files

a

parent 2b12d3b3
......@@ -6,7 +6,7 @@ CFLAGS += -fopenmp
LDFLAGS += -fopenmp
LDLIBS += -lm
all: macaulay_gen lanczos reconstruct_monomials
all: macaulay_gen lanczos reconstruct_monomials test_parser_main test_polynomial_main
lanczos.o: CFLAGS += -flax-vector-conversions
monica: parser.o monica.o
......@@ -19,13 +19,19 @@ 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 tools_parser.o monomials.o reconstruct_monomials.o
reconstruct_monomials: LDLIBS += `pkg-config --libs m4ri`
monomial_new_test: tools.o monomial_new.o monomial_new_test.o
test_parser_main: tools_parser.o monomials.o tools.o parser.o test_parser_main.o
test_polynomial_main: tools_parser.o monomials.o tools.o parser.o polynomial.o test_polynomial_main.o
#test: parser.o tools.o tools_parser.o monomials.o test.o
#test: LDLIBS += `pkg-config --libs m4ri`
parser.o: parser.h
tools.o: tools.h
monomials.o: monomials.h tools.h
polynomial.o: polynomial.h tools.h tools_parser.h monomials.h
tools_parser.o: tools_parser.h monomials.h tools.h
monomial_new.o: tools.h monomial_new.h
macaulay_gen.o: parser.h tools_parser.h tools.h monomials.h
macaulay_gen.o: CFLAGS += `pkg-config --cflags m4ri`
......@@ -33,9 +39,16 @@ macaulay_gen.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`
monomial_new_test.o:tools.h monomial_new.h
test_parser_main.o: tools_parser.h monomials.h tools.h parser.h
test_polynomial_main.o: tools_parser.h monomials.h tools.h parser.h polynomial.h
#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 test
rm -f parser_demo moebius_demo parser_demo_alt monica monica_vector macaulay_gen lanczos reconstruct_monomials test test_parser_main
File added
......@@ -3,30 +3,25 @@
#include "tools_parser.h"
#include "tools.h"
// #include "monomials.h"
#include "monomial_new.h"
// int convert_gray_code(int n)
// {
// return n ^ (n >> 1);
// }
struct polynomial_system {
int n; /* # variables */
int v; /* # linear variables */
char ** variable_name; /* array of strings, size n */
int D; /* max degree of any monomial */
int m; /* # polynomials */
bool ** poly; /* array of polynomials, size m */
struct monomial_ctx_t * mono; /* description of the monomials */
};
#include "enum.h"
struct matrix_system {
int row;
int col;
int n; /* # of monomials */
bool **A; /* A[row*col][n] */
bool **B; /* B[row*1][n] */
struct monomial_ctx_t *mono;
};
// This function still need to be tested
// Give A matrix values with monomial sol input
bool **update_mat(bool **poly, monomial_t sol, struct polynomial_system *poly_sys)
bool **update_mat_A(bool **poly, monomial_t sol, struct polynomial_system *poly_sys)
{
bool **res;
// Aloccate res
// Alloccate res
res = malloc(sizeof(bool *)* poly_sys->m);
for (int i = 0; i < poly_sys->m; i++)
res[i] = malloc(sizeof(bool) * poly_sys->v);
......@@ -55,13 +50,67 @@ bool **update_mat(bool **poly, monomial_t sol, struct polynomial_system *poly_sy
return res;
}
//TODO
monomial_t solve(monomial_t sol, struct polynomial_system *poly_sys)
bool *update_mat_B(bool **poly, monomial_t sol, struct polynomial_system *poly_sys)
{
bool *res;
res = malloc(sizeof(bool)* poly_sys->m);
for (int i = 0; i < poly_sys->m; i++)
{
// check if there is 1 in polynome
bool cur = poly[i][0];
for (int j = 1; j < monomial_rank_degD(poly_sys->n, poly_sys->mono); j++)
{
if (poly[i][j] == 1)
{
if (monomial_gcd(sol, monomial_unrank(poly_sys->mono, j) == sol))
cur ^= 1;
}
}
res[i] = cur;
}
return res;
}
monomial_t bool_to_mono(bool *arr, int length)
{
monomial_t ans = 0;
for (int i = 0; i < length; i++)
{
ans <<= 1;
ans ^= arr[i];
}
return ans;
}
monomial_t test_sol(monomial_t sol, struct polynomial_system *poly_sys)
{
if (sol == -1)
return -1;
return;
bool **new_A = update_mat_A(poly_sys->A, sol, poly_sys);
bool *new_B = update_mat_B(poly_sys->B, sol, poly_sys);
// solve from pdg.c
bool *new_sol = solve(new_A, new_B, poly_sys->m, poly_sys->v);
bool tmp = check_ans(new_A, new_B, poly_sys->m, poly_sys->v, new_sol);
// Free A and B
for (int i = 0; i < poly_sys->m; i++)
free(new_A[i]);
free(new_A);
free(new_B);
// bad solution
if (tmp == 0)
return -1;
return bool_to_mono(tmp, poly_sys->v);
}
monomial_t enum_recur(monomial_t a, int b, int v, struct polynomial_system *poly_sys)
......@@ -73,9 +122,9 @@ monomial_t enum_recur(monomial_t a, int b, int v, struct polynomial_system *poly
monomial_t sol1 = enum_recur(a, 0, v-1, poly_sys);
monomial_t sol2 = enum_recur(a, 1, v-1, poly_sys);
monomial_t tmp = solve(sol1, poly_sys);
monomial_t tmp = test_sol(sol1, poly_sys);
if (tmp == -1)
tmp = solve(sol2, poly_sys);
tmp = test_sol(sol2, poly_sys);
return tmp;
}
......@@ -96,10 +145,6 @@ monomial_t enumerate(struct polynomial_system *poly_sys)
int main()
{
// for (int n = 0; n < 16; n++)
// {
// printf("n = %d, ng = %d\n", n, convert_gray_code(n));
// }
printf("%d\n", 1 << 1);
return 0;
}
\ No newline at end of file
#include "monomial_new.h"
struct polynomial_system {
int n; /* # variables */
int v; /* # linear variables */
char ** variable_name; /* array of strings, size n */
int D; /* max degree of any monomial */
int m; /* # polynomials */
bool ** poly; /* array of polynomials, size m */
bool **A;
bool **B;
struct monomial_ctx_t * mono; /* description of the monomials */
};
bool **update_mat(bool **poly, monomial_t sol, struct polynomial_system *poly_sys);
monomial_t enumerate(struct polynomial_system *poly_sys);
\ No newline at end of file
This diff is collapsed.
#include "monomial_new.h"
#include <string.h> // ffs
struct monomial_ctx_t
{
monomial_t *monomials;
int Nd[MAXDEG];
int D;
int n_var;
int v;
};
monomial_t monomial_one()
......@@ -15,7 +8,7 @@ monomial_t monomial_one()
return (monomial_t)0;
}
monomial_t monomial_var(int n)
monomial_t mkvar(int n)
{
monomial_t res = 1;
res = res << n;
......@@ -37,22 +30,21 @@ int monomial_degree(monomial_t a)
return hamming_weight(a);
}
void monomial_deconstruct(monomial_t a, int *vars, int nb_vars)
{
vars = malloc(sizeof(int) * nb_vars);
for (int i = 0; i < nb_vars; i++)
{
vars[i] = a%2;
a = a >> 1;
}
return ;
}
// void monomial_deconstruct(monomial_t a, int *vars, int nb_vars)
// {
// vars = malloc(sizeof(int) * nb_vars);
// for (int i = 0; i < nb_vars; i++)
// {
// vars[i] = a%2;
// a = a >> 1;
// }
// return ;
// }
void monomial_deconstruct(monomial_t a, int *vars, int nb_vars)
{
vars = malloc(sizeof(int) * nb_vars);
for (int i = 0; i < nb_vars; i++)
vars[nb_vars - i - 1] = a&(1<<i);
vars[nb_vars - i - 1] = (int)(((a&(1<<i))>>i));
return ;
}
......@@ -75,7 +67,7 @@ int monomial_getvar(monomial_t a)
return -1;
}
static int hamming_weight(monomial_t m)
int hamming_weight(monomial_t m)
{
return __builtin_popcountl(m) + __builtin_popcountl(m >> 64);
}
......@@ -96,7 +88,7 @@ static int binary_search(monomial_t target, int lo, int hi, struct monomial_ctx_
assert(0);
}
bool monomial_isvalid(const struct monomial_ctx_t *mono, monomial_t a)
bool monomial_isvalid(struct monomial_ctx_t *mono, monomial_t a)
{
monomial_t mask = 0;
for (int i = 1; i < mono->Nd[1]; i++)
......@@ -106,13 +98,13 @@ bool monomial_isvalid(const struct monomial_ctx_t *mono, monomial_t a)
return 0;
}
int monomial_rank(const struct monomial_ctx_t *mono, monomial_t a)
int monomial_rank(struct monomial_ctx_t *mono, monomial_t a)
{
int d = hamming_weight(a);
return binary_search(a, mono->Nd[d], mono->Nd[d + 1], mono);
}
monomial_t monomial_unrank(const struct monomial_ctx_t *mono, int n)
monomial_t monomial_unrank(struct monomial_ctx_t *mono, int n)
{
return mono->monomials[n];
}
......@@ -123,6 +115,29 @@ int monomial_rank_degD(int D, struct monomial_ctx_t *mono)
}
// enumerate in grlex order
static 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;
}
// enumerate in lex order
static int monomial_enumerate_lex(int degree, monomial_t *out, int n_var)
{
if (degree == 0) {
......@@ -153,9 +168,9 @@ static int monomial_enumerate_lex(int degree, monomial_t *out, int n_var)
}
static int monomial_rank_lex(monomial_t m, struct monomial_ctx_t *mono)
int monomial_rank_lex(monomial_t m, struct monomial_ctx_t *mono)
{
return binary_search(m, mono->Nd[0], mono->Nd[1], mono);
return binary_search(m, mono->Nd[0], mono->Nd[mono->n_var], mono);
}
......@@ -196,6 +211,21 @@ struct monomial_ctx_t *monomial_setup(int n_var, int v, int D)
return mono;
}
void test_lex(struct monomial_ctx_t *mono)
{
monomial_t tmp = -1;
for (int i = 0; i < mono->Nd[mono->n_var]; i++)
{
if (monomial_rank_lex(mono->monomials[i], mono) != i)
printf("Error in lex ranking\n");
if (hamming_weight(mono->monomials[i]) > 3)
printf("Error in hamming_weight lex\n");
if (mono->monomials[i] <= tmp)
printf("Error not in lexicographical order !\n");
tmp = mono->monomials[i];
}
// printf("test_lex :\nEZ\n");
}
static void monomial_setup_list_lex(struct monomial_ctx_t *mono)
{
......@@ -218,10 +248,9 @@ static void monomial_setup_list_lex(struct monomial_ctx_t *mono)
mono->monomials = alloc_array(mono->Nd[mono->n_var], sizeof(*(mono->monomials)), "monomial list");
mono->monomials[0] = 0;
for (int i = 1; i < mono->n_var; i++)
monomial_enumerate(D, mono->monomials + mono->Nd[i-1], i);
// test_grlex_order(mono);
// test_ranking(mono);
for (int i = 1; i <= mono->n_var; i++)
monomial_enumerate_lex(D, mono->monomials + mono->Nd[i-1], i);
test_lex(mono);
log(1, " - Monomials setup: %.1fs\n", wtime() - start);
}
......@@ -236,3 +265,13 @@ struct monomial_ctx_t *monomial_setup_lex(int n_var, int v, int D)
monomial_setup_list_lex(mono);
return mono;
}
void print_monomial(monomial_t m, struct monomial_ctx_t *mono)
{
int *vars = malloc(sizeof(int) * mono->n_var);
monomial_deconstruct(m, vars, mono->n_var);
for (int i = 0; i < mono->n_var; i++)
printf("%d", vars[i]);
printf("\n");
free(vars);
}
......@@ -6,7 +6,16 @@
#include <err.h>
typedef __int128 monomial_t;
struct monomial_ctx_t;
// struct monomial_ctx_t;
// For testing purpose, this is here
struct monomial_ctx_t
{
monomial_t *monomials;
int Nd[MAXDEG];
int D;
int n_var;
int v;
};
/* generic monomials */
......@@ -23,11 +32,17 @@ int monomial_getvar(monomial_t a);
/* ordered monomials */
struct monomial_ctx_t * monomial_setup_lex(int n_var, int D);
struct monomial_ctx_t * monomial_setup_lex(int n_var, int v, int D);
struct monomial_ctx_t * monomial_setup(int n_var, int v, int D);
bool monomial_isvalid(const struct monomial_ctx_t *mono, monomial_t a);
int monomial_rank(const struct monomial_ctx_t *mono, monomial_t a);
monomial_t monomial_unrank(const struct monomial_ctx_t *mono, int i);
int monomial_rank(struct monomial_ctx_t *mono, monomial_t a);
monomial_t monomial_unrank(struct monomial_ctx_t *mono, int n);
int monomial_rank_degD(int D, struct monomial_ctx_t *mono);
bool monomial_isvalid(struct monomial_ctx_t *mono, monomial_t a);
int monomial_rank_lex(monomial_t m, struct monomial_ctx_t *mono);
int hamming_weight(monomial_t m);
void print_monomial(monomial_t m, struct monomial_ctx_t *mono);
\ No newline at end of file
#include "monomial_new.h"
// void test_lex(struct monomial_ctx_t *mono)
// {
// monomial_t tmp = -1;
// for (int i = 0; i < mono->Nd[mono->n_var]; i++)
// {
// if (monomial_rank_lex(mono->monomials[i], mono) != i)
// printf("Erreur dans test_lex dans le rank\n");
// if (hamming_weight(mono->monomials[i]) > 3)
// printf("Erreur dans le hamming_weight\n");
// if (mono->monomials[i] <= tmp)
// printf("Erreur pas dans l'ordre lexicographique !\n");
// tmp = mono->monomials[i];
// }
// printf("test_lex :\nEZ\n");
// }
int main()
{
struct monomial_ctx_t *mono = monomial_setup_lex(6, 2, 3);
free(mono->monomials);
// test_lex(mono);
// for (int i = 0; i < mono->Nd[mono->n_var]; i++)
// print_monomial(mono->monomials[i], mono);
return (0);
}
\ No newline at end of file
......@@ -17,14 +17,6 @@ struct monomial_ctx_t
int v;
};
struct polynomial_system {
int n; /* # variables */
char ** variable_name; /* array of strings, size n */
int D; /* max degree of any monomial */
int m; /* # polynomials */
bool ** poly; /* array of polynomials, size m */
struct monomial_ctx * mono; /* description of the monomials */
};
static inline monomial_t mkvar(int i)
{
......@@ -133,9 +125,10 @@ bool monomial_is_bad(int index, struct monomial_ctx_t *mono)
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");
return mono->Nd[d];
errx(1, "deg in get_Nd too high");
}
void print_monomial(FILE *stream, int monomial, char **variable_name, struct monomial_ctx_t *mono)
......@@ -192,71 +185,186 @@ int create_monomial(int degree, int *vars, struct monomial_ctx_t *mono)
int *known_lin(struct polynomial_system *poly_sys)
// int *known_lin(struct polynomial_system *poly_sys)
// {
// int *res = (int *)malloc(sizeof(int) * poly_sys->mono->Nd[D+1]);
// monomial_t mask = 2**(poly_sys->v + 1) - 1;
// for (int i = 0; i < poly_sys->mono->Nd[D+1]; i++)
// {
// monomial_t tmp = poly_sys->mono->monomials[i];
// tmp &= mask;
// // non linear
// if (tmp == 0)
// res[i] = -1;
// else
// res[i] = ffs(tmp);
// }
// return res;
// }
// enumerate in lex order
static int monomial_enumerate_lex(int degree, monomial_t *out, int n_var)
{
int *res = (int *)malloc(sizeof(int) * poly_sys->mono->Nd[D+1]);
monomial_t mask = 2**(poly_sys->v + 1) - 1;
for (int i = 0; i < poly_sys->mono->Nd[D+1]; i++)
if (degree == 0) {
if (out != NULL)
out[0] = 0;
return 1;
}
int count = 0;
monomial_t top = mkvar(n_var);
monomial_t i = mkvar(n_var-1);
while (i < top) {
if (out != NULL)
out[count] = i;
count++;
if (hamming_weight(i) < degree)
i++;
else if (hamming_weight(i) == degree)
{
monomial_t tmp = poly_sys->mono->monomials[i];
tmp &= mask;
// non linear
if (tmp == 0)
res[i] = -1;
if (i%2 == 0)
i += i^(i & (i-1));
else
i++;
}
else
res[i] = ffs(tmp);
printf("Error in monomial_enumerate_lex\n");
}
return res;
return count;
}
int monomial_rank_convert(int n, int v, struct monomial_ctx_t *mono_1, struct monomial_ctx_t *mono_2)
int monomial_rank_lex(monomial_t m, struct monomial_ctx_t *mono)
{
monomial_t tmp = mono_1->monomials[n];
tmp = tmp >> v;
return monomial_rank(tmp, mono_2);
return binary_search(m, 0, mono->Nd[mono->n_var], mono);
}
struct polynomial_system **new_poly(struct polynomial_system *poly_sys)
int create_monomial_lex(int degree, int *vars, struct monomial_ctx_t *mono)
{
struct monomial_ctx_t *mono_A = monomial_setup(poly_sys->n - poly_sys->v, poly_sys->v, D);
struct monomial_ctx_t *mono_B = monomial_setup(poly_sys->n - poly_sys->v, poly_sys->v, D);
monomial_t term = 0;
for (int i = 0; i < degree; i++)
term += mkvar(vars[i]);
int mon = monomial_rank_lex(term, mono);
// assert(hamming_weight(mono->monomials[mon]) <= 2);
// assert(mon < mono->Nd[3]);
return mon;
}
int *tab = known_lin(poly_sys);
bool **poly_A = (bool **)malloc(sizeof(bool) * poly_sys->m * poly_sys->v);
bool **poly_B = (bool **)malloc(sizeof(bool) * poly_sys->m);
void test_lex(struct monomial_ctx_t *mono)
{
monomial_t tmp = -1;
for (int i = 0; i < mono->Nd[mono->n_var]; i++)
{
if (monomial_rank_lex(mono->monomials[i], mono) != i)
printf("Error in lex ranking\n");
if (hamming_weight(mono->monomials[i]) > mono->D)
printf("Error in hamming_weight lex\n");
if (mono->monomials[i] <= tmp)
printf("Error not in lexicographical order !\n");
tmp = mono->monomials[i];
}
// printf("test_lex :\nEZ\n");
}
static void monomial_setup_list_lex(struct monomial_ctx_t *mono)
{
int D = mono->D;
for (int i = 0; i < poly_sys->m; i++)
{
for (int k = 0; k < poly_sys->v; k++)
{
if (k == 0)
poly_B[i] = (bool *)malloc(sizeof(bool) * poly_sys->mono->Nd[D+1]);
poly_A[i*poly_sys->v + k] = (bool *)malloc(sizeof(bool) * poly_sys->mono->Nd[D+1]);
for (int c = 0; c < poly_sys->mono->Nd[D+1]; c++)
{
if (k == 0)
poly_B[i][c] = 0;
poly_A[i*poly_sys->v + k][c] = 0;
assert(D >= 2);
double start = wtime();
log(0, "Monomials setup\n");
log(1, " - Counting monomials\n");
mono->Nd[0] = 1;
for (int n = 0; n < mono->n_var; n++) {
int binom_n_d = monomial_enumerate_lex(D, NULL, n+1); /* not the smartest way to do the job... */
mono->Nd[n + 1] = mono->Nd[n] + binom_n_d;
// log(2, " - %d monomials of degree %d\n", binom_n_d, d);
}
for (int j = 0; j < poly_sys->mono->Nd[D+1]; j++)
{
if (tab[poly_sys->poly[i][j]] == -1)
log(2, " - %d monomials of degree <= %d\n", mono->Nd[mono->n_var], D);
log(1, " - Building monomial list\n");
mono->monomials = alloc_array(mono->Nd[mono->n_var], sizeof(*(mono->monomials)), "monomial list");