Commit b552e184 authored by Charles Bouillaguet's avatar Charles Bouillaguet
Browse files

for publication

parent f562b4fe
CC = gcc
CFLAGS = -std=c11 -g -Wall -Wextra -O3 -Wno-unused-parameter -Werror -march=native -mavx2
CFLAGS = -std=c11 -g -Wall -Wextra -O3 -Wno-unused-parameter -Wno-maybe-uninitialized -Werror -march=native -mavx2
# OpenMP
CFLAGS += -fopenmp
LDFLAGS += -fopenmp
LDLIBS += -lm
all: parser_demo parser_demo_alt moebius_demo monica monica_packed monica_oblivious
all: parser_demo parser_demo_alt moebius_demo monica monica_vector
monica: parser.o monica.o
monica_packed: parser.o monica_packed.o
monica_oblivious: parser.o monica_oblivious.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
moebius_demo.o: parser.h moebius.h
parser_demo.o: parser.h
parser.o: parser.h
moebius.o: moebius.h
monica.o: parser.h ffs.h
monica_packed.o: parser.h ffs.h
monica_oblivious.o: parser.h ffs.h
monica_vector.o: parser.h ffs.h
clean:
rm -f *.o
rm -f parser_demo moebius_demo parser_demo_alt monica monica_packed monica_oblivious
rm -f parser_demo moebius_demo parser_demo_alt monica monica_vector
# Programs for solving Boolean polynomial systems
This code is in the public domain. Several algorithms are implemented.
Exhaustive search is available elsewhere, in the [libfes-lite](https://gitlab.lip6.fr/almasty/libfes-lite) library.
# Monica
This program implements a striped-down version of the [Crossbred algorithm of Joux-Vitse](https://hal.archives-ouvertes.fr/PRISM/hal-01981516v1). It handles at most 64 variables. At this point, it is not multi-threaded.
Internally, it uses "libfes-style" gray-code enumeration to evaluate polynomials quickly.
There are two versions : `monica.c`, which is as simple as possible, and `monica_vector.c`,
which tries to use the fact that several bits fits into a single word to speed things up.
The most important parameter is `--inner-hybridation`, which specifies how
many variables are not guessed. It may make sense to try one less than the default value.
Example:
```
./monica_vector --inner-hybridation 7 < examples/random_40_quad.in
```
# Moebius
This small piece of code solves systems of arbitrary-degree boolean
polynomials, by "exhaustive search", using the moebius transform algorithm.
This code is in the public domain.
The example program requires a C11-capable compiler (because it uses `aligned_alloc`).
It is multi-threaded, using OpenMP. By default, it will run one thread by "processor",
which means per logical core.
......@@ -19,6 +42,7 @@ Example:
# ./moebius_demo < examples/large.in
```
# Input format
This code parses polynomial systems modulo 2 in the following format:
......@@ -34,17 +58,8 @@ This code parses polynomial systems modulo 2 in the following format:
(It follows that an empty line denote the zero polynomial).
# Parser
The parser is somewhat generic and could be useful in other contexts.
See `parser_demo.c` and `parser_demo_alt.c` for examples of use.
# Monica
This program implements a striped-down the algorithm of Monica Trimoska for quadratic
polynomials. It handles <= 64 variables. At this point, it is not multi-threaded.
Internally, it uses "libfes-style" gray-code enumeration to evaluate polynomials quickly.
There are two versions : `monica.c`, which is a simple as possible, and `monica_packed.c`,
which tries to use the fact that several bits fits into a single word to speed things up.
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -459,11 +459,8 @@ void oblivious_solve(const struct vector_enumeration_state_t *ves, int l, int v,
A[i][j] = ves[i].Ac[j];
}
// TODO : merge elimination with next pivot search (keep stuff in registers)
/* solve Ax == b */
for (int j = 0; j < v; j++) {
/* find a pivot on column j */
/* find a pivot on column 0 */
int j = 0;
for (int i = j + 1; i < l; i++) {
vector swap = A[i][j] & (~A[j][j]);
// conditionnaly swap rows i and j
......@@ -477,16 +474,45 @@ void oblivious_solve(const struct vector_enumeration_state_t *ves, int l, int v,
b[j] ^= xor;
}
/* solve Ax == b */
for (int j = 0; j < v - 1; j++) {
// eliminate everything below row j
for (int i = j + 1; i < l; i++) {
vector add = A[i][j];
vector add = A[j+1][j];
for (int k = j + 1; k < v; k++)
A[j+1][k] ^= add & A[j][k];
b[j+1] ^= add & b[j];
/* interleaving of elimination on column j + pivot search/row permutation on column j+1 */
for (int i = j + 2; i < l; i++) {
vector add = A[i][j];
A[i][j+1] ^= add & A[j][j+1];
vector swap = A[i][j+1] & (~A[j+1][j+1]); // swap rows i <----> j+1 ?
vector xor = swap & (A[i][j+1] ^ A[j+1][j+1]);
A[i][j+1] ^= xor;
A[j+1][j+1] ^= xor;
for (int k = j + 2; k < v; k++) {
A[i][k] ^= add & A[j][k];
vector xor = swap & (A[i][k] ^ A[j+1][k]);
A[i][k] ^= xor;
A[j+1][k] ^= xor;
}
b[i] ^= add & b[j];
vector xorb = swap & (b[i] ^ b[j+1]);
b[i] ^= xorb;
b[j+1] ^= xorb;
}
}
// eliminate everything below row j + consistency check
b[v] ^= A[v][v-1] & b[v-1];
vector inconsistent = b[v];
for (int i = v; i < l; i++) {
b[i] ^= A[i][v-1] & b[v-1];
inconsistent |= b[i];
}
*consistent = ~inconsistent;
// backward substitution
// backward substitution (this is only needed to actually get the solutions...)
for (int i = v-2; i >= 0; --i)
for (int j = i + 1; j < v; j++)
b[i] ^= A[i][j] & b[j];
......@@ -496,30 +522,25 @@ void oblivious_solve(const struct vector_enumeration_state_t *ves, int l, int v,
for (int j = 1; j < v; j++)
full_rank &= A[j][j];
*rank_defect = ~full_rank;
vector inconsistent = b[v];
for (int i = v + 1; i < l; i++)
inconsistent |= b[i];
*consistent = ~inconsistent;
}
/*************************** processing candidate solutions***************************/
/* statistics */
u64 n_special = 0; /* consistent AND rank-defective systems encountered */
u64 n_expensive = 0; /* "expensive" iterations where there are candidates */
u64 n_consistent = 0; /* consistent systems encountered */
u64 n_defective = 0; /* rank-defective systems encountered */
u64 n_expensive = 0; /* "expensive" iterations where there are candidates */
/* evaluate x on the m initial polynomials, with early abort */
bool candidate_solution(u64 x) {
static bool candidate_solution(u64 x) {
for (int i = 0; i < m; i++)
if (eval_poly(&poly[i], n, x))
return 0;
return 1;
}
bool not_implemented_warning;
/* deal with this value of y, which triggers a consistent, rank-defective linear system */
void process_special_case(const struct symbolic_poly_t *sp, int u, int v, u64 y)
......@@ -527,9 +548,12 @@ void process_special_case(const struct symbolic_poly_t *sp, int u, int v, u64 y)
// assemble system
// solve generic
// test solutions
if (!not_implemented_warning)
warnx("potential solutions discarded because the complex case is not yet implemented");
not_implemented_warning = 1;
}
void process_candidates(int v, u64 y, vector consistent, vector defect, const vector *solution)
static void process_candidates(int v, u64 y, vector consistent, vector defect, const vector *solution)
{
if (vector_is_zero(consistent | defect))
return;
......@@ -652,7 +676,8 @@ int main(int argc, char **argv)
struct ffs_t ffs;
ffs_reset(&ffs, u);
ffs_step(&ffs);
for (u64 y = 0; y < (1ull << (u - VLEN)); y++) {
u64 n_iterations = 1ull << (u - VLEN);
for (u64 y = 0; y < n_iterations; y++) {
vector consistent;
vector defect;
vector solution[64];
......@@ -666,6 +691,13 @@ int main(int argc, char **argv)
for (int i = 0; i < l; i++)
update(&ves[i], &ffs, v);
ffs_step(&ffs);
if ((y & 0xffff) == 0) {
double rate = y / (wtime() - start);
double ETA = (n_iterations - y) / rate;
fprintf(stderr, "\ry = %lx / %lx ; ETA=%.0fs ; %.0f iteration/s", y, n_iterations, ETA, rate);
fflush(stderr);
}
}
double stop = wtime();
......
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