Revision 0493cb82f256ff6d122b63d0b5226908965ec87f authored by Paul Zimmermann on 17 February 2014, 16:32:26 UTC, committed by Paul Zimmermann on 17 February 2014, 16:32:26 UTC
1 parent 9caacb8
explain.c
/* This takes a range of row numbers in the post-merge relation matrix,
* and checks the corresponding factorization into rational primes and
* algebraic prime ideals.
*/
/* Usage string:
*
$bindir/filter/explain -sosp $wdir/$name.sos-purged.bin -sosr $wdir/$name.sos-replay.bin -purged $wdir/$name.purged -index $wdir/$name.index -matrix $wdir/$name.small.bin -skip 32 -start 0 -end 10
* Do not forget the skip parameter -- otherwise the output is completely
* useless.
*
* The file name.sos-replay.bin is created by replay.
* The file sos-purged.bin is created by purge in text format, the
* following snippet converts it to binary (yeah, purge should do this).
#include <stdio.h>
main()
{
unsigned int i,p,r;
for( ; scanf("%u %x %x", &i,&p,&r) == 3 ; ) {
fwrite(&p, sizeof(unsigned int), 1, stdout);
fwrite(&r, sizeof(unsigned int), 1, stdout);
}
}
*
* -start defaults to 0
* -end defaults to nrows
*/
#include "cado.h"
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include "portability.h"
#include "utils.h"
/* copy my_malloc_int interface from purge {{{ */
#define BLOCK_SIZE 1000000 /* memory blocks are allocated of that # of int's */
/* mmi_list is a list of blocks, each one of BLOCK_SIZE ints */
int **mmi_list = NULL;
int mmi_current = -1; /* index of current block */
unsigned long mmi_used = BLOCK_SIZE; /* usage of current block */
/* converts an index to a pointer */
static int* mmi_translate(int v)
{
int block = v / BLOCK_SIZE;
return mmi_list[block] + (v % BLOCK_SIZE);
}
/* return the index of an array of n ints */
static int
mmi_get (unsigned long n)
{
int index;
if (mmi_used + n > BLOCK_SIZE) {
/* new block */
mmi_current ++;
mmi_list = (int**) realloc (mmi_list, (mmi_current + 1) * sizeof (int*));
mmi_list[mmi_current] = (int*) malloc (BLOCK_SIZE * sizeof (int));
mmi_used = 0;
}
index = mmi_current * BLOCK_SIZE + mmi_used;
mmi_used += n;
return index;
}
static void
mmi_clear (void)
{
while (mmi_current >= 0)
free (mmi_list[mmi_current--]);
free (mmi_list);
mmi_list = NULL;
mmi_used = BLOCK_SIZE;
}
/* }}} */
typedef int (*sortfunc_t) (const void*, const void*);
int uint32cmp(const uint32_t * a, const uint32_t * b)
{
return (*a>*b)-(*b>*a);
}
struct ideal {
uint32_t p;
uint32_t r;
};
struct ideal * ideal_decode;
int intcmp2(const int * a, const int * b)
{
if (ideal_decode[*a].r == -2u && ideal_decode[*b].r != -2u)
return -1;
if (ideal_decode[*a].r != -2u && ideal_decode[*b].r == -2u)
return 1;
int r = uint32cmp(&ideal_decode[*a].p, &ideal_decode[*b].p);
if (r) return r;
r = uint32cmp(&ideal_decode[*a].r, &ideal_decode[*b].r);
return r;
}
int
main (int argc, char **argv)
{
param_list pl;
param_list_init(pl);
argv++,argc--;
for( ; argc ; ) {
if (param_list_update_cmdline(pl, &argc, &argv)) { continue; }
fprintf (stderr, "Unknown option: %s\n", argv[0]);
abort();
}
const char * sosp = param_list_lookup_string(pl, "sosp");
const char * sosr = param_list_lookup_string(pl, "sosr");
const char * purged = param_list_lookup_string(pl, "purged");
const char * index = param_list_lookup_string(pl, "index");
const char * matrix = param_list_lookup_string(pl, "matrix");
const char * outname = param_list_lookup_string(pl, "out");
int skip=0;
param_list_parse_int(pl, "skip", &skip);
ASSERT_ALWAYS(sosp && sosr && purged && index && matrix);
/* just avoid warnings. */
param_list_lookup_string(pl, "start");
param_list_lookup_string(pl, "end");
if (param_list_warn_unused(pl))
exit(1);
FILE * outfile = stdout;
int pipe = 0;
if (outname) {
outfile = fopen_maybe_compressed2(outname, "w", &pipe,NULL);
}
/* Prepare the list of (a,b) pairs */
unsigned int n_ab;
unsigned int n_ideals;
size_t nb_stored_ints_purged = 0;
struct abpair {
int64_t a;
uint64_t b;
int mmi_index;
};
struct abpair * abs;
{
purgedfile_stream ps;
purgedfile_stream_init(ps);
purgedfile_stream_openfile(ps, purged);
n_ab = ps->nrows;
n_ideals = ps->ncols;
fprintf(stderr, "Allocating room for %u abs: %.2f MB\n",
n_ab,
(ps->nrows * sizeof(struct abpair)) /1048576.0);
abs = malloc(ps->nrows * sizeof(struct abpair));
ASSERT_ALWAYS(abs);
for(int i = 0 ; purgedfile_stream_get(ps, NULL) >= 0 ; i++) {
abs[i].a = ps->a;
abs[i].b = ps->b;
abs[i].mmi_index = mmi_get(ps->nc + 1);
int * ptr = mmi_translate(abs[i].mmi_index);
*ptr++ = ps->nc;
for(int k = 0 ; k < ps->nc ; k++) {
ptr[k] = ps->cols[k];
}
nb_stored_ints_purged += (ps->nc + 1);
}
fprintf(stderr, "Storing (a,b) pairs factorizations: %.2f MB\n",
nb_stored_ints_purged * sizeof(int) /1048576.0);
purgedfile_stream_closefile(ps);
purgedfile_stream_clear(ps);
}
/* preparing the correspondence between ideals appearing in .purged,
* and real ideals
*/
{
fprintf(stderr, "Allocating room for %u ideal descriptons: %.2f MB\n",
n_ideals,
(n_ideals * sizeof(struct ideal)) /1048576.0);
ideal_decode = malloc(n_ideals * sizeof(struct ideal));
ASSERT_ALWAYS(ideal_decode);
FILE * f = fopen(sosp, "r"); ASSERT_ALWAYS(f);
int r = fread(ideal_decode, sizeof(struct ideal), n_ideals, f);
ASSERT_ALWAYS(n_ideals == (unsigned int) r);
fclose(f);
}
/* preparing to decode ideals which appear in the matrix */
unsigned int n_ideals_inmat;
uint32_t * ideals_inmat;
{
struct stat sbuf[1];
if (stat(sosr, sbuf) < 0) { perror(sosr); exit(1); }
n_ideals_inmat = sbuf->st_size / sizeof(uint32_t);
fprintf(stderr, "Allocating room for %u ideal indirections: %.2f MB\n",
n_ideals_inmat,
(n_ideals_inmat * sizeof(uint32_t)) /1048576.0);
ideals_inmat = malloc(n_ideals_inmat * sizeof(uint32_t));
ASSERT_ALWAYS(ideals_inmat);
FILE * f = fopen(sosr, "r"); ASSERT_ALWAYS(f);
int r = fread(ideals_inmat, sizeof(uint32_t), n_ideals_inmat, f);
ASSERT_ALWAYS(n_ideals_inmat == (unsigned int) r);
fclose(f);
}
/* read relation-sets and matrix rows simultaneously */
{
FILE * f = fopen(index, "r"); ASSERT_ALWAYS(f);
unsigned int nrows;
unsigned int ncols;
int r = fscanf(f, "%u %u",&nrows, &ncols);
unsigned int i0 = 0;
unsigned int i1 = nrows;
param_list_parse_uint(pl, "start", &i0);
param_list_parse_uint(pl, "end", &i1);
ASSERT_ALWAYS(r == 2);
ASSERT_ALWAYS(ncols == n_ideals_inmat);
FILE * g = fopen(matrix, "r"); ASSERT_ALWAYS(f);
fprintf(stderr, "Dumping relset info for relsets [%d..%d[ (complete range is [%d..%d[\n",i0,i1,0,nrows);
for(unsigned int i = 0 ; i < i0 ; i++) {
int nc;
r = fscanf(f, "%u", &nc);
ASSERT_ALWAYS(r == 1);
for(int j = 0 ; j < nc ; j++) {
unsigned int k;
r = fscanf(f, "%x", &k);
ASSERT_ALWAYS(r == 1);
}
uint32_t z;
int r = fread(&z, sizeof(uint32_t), 1, g);
ASSERT_ALWAYS(r == 1);
for(int j = 0 ; j < nc ; j++) {
r = fread(&z, sizeof(uint32_t), 1, g);
ASSERT_ALWAYS(r == 1);
}
}
int * fac = NULL;
int fac_alloc = 0;
int fac_size = 0;
for(unsigned int i = i0 ; i < i1 ; i++) {
/* Use the index file to decode the relation-set */
int nc;
fac_size = 0;
r = fscanf(f, "%u", &nc);
ASSERT_ALWAYS(r == 1);
fprintf(outfile, "relset %u\n", i);
for(int j = 0 ; j < nc ; j++) {
unsigned int k;
r = fscanf(f, "%x", &k);
ASSERT_ALWAYS(r == 1);
ASSERT_ALWAYS(k < n_ab);
fprintf(outfile, "%" PRId64 " %" PRIu64 "\n", abs[k].a, abs[k].b);
int * ff = mmi_translate(abs[k].mmi_index);
if (fac_size + *ff > fac_alloc) {
fac_alloc += *ff + 64 + fac_alloc / 2;
fac = realloc(fac, fac_alloc * sizeof(int));
}
memcpy(fac + fac_size, ff + 1, *ff * sizeof(int));
fac_size += *ff;
}
fprintf(outfile, "expected\n");
qsort(fac, fac_size, sizeof(int), (sortfunc_t)&intcmp2);
for(int s = 0; s < fac_size ;) {
int idl = fac[s];
int v = 0;
for( ; s + v < fac_size && fac[s+v] == idl ; v++) ;
ASSERT_ALWAYS((unsigned int) idl < n_ideals);
fprintf(outfile, "%d %d %d\n", v, ideal_decode[idl].p, ideal_decode[idl].r);
s+=v;
}
fprintf(outfile, "recorded\n");
uint32_t z;
int r = fread(&z, sizeof(uint32_t), 1, g);
ASSERT_ALWAYS(r == 1);
nc = z;
fac_size = 0;
if (fac_size + nc > fac_alloc) {
fac_alloc += nc + 64 + fac_alloc / 2;
fac = realloc(fac, fac_alloc * sizeof(int));
}
for(int j = 0 ; j < nc ; j++) {
r = fread(&z, sizeof(uint32_t), 1, g);
z += skip;
ASSERT_ALWAYS(r == 1);
ASSERT_ALWAYS(z < n_ideals_inmat);
int idl = ideals_inmat[z];
fac[fac_size++]=idl;
}
qsort(fac, fac_size, sizeof(int), (sortfunc_t)&intcmp2);
for(int s = 0; s < fac_size ;) {
int idl = fac[s];
int v = 0;
for( ; s + v < fac_size && fac[s+v] == idl ; v++) ;
ASSERT_ALWAYS((unsigned int) idl < n_ideals);
fprintf(outfile, "%d %d %d\n", v, ideal_decode[idl].p, ideal_decode[idl].r);
s+=v;
}
fprintf(outfile, "\n");
}
free(fac);
fclose(f);
fclose(g);
}
if (outname) {
if (pipe)
pclose(outfile);
else
fclose(outfile);
}
mmi_clear();
param_list_clear(pl);
return 0;
}
Computing file changes ...