Revision 132f63c1a343af4280df47cf9cc0453060d506b7 authored by Jérémie Detrey on 10 March 2014, 10:29:14 UTC, committed by Jérémie Detrey on 10 March 2014, 10:31:12 UTC
1 parent b2cdbc1
check_rels.c
/* check_rels: check the factorization of the norm on both alg and rat side.
Can also, in option, check the primality of ideal and correct uncomplete
relation or relation with non prime ideal */
#include "cado.h"
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <inttypes.h>
#include <ctype.h>
#include <time.h>
#include <limits.h> /* for CHAR_BIT */
#include <unistd.h>
#include <assert.h>
#include <fcntl.h> /* for _O_BINARY */
#include "portability.h"
#include "macros.h"
#include "utils_with_io.h"
#include "relation.h"
#include "mpz_poly.h"
static index_t nrels_read, nrels_ok, nrels_err, nrels_completed, nrels_noprime;
static index_t nrels_toolarge;
cado_poly cpoly;
unsigned long lpb[2] = {0, 0};
int verbose = 0;
int abhexa = 0;
int check_primality = 0; /* By default no primality check */
int complete_rels = 0; /* By default, we just check the rels */
/* used for counting time in different processes */
timingstats_dict_t stats;
void
rel_add_prime (earlyparsed_relation_ptr rel, unsigned int side, p_r_values_t p,
exponent_t e)
{
for(weight_t i = 0; i < rel->nb ; i++)
{
if (rel->primes[i].p == p && rel->primes[i].h == side)
{
rel->primes[i].e += e;
return;
}
}
if (rel->nb == rel->nb_alloc)
realloc_buffer_primes_c (rel);
rel->primes[rel->nb] = (prime_t) {.h = side, .p = p, .e = e};
rel->nb++;
}
static inline void
print_error_line(uint64_t num, int64_t a, uint64_t b, int err, int already)
{
if (!already)
{
char *str = (err) ? "Error" : "Warning";
fprintf (stderr, "%s with relation %" PRIu64 " (a,b) = "
"(%" PRId64 ",%" PRIu64 "):\n", str, num, a, b);
}
}
static inline void
factor_nonprime_ideal(earlyparsed_relation_ptr rel, weight_t i)
{
exponent_t e = rel->primes[i].e;
p_r_values_t p = rel->primes[i].p;
unsigned int side = rel->primes[i].h, pr = 2;
modulusul_t m;
do
{
exponent_t e_pr_in_p = 0;
while (p % pr == 0)
{
p = p / pr;
e_pr_in_p++;
}
if (p == 1)
{
rel->primes[i] = (prime_t) {.h= side, .p= pr, .e= e * e_pr_in_p};
break;
}
else if (e_pr_in_p != 0)
rel_add_prime (rel, side, pr, e * e_pr_in_p);
pr = getprime(pr);
modul_initmod_ul (m, p);
} while (!modul_isprime(m));
getprime(0);
modul_clearmod (m);
if (p != 1) //means remaining p is prime
rel->primes[i] = (prime_t) {.h = side, .p = p, .e = e};
}
/* return 0 if everything is ok (factorization, primality, and complete)
* return -1 if a ideal does not divide the norm (primality and completeness
* are not checked)
* return r = r2*2^2 + r1*2 + r0 > 0 if all ideals divides the norm but:
* r0 = 1 means that an ideal was not prime
* r1 = 1 means that the factorization of a norm was not complete
* r2 = 1 means that an ideal was larger that a lpb
* If check_primality = 0, primality checks are not done and all ideals are
* supposed prime (so r0 is always 0).
* If complete_rels != 0, at the end, if r > 0, the relation is correct
* (it is completed, all ideals are below the lpbs and all non-prime are
* factored if check_primality != 0)
*/
int
process_one_relation (earlyparsed_relation_ptr rel)
{
mpz_t norm[2];
mpz_init (norm[0]);
mpz_init (norm[1]);
unsigned int fac_error = 0, need_completion = 0, nonprime = 0, toolarge = 0;
/* If we complete relations, do not print error msg but warning. */
int err = (complete_rels) ? 0 : 1;
/* compute the norm on alg and rat sides */
for(unsigned int side = 0 ; side < 2 ; side++)
{
mpz_poly_ptr ps = cpoly->pols[side];
mpz_poly_homogeneous_eval_siui (norm[side], ps, rel->a, rel->b);
}
/* check for correctness of the factorization of the norms */
for(weight_t i = 0; i < rel->nb ; i++)
{
unsigned int side = rel->primes[i].h;
p_r_values_t p = rel->primes[i].p;
exponent_t e = rel->primes[i].e;
for (int j = 0; j < e; ++j)
{
//if (mpz_fdiv_q_ui (norm[side], norm[side], p) != 0)
if (!mpz_divisible_ui_p (norm[side], p))
{
if (verbose != 0)
{
print_error_line (rel->num, rel->a, rel->b, 1, fac_error);
fprintf (stderr, " given factor %" PRpr " with exponent %u does "
"not divide the norm on side %u\n", p, e, side);
}
fac_error = 1;
}
else
mpz_divexact_ui (norm[side], norm[side], p);
}
}
/* With an error in the factorization of the norm, no need to continue */
if (fac_error)
{
mpz_clear(norm[0]);
mpz_clear(norm[1]);
return -1;
}
/* check primality of all ideals appearing in the relations */
if (check_primality != 0)
{
weight_t len = rel->nb; // no need to check the prime that can be added
for(weight_t i = 0; i < len ; i++)
{
modulusul_t p;
modul_initmod_ul (p, rel->primes[i].p);
if (!modul_isprime(p))
{
if (verbose != 0)
{
unsigned int side = (unsigned int) rel->primes[i].h;
print_error_line (rel->num, rel->a, rel->b, err, nonprime);
fprintf (stderr, " given factor %" PRpr " is not prime on side "
"%u\n", rel->primes[i].p, side);
}
nonprime = 1;
if (complete_rels) // if complete_rels = 1, we factor it
factor_nonprime_ideal(rel, i);
}
modul_clearmod (p);
}
}
/* complete relation if it is asked and necessary */
if (complete_rels)
{
for (unsigned long p = 2; mpz_cmp_ui (norm[0], 1) != 0 ||
mpz_cmp_ui (norm[1], 1) != 0 ; p = getprime(p))
{
for(unsigned int side = 0 ; side < 2 ; side++)
{
exponent_t e = 0;
while (mpz_divisible_ui_p (norm[side], p))
{
e++;
mpz_divexact_ui (norm[side], norm[side], p);
}
if (e != 0)
{
rel_add_prime (rel, side, p, e);
if (verbose != 0)
{
print_error_line (rel->num, rel->a, rel->b, err,
need_completion | nonprime);
gmp_fprintf (stderr, " factorization of the norm on side %u is "
"not complete, need to add %" PRpr "^%u\n",
side, p, e);
}
need_completion = 1;
}
}
}
getprime(0);
}
else
{
for(unsigned int side = 0 ; side < 2 ; side++)
{
if (mpz_cmp_ui (norm[side], 1) != 0)
{
if (verbose != 0)
{
print_error_line (rel->num, rel->a, rel->b, err,
need_completion | nonprime);
gmp_fprintf (stderr, " factorization of the norm on side %u is not "
"complete (%Zu is missing)\n", side, norm[side]);
}
need_completion = 1;
}
}
}
/* check that ideals appearing in the relations are below the lpb */
if (lpb[0] != 0 || lpb[1] != 0)
{
for(weight_t i = 0; i < rel->nb ; i++)
{
p_r_values_t p = rel->primes[i].p;
unsigned int side = rel->primes[i].h;
if (lpb[side] != 0 && p > lpb[side])
{
if (verbose != 0)
{
print_error_line (rel->num, rel->a, rel->b, err,
nonprime | need_completion | toolarge);
fprintf (stderr, " given factor %" PRpr " is greater than lpb "
"(=%lu) on side %u\n", p, lpb[side], side);
}
toolarge = 1;
}
}
}
mpz_clear(norm[0]);
mpz_clear(norm[1]);
return nonprime + 2*need_completion + 4*toolarge;
}
static inline void
print_relation (FILE *outfile, earlyparsed_relation_ptr rel)
{
char buf[1 << 12], *p, *op;
size_t t;
unsigned int i, j;
if (!abhexa)
{
p = d64toa10(buf, rel->a);
*p++ = ',';
p = u64toa10(p, rel->b);
*p++ = ':';
}
else
{
p = d64toa16(buf, rel->a);
*p++ = ',';
p = u64toa16(p, rel->b);
*p++ = ':';
}
for(unsigned int side = 0 ; side < 2 ; side++)
{
*(--p) = ':';
p++;
for (i = 0; i < rel->nb; i++)
{
if (rel->primes[i].h == side && rel->primes[i].e > 0)
{
op = p;
p = u64toa16(p, (uint64_t) rel->primes[i].p);
*p++ = ',';
t = p - op;
for (j = (unsigned int) ((rel->primes[i].e) - 1); j--; p += t)
memcpy(p, op, t);
}
}
}
*(--p) = '\n';
p[1] = 0;
fputs(buf, outfile);
}
/* Callback function called by filter_rels */
static void *
thread_callback (void * context_data, earlyparsed_relation_ptr rel)
{
FILE *outfile = (FILE *) context_data;
nrels_read++;
int ret = process_one_relation (rel);
if (complete_rels)
{
/* if complete_rels != 0, the only two cases where a rel is wrong are
when ret=-1 or ret >= 4 */
if (ret == -1)
nrels_err++;
else if (ret >= 4)
{
nrels_err++;
nrels_toolarge++;
}
else
{
nrels_ok++;
if (ret % 2) // rel contained a non prime ideal
nrels_noprime++;
if (ret >= 2) // factorization of a norm of rel was no complete
nrels_completed++;
print_relation (outfile, rel);
}
}
else
{
/* if complete_rels = 0, the only case where a rel is ok is ret = 0 */
if (ret == 0)
nrels_ok++;
else
{
nrels_err++;
if (ret > 0)
{
if (ret % 2) // rel contain a non prime ideal
nrels_noprime++;
ret /= 2;
if (ret % 2) // factorization of a norm of rel is no complete
nrels_completed++;
ret /= 2;
if (ret % 2) // an ideal was larger than a lpb
nrels_toolarge++;
}
}
}
return NULL;
}
static void
usage(const char *argv0)
{
fprintf (stderr, "Usage: %s [options] ", argv0);
fprintf (stderr, "[ -filelist <fl> [-basepath <dir>] | file1 ... filen ]\n");
fprintf (stderr, "Mandatory command line options:\n");
fprintf (stderr, " -poly <file> - polynomials file\n");
fprintf (stderr, "\nOther command line options:\n");
fprintf (stderr, " -abhexa - read and write a and b as hexa "
"(instead of decimal)\n");
fprintf (stderr, " -check_primality - check primality of ideal "
"(by default no checking)\n");
fprintf (stderr, " -complete <file> - write rels in file. If possible "
"incorrect rels are corrected\n");
fprintf (stderr, " -lpb0 <l> - chech that ideals on side 0 are "
"below 2^l\n");
fprintf (stderr, " -lpb1 <l> - chech that ideals on side 1 are "
"below 2^l\n");
fprintf (stderr, " -v - more verbose output\n");
exit (1);
}
int
main (int argc, char * argv[])
{
char * argv0 = argv[0];
param_list pl;
param_list_init(pl);
argv++,argc--;
FILE *outfile = NULL;
param_list_configure_switch(pl, "abhexa", &abhexa);
param_list_configure_switch(pl, "v", &verbose);
param_list_configure_switch(pl, "check_primality", &check_primality);
#ifdef HAVE_MINGW
_fmode = _O_BINARY; /* Binary open for all files */
#endif
for( ; argc ; ) {
if (param_list_update_cmdline(pl, &argc, &argv)) { continue; }
/* Since we accept file names freeform, we decide to never abort
* on unrecognized options */
break;
}
/* Update parameter list at least once to register argc/argv pointers. */
param_list_update_cmdline (pl, &argc, &argv);
/* print command-line arguments */
param_list_print_command_line (stdout, pl);
fflush(stdout);
param_list_parse_ulong(pl, "lpb0", &lpb[0]);
param_list_parse_ulong(pl, "lpb1", &lpb[1]);
lpb[0] = (lpb[0] == 0) ? 0 : 1UL << lpb[0];
lpb[1] = (lpb[1] == 0) ? 0 : 1UL << lpb[1];
const char * polyfilename = param_list_lookup_string(pl, "poly");
const char *outfilename = param_list_lookup_string(pl, "complete");
const char * filelist = param_list_lookup_string(pl, "filelist");
const char * basepath = param_list_lookup_string(pl, "basepath");
const char * path_antebuffer = param_list_lookup_string(pl, "path_antebuffer");
set_antebuffer_path (argv0, path_antebuffer);
if (param_list_warn_unused(pl))
usage(argv0);
if (basepath && !filelist)
{
fprintf(stderr, "-basepath only valid with -filelist\n");
usage(argv0);
}
if (outfilename)
complete_rels = 1;
if ((filelist != NULL) + (argc != 0) != 1)
{
fprintf(stderr, "Provide either -filelist or freeform file names\n");
usage(argv0);
}
if (!polyfilename)
usage(argv0);
cado_poly_init (cpoly);
if (!cado_poly_read (cpoly, polyfilename))
{
fprintf (stderr, "Error reading polynomial file\n");
exit (EXIT_FAILURE);
}
char ** files = filelist ? filelist_from_file(basepath, filelist, 0) : argv;
nrels_read = nrels_err = nrels_ok = nrels_completed = nrels_noprime = 0;
nrels_toolarge = 0;
if (complete_rels)
{
outfile = fopen_maybe_compressed (outfilename, "w");
printf("# Correct relations will be written in %s\n", outfilename);
}
printf ("# Verbose output: %s\n", verbose ? "yes" : "no");
printf ("# Correct wrong relations if possible: %s\n",
complete_rels ? "yes" : "no");
printf ("# Check primality of ideal: %s\n", check_primality ? "yes" : "no");
timingstats_dict_init(stats);
filter_rels(files, (filter_rels_callback_t) &thread_callback, (void*)outfile,
EARLYPARSE_NEED_PRIMES |
(abhexa ? EARLYPARSE_NEED_AB_HEXA : EARLYPARSE_NEED_AB_DECIMAL),
NULL, stats);
printf("Number of read relations: %" PRid "\n", nrels_read);
if (complete_rels)
{
printf("Number of deleted relations: %" PRid "\n", nrels_err);
printf(" among which %" PRid " contained an ideal larger than a lpb\n",
nrels_toolarge);
printf("Number of keeped relations: %" PRid "\n", nrels_ok);
printf(" among which %" PRid " were completed\n"
" and %" PRid " contained at least one "
"non-primes ideal\n", nrels_completed, nrels_noprime);
}
else
{
printf("Number of correct relations: %" PRid "\n", nrels_ok);
printf("Number of wrong relations: %" PRid "\n", nrels_err);
printf(" among which %" PRid " were not complete\n"
" and %" PRid " contained an ideal larger than a lpb\n"
" and %" PRid " contained at least one "
"non-primes ideal\n", nrels_completed, nrels_toolarge,
nrels_noprime);
}
if (filelist)
filelist_clear(files);
if (complete_rels)
fclose_maybe_compressed (outfile, outfilename);
param_list_clear(pl);
cado_poly_clear (cpoly);
timingstats_dict_add_mythread(stats, "main");
timingstats_dict_disp(stats);
timingstats_dict_clear(stats);
return (nrels_err) ? EXIT_FAILURE : EXIT_SUCCESS ;
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...