Revision e4f0108028701f575c3bf69e567316e3150f7dd1 authored by Paul Zimmermann on 13 October 2010, 08:19:28 UTC, committed by Paul Zimmermann on 13 October 2010, 08:19:28 UTC
[bwc] added "krylov" or "mksol" in ETA log lines git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/cado-nfs/trunk@352 3eaf19af-ecc0-40f6-b43f-672aa0c71c71
1 parent ebbb4f9
freerel.c
/*
* Program: free relations
* Author : F. Morain
* Purpose: creating free relations in a suitable format
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gmp.h>
#include "cado.h"
#include "mod_ul.c"
#include "utils.h"
#include "sieve/fb.h"
#include "hashpair.h"
// The best reference for free relations is Huizing (Exp. Math. 1996,
// Section 4). Handling them correctly in the factorization is also detailed.
// If f(X)=cK*X^dK+... and g(X)=cL*X^dL+... are the two polynomials used,
// then p is a free prime iff f(X) and g(X) split completely mod p into
// dK (resp. dL) *distinct* linear factors. In other words, we have to check
// that p does not divide the discriminants of f and g, and p doesn't divide
// cK*cL.
//
// TODO: make this definition go into the actual code, which is not the
// case currently.
// When p1 == p2, sort r's in increasing order
int
compare_ul2(const void *v1, const void *v2)
{
unsigned long w1 = *((unsigned long*) v1);
unsigned long w2 = *((unsigned long*) v2);
if(w1 > w2)
return 1;
if(w1 < w2)
return -1;
w1 = *(1 + (unsigned long*) v1);
w2 = *(1 + (unsigned long*) v2);
return (w1 >= w2 ? 1 : -1);
}
// find free relations by inspection.
int
findFreeRelations(hashtable_t *H, cado_poly pol, int nprimes)
{
unsigned long *tmp = (unsigned long *)malloc((2+(nprimes<<1)) * sizeof(unsigned long));
unsigned int i, j, k, ntmp = 0;
int pdeg, nfree;
for(i = 0; i < H->hashmod; i++)
if(H->hashcount[i] > 0){
tmp[ntmp++] = (unsigned long) GET_HASH_P(H,i);
tmp[ntmp++] = GET_HASH_R(H,i);
}
qsort(tmp, (ntmp>>1), 2 * sizeof(unsigned long), compare_ul2);
// add a sentinel
tmp[ntmp] = 0;
tmp[ntmp+1] = 0;
// now, inspect
for(i = 0; ; i += 2)
if(((long)tmp[i+1]) >= 0)
break;
j = i+2;
pdeg = 1;
nfree = 0;
while(1){
if(tmp[j] == tmp[i]){
// same prime
if(((long)tmp[j+1]) >= 0)
// another algebraic factor
pdeg++;
}
else{
// new prime
if(pdeg == pol[0].degree){
// a free relation
printf("%lu,0:%lx:", tmp[i], tmp[i]);
for(k = i; k < i+(pdeg<<1); k += 2){
printf("%lx", tmp[k+1]);
if(k < i+(pdeg<<1)-2)
printf(",");
}
printf("\n");
nfree++;
}
i = j;
if(tmp[j] == 0)
break;
pdeg = 1;
}
j += 2;
}
free(tmp);
return nfree;
}
/* On largeFreeRelations:
*
commit 7cb6c0e177ed7ba86c1768942bfb42f2e575522e
Author: morain <morain@de4796e2-0b1e-0410-91b4-86bc23549317>
Date: Tue Feb 12 12:41:46 2008 +0000
* Cleaned purge.c
* added a new option to freerel to look for (small and large) free relations
by inspection of all the relations found. Quite experimental at the time
being. Do not use it.
git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/cado/trunk@801 de4796e2-0b1e-0
*/
void
largeFreeRelations (cado_poly pol, char **fic, int verbose)
{
hashtable_t H;
int Hsizea, nprimes_alg = 0, nfree = 0;
int need64 = (pol->lpbr > 32) || (pol->lpba > 32);
ASSERT(fic != NULL);
/* The number of algebraic large primes is about 1/2*L/log(L)
where L is the algebraic large prime bound, i.e., L=2^lpba.
However since we store separately primes p and the corresponding root r
of f mod p, the number of (p,r) pairs is about L/log(L). */
Hsizea = (1 << pol[0].lpba) / ((int)((double) pol[0].lpba * log(2.0)));
hashInit (&H, Hsizea, verbose, need64);
if (verbose)
fprintf (stderr, "Scanning relations\n");
relation_stream rs;
relation_stream_init(rs);
// rs->sort_primes = 1;
// rs->reduce_mod2 = 1;
for( ; *fic ; fic++) {
relation_stream_openfile(rs, *fic);
if (verbose)
fprintf (stderr, "Adding file %s\n", *fic);
for( ; relation_stream_get(rs, NULL) >= 0 ; ) {
if (rs->rel.b == 0) {
fprintf(stderr, "Ignoring already found free relation...\n");
continue;
}
relation_compress_rat_primes(&rs->rel);
relation_compress_alg_primes(&rs->rel);
reduce_exponents_mod2(&rs->rel);
computeroots(&rs->rel);
for(int j = 0; j < rs->rel.nb_ap; j++){
int h = hashInsert(&H, rs->rel.ap[j].p, rs->rel.ap[j].r);
nprimes_alg += (H.hashcount[h] == 1); // new prime
}
}
relation_stream_closefile(rs);
if (verbose)
hashCheck (&H);
}
relation_stream_clear(rs);
if (verbose)
fprintf (stderr, "nprimes_alg = %d\n", nprimes_alg);
nfree = findFreeRelations(&H, pol, nprimes_alg);
if (verbose)
fprintf (stderr, "Found %d usable free relations\n", nfree);
hashFree(&H);
}
int
countFreeRelations(int *deg, char *roots)
{
FILE *ifile = fopen(roots, "r");
char str[1024], str0[128], str1[128], str2[128], *t;
int nfree = 0, nroots;
*deg = -1;
// look for the degree
while(1){
if(!fgets(str, 1024, ifile))
break;
sscanf(str, "%s %s %s", str0, str1, str2);
if((str0[0] == '#') && !strcmp(str1, "DEGREE:")){
*deg = atoi(str2);
break;
}
}
if(*deg == -1){
fprintf(stderr, "No degree found, sorry\n");
exit(1);
}
fprintf(stderr, "DEGREE = %d\n", *deg);
while(1){
if(!fgets(str, 1024, ifile))
break;
// format is "11: 2,6,10"
sscanf(str, "%s %s", str0, str1);
nroots = 1;
for(t = str1; *t != '\0'; t++)
if(*t == ',')
nroots++;
if(nroots == *deg)
nfree++;
}
fclose(ifile);
return nfree;
}
void
addFreeRelations(char *roots, int deg)
{
factorbase_degn_t *fb = fb_read(roots, 1.0, 0), *fbptr;
fbprime_t p;
int i;
for(fbptr = fb; fbptr->p != FB_END; fbptr = fb_next(fbptr)){
p = fbptr->p;
if(fbptr->nr_roots == deg){
// free relation, add it!
// (p, 0) -> p-0*m
printf("%d,0:%lx:", p, (long)p);
for(i = 0; i < fbptr->nr_roots; i++){
printf("%x", fbptr->roots[i]);
if(i < fbptr->nr_roots-1)
printf(",");
}
printf("\n");
}
}
}
static void
smallFreeRelations(char *fbfilename)
{
int deg;
int nfree = countFreeRelations (°, fbfilename);
fprintf (stderr, "# Free relations: %d\n", nfree);
fprintf (stderr, "Handling free relations...\n");
addFreeRelations (fbfilename, deg);
}
static void
usage (char *argv0)
{
fprintf (stderr, "Usage: %s -fb xxx.roots\n", argv0);
fprintf (stderr, "or %s [-v] -poly xxx.poly xxx.rels1 xxx.rels2 ... xxx.relsk\n", argv0);
exit (1);
}
int
main(int argc, char *argv[])
{
char *fbfilename = NULL, *polyfilename = NULL, **fic;
char *argv0 = argv[0];
cado_poly cpoly;
int nfic = 0;
int verbose = 0;
while (argc > 1 && argv[1][0] == '-')
{
if (argc > 2 && strcmp (argv[1], "-fb") == 0)
{
fbfilename = argv[2];
argc -= 2;
argv += 2;
}
else if (argc > 1 && strcmp (argv[1], "-v") == 0)
{
verbose ++;
argc --;
argv ++;
}
else if (argc > 2 && strcmp (argv[1], "-poly") == 0)
{
polyfilename = argv[2];
argc -= 2;
argv += 2;
}
else
usage (argv0);
}
fic = argv+1;
nfic = argc-1;
if (polyfilename == NULL || ((fbfilename == NULL) && (nfic == 0)))
usage (argv0);
cado_poly_init(cpoly);
if (!cado_poly_read (cpoly, polyfilename))
{
fprintf (stderr, "Error reading polynomial file\n");
exit (EXIT_FAILURE);
}
/* check that n divides Res(f,g) [might be useful to factor n...] */
check_polynomials (cpoly);
#if 0
if (mpz_cmp_ui (cpoly->g[1], 1) != 0)
{
fprintf (stderr, "Error, non-monic linear polynomial not yet treated");
fprintf (stderr, " (more theory needed)\n");
exit (1);
}
#endif
if (nfic == 0)
smallFreeRelations(fbfilename);
else
largeFreeRelations(cpoly, fic, verbose);
cado_poly_clear (cpoly);
return 0;
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...