mersenneforum.org - Msieve with GNFS support

[QUOTE=R.D. Silverman;132988]I'd like to try out the msieve solver, but my code emits and works
with the CWI format. Is there as CWI --> msieve relation converter???[/QUOTE]

It's a simple translation. The only annoying thing is that the GGNFS format (which is the format msieve uses) has a print bound of about 1000 if I remember correctly, so you need to refactor the relations. Here's the code I use. It was originally written by Jason, but has modifications to speed it up a bit. It calls routines from the GGNFS package, but it's simple enough to change to use your own routines.

Greg

[CODE]/**************************************************************/
/* cwi2gg.c */
/* Copyright 2004, Jason Papadopoulos. */
/**************************************************************/
/* This file is part of GGNFS.
*
* GGNFS is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* GGNFS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GGNFS; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <signal.h>
#include <time.h>
#include <gmp.h>
#include <ctype.h>
#include "ggnfs.h"

#define DEFAULT_FB_NAME "default.fb"

#define USAGE \
"[OPTIONS]\n"\
"--help : show this help and exit\n"\
"-a1 : algebraic factors come first in input file (default)\n"\
"-r1 : rational factors come first in input file\n"\
"-fb <filename> : use factor base in file <filename>\n"\
"-i <filename> : convert CWI siever relations in <filename>\n"\
"-j <filename> : use job file <filename>\n"

/*******************************************************/
int parseJobFile(nfs_sieve_job_t *job, char *fName)
/*******************************************************/
{ FILE *fp;
char token[256], value[256], thisLine[1024];
int t;

if (!(fp = fopen(fName, "r"))) {
printf("Error opening %s for read!\n", fName);
return -1;
}

readPoly(fp, &job->FB); /* read n, m, poly, etc. */
job->nDigits = mpz_sizeinbase(job->FB.n, 10);
rewind(fp); /* Back up to the beginning of the file. */

/* Initialize fields to empty, or defaults which are independent of
the number being factored. */
job->type = SIEVE_CLASSICAL;
job->a0 = job->a1 = job->b0 = job->b1 = 0;
job->qIndex0 = job->qIndex1 = 0;
job->sieveArea = 0.0;
job->FB.rfb_size = job->FB.afb_size = -1;
job->FB.rLim = job->FB.aLim = -1;
job->FB.maxP_r = job->FB.maxP_a = -1;
job->FB.maxLP = job->FB.maxLPA = -1;
job->FB.MFB_r = job->FB.MFB_a = -1;
job->rfb_lambda = job->afb_lambda = -1;

while (!(feof(fp))) {
thisLine[0] = 0;
fgets(thisLine, 1023, fp);
if ((sscanf(thisLine, "%255s %255s", token, value)==2) && (thisLine[0] != '#')) {
if (strncmp(token, "a0:",3)==0) {
job->a0 = atol(value);
} else if (strncmp(token, "a1:",3)==0) {
job->a1 = atol(value);
} else if (strncmp(token, "b0:",3)==0) {
job->b0 = atol(value);
} else if (strncmp(token, "b1:",3)==0) {
job->b1 = atol(value);
} else if (strncmp(token, "q0:",3)==0) {
job->qIndex0 = atol(value);
} else if (strncmp(token, "q1:",3)==0) {
job->qIndex1 = atol(value);
} else if (strncmp(token, "sieveArea:",10)==0) {
job->sieveArea = atof(value);
} else if (strncmp(token, "rlambda:",8)==0) {
job->rfb_lambda = atof(value);
} else if (strncmp(token, "alambda:",8)==0) {
job->afb_lambda = atof(value);
} else if (strncmp(token, "fbname:", 7)==0) {
strncpy(job->fbName, value, MAXFNAMESIZE);
} else if (strncmp(token, "rfb_size:", 9)==0) {
job->FB.rfb_size = atol(value);
} else if (strncmp(token, "afb_size:", 9)==0) {
job->FB.afb_size = atol(value);
} else if (strncmp(token, "rlim:", 5)==0) {
job->FB.rLim = atol(value);
} else if (strncmp(token, "alim:", 5)==0) {
job->FB.aLim = atol(value);
} else if (strncmp(token, "lpbr:", 5)==0) {
t = atoi(value); t = MIN(t, 31);
job->FB.maxP_r = (1<<t);
} else if (strncmp(token, "lpba:", 5)==0) {
t = atoi(value); t = MIN(t, 31);
job->FB.maxP_a = (1<<t);
} else if (strncmp(token, "mfbr:", 5)==0) {
job->FB.MFB_r = atoi(value);
} else if (strncmp(token, "mfba:", 5)==0) {
job->FB.MFB_a = atoi(value);
} else if (strncmp(token, "npr:", 4)==0) {
job->FB.maxLP = atoi(value);
} else if (strncmp(token, "npa:", 4)==0) {
job->FB.maxLPA = atoi(value);
}
}
}

if (job->FB.maxLP==-1)
job->FB.maxLP = 3;
if (job->FB.maxLPA==-1)
job->FB.maxLPA = 3;

fclose(fp);
return 0;
}

/***************************************************/
int tryMakeFB(nfs_sieve_job_t *J)
/***************************************************/
/* Try to create a factor base from the info in J. */
/* Return 0 on success, nonzero on failure. */
/***************************************************/
{ int res;

res = createFB(&J->FB, J->fbName);
return res;
}

/**************************************************/
void convertRelations(char *fName, nfs_sieve_job_t *job,
s32 alg_factors_first)
{
FILE *relfile, *outfile;
nfs_fb_t *FB = &job->FB;
mpz_t rRes, aRes;
char buf[512];
s32 rbound = FB->rLim + 1;
s32 abound = FB->aLim + 1;
s32 num_bounds = 0;
s32 b0=0;
s32 pFacts[1000];
s32 numpFacts;

outfile = fopen("spairs.add", "w+a");
if (outfile==NULL) {
printf("cannot open output file\n");
exit(-1);
}
relfile = fopen(fName, "r");
if (relfile==NULL) {
printf("cannot open relation file %s\n", fName);
exit(-1);
}
mpz_init(rRes);
mpz_init(aRes);

fprintf(outfile, "# relations converted from CWI format\n");

fgets(buf, sizeof(buf), relfile);
while (!feof(relfile)) {
s32 rfactors[64];
s32 afactors[64];
s32 rnum, anum;
long long a;
s32 b, i, last_p;
char *tmp, *next_field;

if (buf[0] == '#') {
tmp = strstr(buf, "Prime divisors over ");
if (tmp) {
if (num_bounds == 0) {

if (alg_factors_first)
abound = strtol(tmp+20, NULL, 10);
else
rbound = strtol(tmp+20, NULL, 10);

num_bounds++;
}
else {

if (alg_factors_first)
rbound = strtol(tmp+20, NULL, 10);
else
abound = strtol(tmp+20, NULL, 10);
}
}
fgets(buf, sizeof(buf), relfile);
continue;
}

if (alg_factors_first) {
anum = buf[2];
rnum = buf[3];
}
else {
rnum = buf[2];
anum = buf[3];
}

if (isdigit(rnum)) {
rnum = rnum - '0';
}
else if (isalpha(rnum) && isupper(rnum)) {
rnum = rnum - 'A' + 10;
}
else if (isalpha(rnum)) {
rnum = rnum - 'a' + 36;
}
else {
printf("cannot get number of rational factors\n");
exit(-1);
}

if (isdigit(anum)) {
anum = anum - '0';
}
else if (isalpha(anum) && isupper(anum)) {
anum = anum - 'A' + 10;
}
else if (isalpha(anum)) {
anum = anum - 'a' + 36;
}
else {
printf("cannot get number of algebraic factors\n");
exit(-1);
}

a = strtoll(buf+5, &next_field, 10);
tmp = next_field + 1;
b = strtol(tmp, &next_field, 10);
tmp = next_field + 1;

if (b == 0) {
fgets(buf, sizeof(buf), relfile);
b0++;
continue;
}

if (alg_factors_first) {
for (i = 0; i < anum; i++) {
afactors[i] = strtol(tmp, &next_field, 10);
tmp = next_field + 1;
}
for (i = 0; i < rnum; i++) {
rfactors[i] = strtol(tmp, &next_field, 10);
tmp = next_field + 1;
}
}
else {
for (i = 0; i < rnum; i++) {
rfactors[i] = strtol(tmp, &next_field, 10);
tmp = next_field + 1;
}
for (i = 0; i < anum; i++) {
afactors[i] = strtol(tmp, &next_field, 10);
tmp = next_field + 1;
}
}

if (rnum>0) rbound = MIN(rbound, rfactors[rnum-1]);
if (anum>0) abound = MIN(abound, afactors[anum-1]);

mpz_mul_si(rRes, FB->y1, a);
mpz_addmul_ui(rRes, FB->y0, b);
mpz_abs(rRes, rRes);

mpz_evalF(aRes, a, b, FB->f);
mpz_abs(aRes, aRes);

for (i=0; i<rnum; i++) {
if (!mpz_divisible_ui_p(rRes, rfactors[i])) {
printf("A rational non-factor was found! Perhaps use -a1 or -r1. \n");
} else
mpz_divexact_ui(rRes, rRes, rfactors[i]);
}

for (i=0; i<anum; i++) {
if (!mpz_divisible_ui_p(aRes, afactors[i])) {
printf("An algebraic non-factor was found! Perhaps use -a1 or -r1. \n");
} else
mpz_divexact_ui(aRes, aRes, afactors[i]);
}

for (i = 0; mpz_cmp_ui(rRes,1) && i < 5000 && FB->rfb[2*i] < rbound; i++) {
if (mpz_divisible_ui_p(rRes, FB->rfb[2*i])) {
do mpz_divexact_ui(rRes, rRes, FB->rfb[2*i]); while (mpz_divisible_ui_p(rRes, FB->rfb[2*i]));
rfactors[rnum++] = FB->rfb[2*i];

if ( mpz_probab_prime_p(rRes,6)) {
rfactors[rnum++] = mpz_get_si(rRes);
mpz_set_ui(rRes,1);
}
}
}

if (mpz_cmp_ui(rRes,1)) {
numpFacts = factor(pFacts, rRes, 0);
if (numpFacts>0) {
for (i=0; i<numpFacts; i++) rfactors[rnum++] = pFacts[i];
} else printf("Problem with factor()!\n");
}

last_p = -1;
for (i = 0; mpz_cmp_ui(aRes,1) && i < 5000 && FB->afb[2*i] < abound; i++) {
if (FB->afb[2*i] != last_p &&
mpz_divisible_ui_p(aRes, FB->afb[2*i])) {
do mpz_divexact_ui(aRes,aRes,FB->afb[2*i]); while (mpz_divisible_ui_p(aRes, FB->afb[2*i]));
afactors[anum++] = FB->afb[2*i];

if ( mpz_probab_prime_p(aRes,6)) {
afactors[anum++] = mpz_get_si(aRes);
mpz_set_ui(aRes,1);
last_p = -1;
}
}
last_p = FB->afb[2*i];
}

if (mpz_cmp_ui(aRes,1)) {
numpFacts = factor(pFacts, aRes, 0);
if (numpFacts>0) {
for (i=0; i<numpFacts; i++) afactors[anum++] = pFacts[i];
} else printf("Problem with factor()!\n");
}

//if (mpz_get_ui(aRes)>1) printf("Algebraic residue remaining: %lu\n",mpz_get_ui(aRes));
//if (mpz_get_ui(rRes)>1) printf("Rational residue remaining: %lu\n",mpz_get_ui(rRes));

fprintf(outfile, "%d,%d", a, b);
for (i = 0; i < rnum; i++) {
if (i == 0)
fprintf(outfile, ":%lx", rfactors[i]);
else
fprintf(outfile, ",%lx", rfactors[i]);
}
for (i = 0; i < anum; i++) {
if (i == 0)
fprintf(outfile, ":%lx", afactors[i]);
else
fprintf(outfile, ",%lx", afactors[i]);
}
fprintf(outfile, "\n");

fgets(buf, sizeof(buf), relfile);
}

fclose(relfile);
fclose(outfile);
printf("Ignored %d relations with b=0.\n",b0);
}

/**************************************************/
int main(int argC, char *args[])
{ char fbname[MAXFNAMESIZE];
char jobfile[MAXFNAMESIZE];
char relfile[MAXFNAMESIZE];
int i, fbLoaded=0;
nfs_sieve_job_t Job;
s32 alg_factors_first = 1;

fbname[0]=jobfile[0]=relfile[0]=0;
strncpy(Job.fbName, DEFAULT_FB_NAME, MAXFNAMESIZE);

/***************************/
/* Parse command-line args */
/***************************/
for (i=1; i<argC; i++) {
if (strcmp(args[i], "-fb")==0) {
if ((++i)<argC) {
strcpy(fbname, args[i]);
}
}
else if (strcmp(args[i], "-j")==0) {
if ((++i)<argC) {
strcpy(jobfile, args[i]);
}
}
else if (strcmp(args[i], "-i")==0) {
if ((++i)<argC) {
strcpy(relfile, args[i]);
}
}
else if (strcmp(args[i], "-a1")==0) {
alg_factors_first = 1;
}
else if (strcmp(args[i], "-r1")==0) {
alg_factors_first = 0;
}
else if (strcmp(args[i], "--help")==0) {
printf("USAGE: %s %s\n", args[0], USAGE);
return 0;
}
}
if (!(jobfile[0]) || !(relfile[0])) {
printf("USAGE: %s %s\n", args[0], USAGE);
return 0;
}

/***********************************/
/* Initialize and get factor base. */
/***********************************/
initFB(&Job.FB);

/******************/
/* Get parameters */
/******************/
if (parseJobFile(&Job, jobfile)) {
printf("Some error occured parsing job file. Cannot continue.\n");
exit(-1);
}
if (loadFB(fbname, &Job.FB)) {
fbLoaded=0;
fprintf(stderr, "Could not load factor base from %s. "
"Attempting to create...\n", fbname);
if (tryMakeFB(&Job)) {
fprintf(stderr, "Factor base creation failed! Aborting.\n");
exit(-1);
}
}

convertRelations(relfile, &Job, alg_factors_first);
return 0;
}
[/CODE]