/* ppmpqs without maxbig, based on lip.c */ /* mpqs restartfile [infile logfile [outputdirectory/ [mail-cmd]]] */ #include #include #include #include #include #include "lip.h" #if defined(sparc) && defined(__svr4__) #include #endif /* SOLARIS */ #ifdef NO_SQRTKN #include /* need SEEK_SET */ #endif #define NBITS (30) #define RADIX (1< besttryval) { besttry = try; besttryval = tryval; } } return(besttry); } long root(u,p,lit) /* returns square root of u modulo prime p */ /* dirty, could make this more efficient, needed only in initialization */ long u,p,lit; { long e, q, n, v, w; if (u <= 1) return(u); if (lit) { /* separate code for prime powers (smallprimes) */ for (n=0; n>= 1; } /* p = 1 + q*(2^(e+1)), q odd */ w = zexpmods(u, (q-1)/2, p); v = zmulmods(u, w, p); /* u^((q+1)/2) */ if (e != 0) { long r, yroots[NBITS]; /* yroots[k] = (2^k) root of -1 */ w = zmulmods(v, w, p); /* u^q */ for (n=5; ; n+=6) { long k; if (n >= 1000) { printf("root - can't find non-residue for %ld \n", p); exit(0); } if (zjacobis(n, p) != -1) continue; yroots[e] = zexpmods(n, q, p); for (k=e; k>0; k--) { yroots[k-1] = zmulmods(yroots[k], yroots[k], p); } if (yroots[0] == p-1) break; } r = e+1; while (w != 1) { /* v^2 == u*w, w^(2^(r-1)) == 1 */ long k = 1, z = w; while (z != p-1) { z = zmulmods(z, z, p); if (k++ >= r) { printf("Error in root, perhaps a non-residue:\n"); printf("p=%ld, u=%ld, v=%ld, w=%ld, r=%ld, z=%ld, n=%ld\n", p, u, v, w, r, z, n); exit(0); } } v = zmulmods(v, yroots[k], p); w = zmulmods(w, yroots[k-1], p); r=k; /* printf("root - u = %ld, p=%ld, v=%ld, w=%ld, r=%ld\n", u, p, v, w, r); */ } /* while w */ } /* if e */ if (v > p/2) v = p-v; if (zmulmods(v,v,p) != u) { printf("Error in root: %ld %ld %ld \n", u, p, v); exit(0); } return(v); } createfactorbase() { long lit, p, knmodp, i = 1; uchar addlog; long exp_addlog; zpstart(); /* so that next prime will be 3 */ if ((kn[1]%8) == 1) { primbase[1] = 2; sqrtkn[1] = root(kn[1],smallprimes[2],(long)1); i = 2; nbsmallps = 1; } else nbsmallps = 0; while (i <= sizefacbase) { lit=0; p = zpnext(); knmodp = zsmod(kn, p); if (zjacobis(knmodp, p) != -1 || p == multiplier) { primbase[i] = p; if (p <= MAXSMALLPRIME) { lit =1; nbsmallps ++; p = smallprimes[p]; knmodp = zsmod(kn, p); } sqrtkn[i] = root(knmodp,p,lit); i++; } } nbsmallpsp = nbsmallps + 1; lastprime = primbase[sizefacbase]; for (i=1; i<=nbsmallps; i++) smlogp[i]=(log((double)primbase[i])/1.0); logcnt = 0; if (i<=sizefacbase) { firstlog = (uchar)(log((double)primbase[i])/1.0); addlog = firstlog; exp_addlog = (long)(exp((double)(addlog+1))/1.0); } i++; for (; i<=sizefacbase; i++) { if (primbase[i] > exp_addlog ) { nextlog[logcnt++] = primbase[i]; addlog ++; /* step > 1 possible error, but unlikely */ /* to affect many logs, at most a few */ exp_addlog = (long)(exp((double)(addlog+1))/1.0); } } nextlog[logcnt] = lastprime+1; } void findsievebound() /* better to determine experimentally */ /* this one gives a rather low result (which is better, */ /* because it gives more information. In the course of the */ /* computation the bound should be adapted. */ { double lkn, loglast; lkn = log((double)kn[kn[0]]) + (double)(kn[0]-1)*log((double)RADIX); loglast = 2 * log((double)maxbig); sievebound=(uchar)(((lkn+log((double)sievelength))/2.0-loglast)/1.0-4.0); smallsievebnd= sievebound - (uchar)SMALLCORRECT; } long dumpdata() { /* dump in decimal notation, and not using fwrite, so */ /* that different machines can read each others restart file */ /* slow, but saves more time that it wastes */ /* put on separate lines, to make it easier to change RESTART */ /* using some editor (long line can cause problems) */ register long i; if (fprintf(restart,"%ld",sizefacbase) < 0) return(1); if (fprintf(restart," %ld\n",identification) < 0) return(1); zfwriteln(restart,n); if (fprintf(restart,"%ld\n",multiplier) < 0) return(1); zfwriteln(restart,sqkn); if (fprintf(restart,"%ld",sievelength) < 0) return(1); if (fprintf(restart," %ld",(long)sievebound) < 0) return(1); if (fprintf(restart," %ld\n",maxbig) < 0) return(1); if (fprintf(restart," %ld\n",nbsmallps) < 0) return(1); for (i=1; i<=nbsmallps; i++) { if (fprintf(restart," %ld",(long)smlogp[i]) < 0) return(1); if (i%5 == 0) fprintf(restart,"\n"); } if (fprintf(restart," %ld\n",(long)firstlog) < 0) return(1); if (fprintf(restart," %ld\n",logcnt) < 0) return(1); for (i=0; i<=logcnt; i++) { if (fprintf(restart," %ld",nextlog[i]) < 0) return(1); if (i%5 == 0) fprintf(restart,"\n"); } if (i%5 != 1) fprintf(restart,"\n"); zfwriteln(restart,startd); for (i=1; i<=sizefacbase; i++) { #ifdef USEHEX long gap = (i <= 2) ? primbase[i] : ((primbase[i]-primbase[i-1]) >> 1); #else long gap = primbase[i]; #endif if (fprintf(restart," %ld",gap) < 0) return(1); if (i%5 == 0) fprintf(restart,"\n"); } for (i=1; i<=sizefacbase; i++) { if (fprintf(restart,SPNUM,sqrtkn[i]) < 0) return(1); if (i%5 == 0) fprintf(restart,"\n"); } fflush(restart); #ifdef NO_SQRTKN free (sqrtkn); #else fclose(restart); #endif return(0); } makestartd() { long zcompare(); if (unique < 0) { zintoz(-unique,&altstartd); zlshift(altstartd,NBITS,&altstartd); if (zcompare(altstartd,startd) >= 0) { printf("negative unique process ID not allowed, exit\n"); exit(0); } zsub(startd,altstartd,&startd); startd[1] |= 3; } else { zintoz(unique,&altstartd); zlshift(altstartd,NBITS,&altstartd); zadd(altstartd,startd,&startd); startd[1] |= 3; } /* no overlap is almost guaranteed in this way, unless machines are */ /* VERY fast, and run for a LONG time. resulting values are close */ /* enough to optimal for our range of interest (10^90 - 10^110) */ } getmemory() { long bytefacbase = (sizefacbase+1) * sizeof(long); if (((primbase = (long*)malloc(bytefacbase))==NULL)|| ((sqrtkn=(long*)malloc(bytefacbase))==NULL) || ((roots=(roots_t*)malloc((sizefacbase+1)*sizeof(roots_t)))==NULL)) { printf("unable to get memory, exit\n"); exit(0); } notallocated = 0; #if 0 printf("primbase = %lx, sqrtkn = %lx, roots = %lx\n", (long)primbase, (long)sqrtkn, (long)roots); #endif } timestamp () { long tp = time(NULL); timenow = getutime(); printf(" after %10.2lf seconds at %s",timenow-timestart,ctime(&tp)); fflush(stdout); } long getandalloc() { register long i; long dummy; if (how == SLAVE) { printf("Attempt to read restart file %s", restartfn); timestamp(); if ((restart=fopen(restartfn,"r"))==NULL) return(1); if (fscanf(restart,"%ld",&sizefacbase) == EOF) goto wrong; if (fscanf(restart,"%ld",&dummy) == EOF) goto wrong; if (dummy != identification) goto wrong; zfread(restart,&dum); if (zcompare(n,dum) != 0) goto wrong; } else { printf("impossible call to getandalloc, exit"); timestamp(); exit(0); } /*successfully opened file, and contains the right n */ getmemory(); if (fscanf(restart,"%ld",&multiplier)==EOF) goto wrong; zfread(restart,&sqkn); if (fscanf(restart,"%ld",&sievelength)==EOF) goto wrong; if (fscanf(restart,"%ld",&dummy)==EOF) goto wrong; sievebound = (uchar)dummy; smallsievebnd= sievebound - (uchar)SMALLCORRECT; if (fscanf(restart,"%ld",&maxbig)==EOF) goto wrong; if (fscanf(restart,"%ld\n",&nbsmallps)==EOF) goto wrong; nbsmallpsp= nbsmallps+1; for (i=1; i<=nbsmallps; i++) { if (fscanf(restart,"%ld",&dummy)==EOF) goto wrong; } if (fscanf(restart,"%ld\n",&dummy)==EOF) goto wrong; firstlog = (uchar)dummy; if (fscanf(restart,"%ld\n",&logcnt)==EOF) goto wrong; for (i=0; i<=logcnt; i++) { if (fscanf(restart,"%ld",&nextlog[i])==EOF) goto wrong; } zfread(restart,&startd); for (i=1;i<=sizefacbase;i++) { if (fscanf(restart,"%ld",&primbase[i]) == EOF) goto wrong; #ifdef USEHEX if (i > 2) primbase[i] = primbase[i-1] + 2*primbase[i]; #endif } #ifdef NO_SQRTKN sqrtknbase = ftell (restart); #else for (i=1;i<=sizefacbase;i++) if (fscanf(restart,NOSPNUM,&sqrtkn[i]) == EOF) goto wrong; fclose(restart); #endif for (i=1; i<=nbsmallps; i++) smlogp[i]=(log((double)primbase[i])/1.0); lastprime = primbase[sizefacbase]; zsmul(n,multiplier,&kn); if ((multiplier>1) && (multiplier <= MAXSMALLPRIME)) smallprimes[multiplier] = multiplier; makestartd(); return(0); wrong: fclose(restart); if (how == SLAVE) { unlink(restartfn); printf("Attempt failed, restart file removed, creating %s", restartfn); timestamp(); } else printf("Exit: unsuccessful restart"); timestamp(); return(1); } setupmailer() { scanf("%ld",&tobemailed); printf("the results will be written to the file %s\n",nameoutput); if (tobemailed) { sprintf(callsendit,mail_cmd, nameoutput); printf("As soon as it contains %ld relations,\n",MAILBATCH); printf("this file will be mailed using the command:\n"); printf(" %s\n",callsendit); } else { printf("Please process %s every now and then:\n",nameoutput); #ifdef VMS printf(" $ rename %s MAILOUT\n", nameoutput); printf(" (wait a moment, maybe process is writing to MAILOUT)\n $ "); printf(mail_cmd, "MAILOUT"); printf ("\n"); #else printf(" mv %s MAILOUT\n",nameoutput); printf(" (wait a moment, maybe process is writing to MAILOUT)\n "); printf(mail_cmd, "MAILOUT\n"); #endif } fflush(stdout); } setup () { long bytefacbase, dummy; #if BUILD_RESTART char restartname[150]; #endif zstart(); /* set up long integer arithmetic */ if (how != STARTER) { scanf("%ld %ld",&identification,&unique); sprintf(parid,"par%ld#",identification); sprintf(pprid,"ppr%ld#",identification); sprintf(fulid,"ful%ld#",identification); sprintf(paridc,"par%ldc#",identification); sprintf(ppridc,"ppr%ldc#",identification); sprintf(fulidc,"ful%ldc#",identification); sprintf(nameoutput,"%sOUT%07ld",file_prefix,unique); zrstarts(abs(unique)+1); } if (how == SLAVE) { zread(&n); if (getandalloc() == 0) { printf("successful restart from file %s", restartfn); timestamp(); scanf("%ld",&dummy); /* multiplier does not change */ scanf("%ld",&dummy); /* sizefacbase does not change */ scanf("%ld",&sievelength); /* sievelength may be changed */ scanf("%ld",&dummy); sievebound = (uchar)dummy; /* sievebound may be changed */ scanf("%ld",&maxbig); /* maxbig may be changed */ zread(&startd); makestartd(); } else { scanf("%ld %ld",&multiplier,&sizefacbase); zsmul(n,multiplier,&kn); zsqrt(kn,&sqkn,&dum); if ((multiplier>1) && (multiplier <= MAXSMALLPRIME)) smallprimes[multiplier] = multiplier; if (notallocated) getmemory(); scanf("%ld %ld %ld",&sievelength,&dummy,&maxbig); sievebound = (uchar)dummy; printf("start creation new restart data"); timestamp(); createfactorbase(); zread(&startd); #if BUILD_RESTART #ifdef VMS /* VMS has versioning, don't need to make unique */ strcpy(restartname, restartfn); #else sprintf(restartname,"%s%07ld",restartfn, unique); #endif /*VMS*/ if ((restart = fopen(restartname,"w")) == NULL) { printf("Cannot create restart file %s, non fatal\n", restartname); } else { if (dumpdata()) { printf( "Cannot dump restart data, non fatal\n"); unlink(restartname); } else { #ifndef VMS link(restartname, restartfn); unlink(restartname); #endif /*VMS */ printf( "Created restart file %s", restartfn); timestamp(); } } fflush(stdout); #endif /* BUILD_RESTART */ makestartd(); } smallsievebnd= sievebound - (uchar)SMALLCORRECT; printf("multiplier: %ld\n",multiplier); printf("sizefacbase: %ld\n",sizefacbase); printf("sievelength: %ld\n",sievelength); printf("sievebound: %ld\n",(long)sievebound); printf("maxbig: %ld\n",maxbig); printf("startd:"); zwrite(startd); printf("\n"); fflush(stdout); setupmailer(); started = 1; return; } if (how == STARTER) { /*includes FIREFLYSTART */ /* reads n, n's ID, sizefacbase, and maxbig */ /* builds RESTART, prints restartvalues, and stops */ printf("start up, create restart file in this directory\n"); printf("give following data:\n"); printf(" number to be factored,\n"); printf(" its unique identification,\n"); printf(" size of factor base,\n"); printf(" and bound for largest prime in partials.\n"); fflush(stdout); zread(&n); scanf("%ld",&identification); if ((n[1] & 1) == 0) { printf("input is even, exit\n"); exit(0); } scanf("%ld %ld",&sizefacbase,&maxbig); printf("begin computing restart data for factorization of\n"); zwrite(n); printf("\n"); printf("identification: %ld\n",identification); printf("size factor base: %ld\n",sizefacbase); printf("primes in partial factorizations at most: %ld\n",maxbig); fflush(stdout); multiplier = findmultiplier((long)100); printf("multiplier: %ld\n",multiplier); fflush(stdout); zsmul(n,multiplier,&kn); zsqrt(kn,&sqkn,&dum); if ((kn[1] % 4) != 1) { printf("wrong multiplier selected, exit\n"); exit(0); } if ((multiplier>1) && (multiplier <= MAXSMALLPRIME)) smallprimes[multiplier] = multiplier; bytefacbase = (sizefacbase+1) * sizeof(long); if (((primbase = (long*)malloc(bytefacbase)) == NULL) || ((sqrtkn = (long*)malloc(bytefacbase)) == NULL)) { printf("unable to get memory, exit\n"); exit(0); } createfactorbase(); /* puts data in primbase, r1, and log info */ /* also sets lastprime */ printf("factor base created\n"); fflush(stdout); if ((dummy = smalldiv()) > 0 ) { printf("small factor: %ld, exit\n",dummy); exit(0); } printf("largest prime in factor base: %ld\n",lastprime); #ifdef SIEVEL if ((sievelength = SIEVEL*lastprime) > 400) #else if ((sievelength = 2*lastprime) > 400) /* results in 4*lastprime */ #endif sievelength = 400*((sievelength+200)/400); printf("sieve length: %ld\n",sievelength); findsievebound(); printf("sievebound: %ld\n",(long)sievebound); fflush(stdout); z2div(kn,&dum); zsqrt(dum,&startd,&dum); zintoz(sievelength,&dum); zdiv(startd,dum,&dum,&startd); zsqrt(dum,&startd,&dum); startd[1] |= 3; if ((startd[0] > 1 ) || (startd[1] > 4)) zsadd(startd,(long)-4,&startd); printf("polynomial generator starts at: "); zwriteln(startd); fflush(stdout); if ((restart = fopen(restartfn,"w")) == NULL) { printf("Could not open restart file, exit\n"); exit(0); } if (dumpdata()) { printf("Could not write restart file, exit\n"); unlink(restartfn); exit(0); } printf("restart file created, done\n"); printf("\n"); printf("Give each machine the following input:\n"); printf("--------- ><8 --------- \n"); printf("%ld %ld i\n",SLAVE,identification); zwriteln(n); printf("%ld %ld %ld %ld %ld\n",multiplier,sizefacbase, sievelength,(long)sievebound,maxbig); zwriteln(startd); printf("0/1\n"); printf("address\n"); printf("\n"); printf("--------- ><8 --------- \n"); printf( "The second to last line contains a zero or a one; a one if\n"); printf("the results should be mailed to address (last line), a\n"); printf("zero if the owner of that machine will mail the\n"); printf( "results. The mpqs-process will create a file RESTART which\n"); printf("can be used by other processes running in the same\n"); printf("directory; this file will not be created if it already\n"); printf("exists for this n (with the same identification).\n"); printf(" The i on the first line should be an integer and at\n"); printf( "least zero. It should be different from any i given to other\n"); printf("processes. Never use any i more than once, even on the\n"); printf("same machine, to avoid duplicating curves/relations.\n"); exit(0); } /* how is not STARTER, SLAVE or FIREFLY: */ printf("wrong restart value, exit\n"); exit(0); } findnextpolynomial(a,b,d2inv) verylong *a, *b, *d2inv; { static verylong t1=0, t2=0, t3=0, t2kn=0, t4=0, ro1=0, ro2=0, pp=0; register long i, p, ap, bp, bpc; long rem1, mod1 = 1033715657; /* 13*23*37*41*43*53 */ long rem2, mod2 = 1035563755; /* 5*17*29*61*67*97 */ long rem3, mod3 = 1036405977; /* 3*19*47*59*79*83 */ long rem4, mod4 = 1039058713; /* 7*11*31*67*73*89 */ i = 1; zcopy(startd,&altstartd); interrupted = 1; startd[1] |= 3; /* to be sure */ rem1 = zsmod(startd, mod1); rem2 = zsmod(startd, mod2); rem3 = zsmod(startd, mod3); rem4 = zsmod(startd, mod4); while (i > 0) { zsadd(startd,(long)4,&startd); rem1 += 4; rem2 += 4; rem3 += 4; rem4 += 4; if ( zjacobis(rem1,mod1) == 0 || zjacobis(rem2,mod2) == 0 || zjacobis(rem3,mod3) == 0 || zjacobis(rem4,mod4) == 0) continue; /* Reject d if divisible by odd prime < 100 */ zrshift(startd,(long)2,&t1); /* t1 = (d-3)/4 */ zexpmod(kn,t1,startd,&t2); /* t2 = 1/sqrt(kn) mod d if d good */ zmul(t2,kn,&t2kn); /* t2kn = t2*kn = kn^((d+1)/4) mod d */ zmul(t2,t2kn,&t3); zdiv(t3,startd,&t4,&t1); /* kn^((d-1)/2) mod d */ if ((t1[0] != 1) || (t1[1] != 1)) continue; /* nonresidue or nonprime*/ if (zprime(startd, (long)1, 2) > 0) i = 0; /* insurance */ } interrupted = 0; zmod(t4,kn,d2inv); if ((*d2inv)[1]&1) zadd(*d2inv, kn, d2inv); z2div(*d2inv,d2inv); zsub(kn,*d2inv,d2inv); /* d2inv = 1/2d mod kn */ zsq(startd,a); zmod(t2kn,*a,&t2kn); zmod(t3,*a,&t3); /* Compute b = t2*kn*(t2^2*kn - 3)/2 = sqrt(kn) mod a */ zsadd(t3,(long)-3,&t3); zmulmod(t2kn,t3,*a,b); if ((*b)[1] & 1) zadd(*a,*b,b); z2div(*b,b); if (((*b)[1] & 1) == 0) zsub(*a,*b,b); /* Odd square root */ zadd(sqkn,*b,&t1); zdiv(t1,*a,&t1,&t2); z2div(t1,&t1); roots[0].r1 = -t1[1]; zsub(sqkn,*b,&t1); zdiv(t1,*a,&t1,&t2); z2div(t1,&t1); roots[0].r2 = t1[1]; #ifdef NO_SQRTKN fseek (restart, sqrtknbase, SEEK_SET); /* Skip over beginning of table. */ for (i=1; i < nbsmallpsp; i++) fscanf(restart,NOSPNUM,&sqrtkni); #endif for (i=nbsmallpsp; i<=sizefacbase; i++) { #ifdef NO_SQRTKN fscanf(restart,NOSPNUM,&sqrtkni); #endif zintoz( (p = primbase[i]), &pp); bp = zsmod(*b,p); ap = zsmod(startd,p); ap = zmulmods(ap,ap,p); /* a mod p (a = d^2) */ if (ap != 0) { bpc = bp; #ifdef NO_SQRTKN if ((bp += sqrtkni) >= p) bp -= p; #else if ((bp += sqrtkn[i]) >= p) bp -= p; #endif ap = zinvodds(2*ap,p); roots[i].r1 = zmulmods(p-ap, bp, p); #ifdef NO_SQRTKN if ((bp = sqrtkni-bpc) < 0) bp += p; #else if ((bp = sqrtkn[i]-bpc) < 0) bp += p; #endif roots[i].r2 = zmulmods(ap, bp, p); } else { zsq(*b,&t1); if (zcompare(kn,t1) >= 0) zsub(kn,t1,&t1); else { zsub(t1,kn,&t1); bp = p-bp; } zrshift(t1,(long)2,&t1); zdiv(t1,*a,&t2,&t1); zintoz(zsmod(t2,p),&ro1); zintoz(zinvodds(bp,p),&ro2); zmulmod(ro1,ro2,pp,&ro1); roots[i].r2 = (roots[i].r1 = ro1[1]); } } } checkroot(a,b) /* This one could be done much nicer, but it is not important */ /* if CHECK has a reasonably large value */ verylong a, b; { long i, misser, p, minc; static verylong r=0, c=0; misser = 0; #ifdef NO_SQRTKN fseek (restart, sqrtknbase, SEEK_SET); #endif for (i=1; i<=sizefacbase; i++) { p = primbase[i]; if (p < MAXSMALLPRIME) p = smallprimes[p]; #ifdef NO_SQRTKN fscanf(restart,NOSPNUM,&sqrtkni); if (zsmod(kn,p) != zmulmods(sqrtkni, sqrtkni, p)) #else if (zsmod(kn,p) != zmulmods(sqrtkn[i], sqrtkn[i], p)) #endif { printf("Wrong at primbase[%ld]=%ld\n",i,p); fflush(stdout); misser ++; } } if (misser == 0) printf("squareroots correct\n"); else printf("%ld mistakes in the squareroots\n",misser); fflush(stdout); zsq(b,&r); if (zcompare(kn,r) >= 0) { zsub(kn,r,&r); minc = -1; } else { zsub(r,kn,&r); minc = 1; printf("b is very large!\n"); } if ((r[1]&3) != 0) { printf("kn-b^2 is not 0 modulo 4, error, exit\n"); exit(0); } zrshift(r,(long)2,&r); zdiv(r,a,&c,&r); if ((r[0] != 1 ) || (r[1] != 0)) { printf("kn-^2 not divisible by a, error, exit\n"); exit(0); } misser = 0; for (i=nbsmallpsp; i<=sizefacbase; i++) { long ap, bp, cp, res1, res2; long root1 = roots[i].r1, root2 = roots[i].r2; p = primbase[i]; if (p <= MAXSMALLPRIME) p = smallprimes[p]; ap = zsmod(a,p); bp = zsmod(b,p); cp = minc*zsmod(c,p); res1 = zmulmods(ap, root1, p) + bp; res2 = zmulmods(ap, root2, p) + bp; if (res1 >= p) res1 -= p; if (res2 >= p) res2 -= p; res1 = zmulmods(res1, root1, p) + cp; res2 = zmulmods(res2, root2, p) + cp; if (res1 != 0 && res1 != p) { misser ++; printf("wrong first root at i=%ld p=%ld\n",i,p); printf("ap=%ld, bp=%ld, cp=%ld, root=%ld\n", ap,bp,cp,root1); #ifdef NO_SQRTKN printf("sqrtkn mod p = %ld, (2a)^-1 mod p = %ld\n",sqrtkni, zinvodds(2*ap,p)); #else printf("sqrtkn mod p = %ld, (2a)^-1 mod p = %ld\n",sqrtkn[i], zinvodds(2*ap,p)); #endif fflush(stdout); } if (root1 == root2) printf("double root at p=%ld\n",p); if (res2 != 0 && res2 != p) { misser ++; printf("wrong second root at i=%ld p=%ld\n",i,p); printf("ap=%ld, bp=%ld, cp=%ld, root=%ld\n",ap,bp,cp,root2); #ifdef NO_SQRTKN printf("sqrtkn mod p = %ld, (2a)^-1 mod p = %ld\n", sqrtkni,zinvodds(2*ap,p)); #else printf("sqrtkn mod p = %ld, (2a)^-1 mod p = %ld\n", sqrtkn[i],zinvodds(2*ap,p)); #endif fflush(stdout); } } if (misser == 0) printf("roots correct\n"); else { printf("%ld mistakes in the roots, exit\n",misser); exit(0); } fflush(stdout); } #ifndef RISC #if !(ASM_sieveinit) sieveinit(lng, inval) /* use that sizeof(long) = 2^LOGLONGSIZE*/ long lng, inval; { register long val = inval, i = (lng + (1<> LOGLONGSIZE, *psiev = (long*)sieve; while (i >= 100) { /* Do 100 longwords at at time */ i -= 100; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; *psiev++ = val; } while (i > 0) {*psiev++ = val; i--;} } #endif /*not RISC*/ #endif /* ASM_sieveinit */ dumpfac(w,wid,widc,ind,fac,bp,bp2,qi) long *w, ind; verylong fac; long bp, bp2; verylong qi; char wid[], widc[]; { register long i; int status; if ((mailit = fopen(nameoutput,"a+")) == NULL) { printf("Cannot open file %s.\n",nameoutput); printf("results to standard output\n"); mailit = stdout; closeit = 0; } (*w) ++; fprintf(mailit,"%c%c%c %ld from %d\n",wid[0],wid[1],wid[2],*w,unique); fprintf(mailit,"%s%ld %ld %ld",wid,bp,bp2,ind); for (i=1;i<=ind;i++) { if (i%10 ==0) fprintf(mailit,"\n%s",widc); fprintf(mailit," %ld",fac[i]); } fprintf(mailit,"\n%s 0 ",widc); zfwrite(mailit,qi,LINELENGTH,"",widc); fprintf(mailit,"\n"); fflush(mailit); if (closeit) fclose(mailit); if (++outputcounter == MAILBATCH && tobemailed) { outputcounter = 0; printf("executing %s",callsendit); timestamp(); status = system(callsendit); #ifdef VMS if ((status & 01) != 1) #else if (status) #endif printf("mail failed '%s' (status = %d)\n", callsendit, status); else /* delete only if mail worked */ unlink(nameoutput); } numfac += ind/2; } long try_ecm( verylong a, verylong *f1, verylong *f2) { static long *r = 0; ecc ++; if (zecm(a,f1,0,10,100,1,0)) { zdiv(a,*f1,f2,&r); ecs ++; return(1); } return(0); } long is_sq(n) long n; { /* Test whether n is a perfect square. If so, return its square root. If not, return -1. Optimized for case where most arguments will be non-squares. (Peter Montgomery) */ static int first = 1; static unsigned char squtab[667]; /* Quadratic residues for primes <= 53 */ register long nn = n; #define MSK371 1 #define MSK517 2 #define MSK559 4 #define MSK589 8 #define MSK615 16 #define MSK629 32 #define MSK667 64 if (first) { register long i; first = 0; for (i = 0; i < 667; i++) squtab[i] = 0; for (i = 0; i < 334; i++) { register long isq = i*i; squtab[isq%371] |= MSK371; squtab[isq%517] |= MSK517; squtab[isq%559] |= MSK559; squtab[isq%589] |= MSK589; squtab[isq%615] |= MSK615; squtab[isq%629] |= MSK629; squtab[isq%667] |= MSK667; /* Initialize quadratic residues */ } } if (nn < 0) return -1; if (nn < 1) return nn; while ((nn & 3) == 0 ) nn >>= 2; if ( (nn & 7) == 1 && (squtab[nn%371] & MSK371) && (squtab[nn%517] & MSK517) && (squtab[nn%559] & MSK559) && (squtab[nn%589] & MSK589) && (squtab[nn%615] & MSK615) && (squtab[nn%629] & MSK629) && (squtab[nn%667] & MSK667)) { /* Checks modulo all primes <= 53 */ long root = zsqrts(n); if (n == root*root) return root; return -1; } else { return -1; } #undef MSK371 #undef MSK517 #undef MSK559 #undef MSK589 #undef MSK615 #undef MSK629 #undef MSK667 } long squf( verylong n, verylong *f1, verylong *f2) { #define MSD 30 static long *t1=0, *t2=0; register long iter, nsqroot, Qprev, Qnow, Pnow, iterbnd=50000; register long den_bound, nsmallden=0, donethat=0; long smalldens[MSD]; int phase = 1; if (zsqrt(n,&t1,&t2)) { zcopy(t1,f1); zcopy(t1,f2); return(1); } Qnow = t2[1]; nsqroot = t1[1]; Pnow = nsqroot; Qprev = 1; den_bound = zsqrts(2*nsqroot+1); for (iter = 1; iter = 2*Qprev) { long quot = (nsqroot + Pprev)/Qprev; Pnow = quot * Qprev -Pprev; Qnow = Qback2 + quot*(Pprev-Pnow); } else { Pnow = Qprev - Pprev; Qnow = Qback2+Pprev-Pnow; } den = Qprev; if ((den & 1) == 0) den >>= 1; if (Pnow == Pprev) { zintoz(den,f1); if (zsdiv(n,den,f2)) { /* printf("error 3 in squf\n"); */ goto doecm; } return(iter); } if (phase && den1) { if (nsmallden 0)) { register long j = nsmallden; smalldens[0] = root; while (smalldens[j] != root) j--; if (root == 1) { goto doecm; } else if (j == 0 ) { phase = 0; if (donethat) { goto doecm; } iterbnd += iter; donethat = 1; Pnow = nsqroot - (nsqroot-Pnow)%root; zintoz(Pnow,&t1); zsq(t1,&t2); zsub(n,t2,&t2); if (zsdiv(t2,root,&t1)) { /* printf("error in squf\n"); */ goto doecm; } Qnow = t1[1]; Qprev = root; } } } doecm: return(try_ecm(n,f1,f2)); } long factor_attempt( verylong a, verylong *f1, verylong *f2) { register long i; i=squf(a,f1,f2); if (i) { if (zcompare(a,*f1) && zcompare(a,*f2)) { if (zcompare(*f1,*f2)>0) zswap(f1,f2); return(1); } } return(0); } sievget(a,b,d2inv,sign,i,rshift,value) verylong a, b, d2inv; long sign,i,rshift; uchar value; { static verylong qisq=0, t1=0, t2=0, t3=0; long factors[MAXFACNB], index = 1; /* Set P(X) = ((2*a*x + b)^2 - kn)/4*d^2, where a = d^2 (see findnextpolynomial). sievget is called when siever detects that P(i) is highly composite. This routine computes P(i) explicitly and possibly outputs a relation containing qisq = (2*a*i + b)/2d and the factors of P(i) == qisq mod kn. */ closeit = 1; zsmul(a,2*(i > 0 ? i : -i),&t1); if (i > 0) zadd(t1,b,&t1); else if (zcompare(t1,b) >= 0) zsub(t1,b,&t1); else zsub(b,t1,&t1); /* t1 = +- (2*a*i + b) */ zmulmod(t1,d2inv,kn,&qisq); zmulmod(qisq,qisq,kn,&t1); if (i>=roots[0].r1 && i<=roots[0].r2) { factors[1]= -1; zsub(kn,t1,&t1); } else factors[1]=1; { register long j; double contrib = 0.0; for (j=1;j<=nbsmallps;j++) { long pj = primbase[j]; if (!(zsdiv(t1,pj,&t2))) { long cntr = 1; zcopy(t2,&t1); cntr=1; while (!(zsdiv(t1,pj,&t2))) { cntr ++; zcopy(t2,&t1); } contrib += (double)cntr*smlogp[j]; factors[++index] = cntr; factors[++index] = j; } } /* for j */ if ((double)value+contrib < (double)sievebound) return; } { register long j; register long offset = i - rshift; /* negative */ register roots_t *loc_roots = roots; register row loc_primbase = primbase; register long loc_sizefacbase = sizefacbase; for (j = nbsmallpsp; j <= loc_sizefacbase; j++) { register long imodp = offset; register long pj = loc_primbase[j]; /* Test whether roots[j].r1 - offset or roots[j].r1 + roots[j].r2 - offset is divisible by pj. Recall offset < 0 and 0 <= roots[j].r1 <= roots[j].r1 + roots[j].r2 < p . Ensure that first argument to % is nonnegative. */ if (imodp <= -(pj << 3)) { imodp = pj - 1 - (-1 - imodp)%pj; } else { /* quotient at most 8, often much less */ do {} while ((imodp += pj) < 0); } imodp -= loc_roots[j].r1; if (imodp == 0 || imodp == loc_roots[j].r2) { if (!(zsdiv(t1,pj,&t2))) { long cntr = 1; zcopy(t2,&t1); while (!(zsdiv(t1,pj,&t2))) { cntr ++; zcopy(t2,&t1); } if (index >= MAXFACNB-2) { printf("Increase MAXFACNB\n"); return; } factors[++index] = cntr; factors[++index] = j; } else { printf("Sievget - bad root, p = %ld\n", pj); } } /* if imodp */ } /* for j */ } { long ival= (long)value; double t1value; int lucky; if (t1[0] > 2) { lac ++; lucky = 0; /* Cofactor too large */ } else if (t1[0] == 1) { long pj = t1[1]; if (pj == 1) { dumpfac(&ful,fulid,fulidc,index,factors,1L,1L,qisq); if (ival>maxful) maxful = ival; if (ivalmaxpar) maxpar = ival; if (ival1 || t2[0]>1 ) { woc ++; lucky = 0; } else { /* now they are both small enough */ if (t2[1] <= lastprime) { printf("Sievget - missed f = %ld\n", t2[1]); } /* the pp case */ dumpfac(&ppr,pprid,ppridc,index, factors,t2[1],t3[1],qisq); if (ival>maxppr) maxppr = ival; if (ivalmaxbad) maxbad = ival; if (ival= lower_bound. */ /* Should be optimized for case where very few elements are >= lower_bound. */ long start; uchar array[], lower_bound; { register uchar *parray = &array[start]; register uchar bound = lower_bound; if (sizeof(long) == 4*sizeof(uchar) & FIRST_GE_SENTINEL >= 4) { /* Advance to longword boundary (N.B. array was allocated by malloc) */ do { if (*parray >= bound) return (parray-array); } while ((++parray - array)&3); /* Check four bytes at a time. */ /* If something interesting found, complete search one byte at a time. */ /* This attempts to cater to both big-endian and little-endian machines. */ { register long long_bound = (0x01010101)*(bound - 1); register long *plarray = (long*)(parray); while ( (((long_bound - plarray[0]) | (long_bound - plarray[1]) | (long_bound - plarray[2]) | (long_bound - plarray[3])) & 0x80808080) == 0) plarray += 4; parray = (uchar*)(plarray); /* printf("first %ld\n",(long)(parray-1-array)); printf("plarray = %lx, parray = %x %x %x %x, bound=%x\n", plarray[-1], parray[0], parray[1], parray[2], parray[3],bound); fflush(stdout); */ } } while (*parray++ < bound) {}; /* null body */ /* printf("final %ld\n",(long)(parray-1-array)); fflush(stdout); */ return (parray-1-array); } #endif /* ASM_first_ge */ /* sievget_check checks for elements of sieve which are >= smallsievebnd */ #define sievget_check(direction,shifter) {i = -1; \ sieve[rsl] = smallsievebnd; \ while ((i = first_ge(sieve, i+1, smallsievebnd)) < rsl) {\ sievget(a,b,d2inv,(long)(direction),i+ishift,\ (long)(shifter),sieve[i]);}} /* Procedures bigp_siever, lowp_siever, medp_siever all have four arguments (ibeg, iend, rsl, logp). Each sieves for all primes p = primbase[i], where ibeg <= i <= iend, by adding logp to all appropriate entries of the sieve. All assume the new roots are roots[i].r1 and roots[i].r1 + roots[i].r2, where 0 <= roots[i].r1 <= roots[i].r1 + roots[i].r2 <= p-1 The difference is the range of values of p assumed by each procedure: bigp_siever -- p >= rsl (at most one hit per root) lowp_siever -- p <= rsl/2 (at least two hits per root) medp_siever -- rsl/2 <= p <= rsl (one or two hits per root) */ #if !(ASM_bigp_siever) bigp_siever(ibeg, iend, rsl, logp) long ibeg; /* lower bound on i for sieving */ long iend; /* upper bound on i for sieving */ long rsl; /* width of sieve interval */ uchar logp; /* logarithm of these primes */ { register long loc_rsl = rsl, i; register uchar loc_logp = logp; register urow sieve_rsl = sieve + loc_rsl; register roots_t *rootsi = roots + ibeg; register row pi = primbase + ibeg; for (i = ibeg - iend - 1; i < 0; i++) { register long root1 = rootsi->r1 - loc_rsl; /* root relative to end of sieve interval */ #if 0 if ( rootsi->r1< 0 || rootsi->r2 < 0 || rootsi->r1 + rootsi->r2 >= *pi || *pi < loc_rsl) { printf("Bad r1 = %ld, r2 = %ld, bigp = %ld\n", rootsi->r1, rootsi->r2, *pi); } #endif if (root1 < 0) { register long root2 = root1 + rootsi->r2; sieve_rsl[root1] += loc_logp; root1 += *pi; /* Advance smaller root by p */ if (root2 < 0) { sieve_rsl[root2] += loc_logp; /* Difference unchanged */ } else { rootsi->r2 = root1 - root2; /* New difference */ root1 = root2; /* New smaller root */ } } rootsi->r1 = root1; rootsi++; pi++; } /* for i */ } #endif /* ASM_bigp_siever */ #if !(ASM_lowp_siever) lowp_siever(ibeg, iend, rsl, logp) long ibeg; /* lower bound on i for sieving */ long iend; /* upper bound on i for sieving */ long rsl; /* width of sieve interval */ uchar logp; /* logarithm of these primes */ { long loc_rsl = rsl; register long i; register urow sieve_rsl = sieve + loc_rsl, sieve1; register uchar loc_logp = logp; register roots_t *rootsi = roots + ibeg; register row pi = primbase + ibeg; for (i = ibeg - iend - 1; i < 0; i++) { register long rootdif = rootsi->r2; register long root2 = rootsi->r1 - loc_rsl + rootdif; /* larger root relative to end of sieve interval. */ /* Will always be negative initially */ register long p = *pi++; #if 0 if ( rootsi->r1 < 0 || rootdif < 0 || rootsi->r1 + rootdif >= p || p+p > loc_rsl) { printf("Bad r1 = %ld, r2 = %ld, lowp = %ld\n", rootsi->r1, rootdif, p); } #endif sieve1 = sieve_rsl - rootdif; /* Note sieve1 = sieve + (rsl-p) + (p-rootdif) >= sieve, as required by ANSI C. */ do { sieve1[root2] += loc_logp; sieve_rsl[root2] += loc_logp; root2 += p; sieve1[root2] += loc_logp; sieve_rsl[root2] += loc_logp; root2 += p; } while (root2 < -p); if (root2 < 0) { sieve1[root2] += loc_logp; sieve_rsl[root2] += loc_logp; root2 += p; } root2 -= rootdif; /* Replace root2 by smaller root */ if (root2 < 0) { /* Extra hit for smaller root */ sieve_rsl[root2] += loc_logp; root2 += rootdif; /* Restore previous value */ rootsi->r2 = p - rootdif; /* New difference */ } rootsi->r1 = root2; /* New smaller root */ rootsi++; } /* for i */ } #endif /* ASM_lowp_siever */ #if !(ASM_medp_siever) medp_siever(ibeg, iend, rsl, logp) long ibeg; /* lower bound on i for sieving */ long iend; /* upper bound on i for sieving */ long rsl; /* width of sieve interval */ uchar logp; /* logarithm of these primes */ { register long loc_rsl = rsl, i; register uchar loc_logp = logp; register urow sieve_rsl = sieve + loc_rsl; register roots_t *rootsi = roots + ibeg; register row pi = primbase + ibeg; for (i = ibeg - iend - 1; i < 0; i++) { register long root1 = rootsi->r1 - loc_rsl; register long root2 = root1 + rootsi->r2; register long p = *pi++; /* roots relative to end of sieve interval */ /* This code resembles that in bigp_siever, except we sieve at least once per root. */ #if 0 if ( rootsi->r1 < 0 || rootsi->r2 < 0 || rootsi->r1 + rootsi->r2 >= p || p+p < loc_rsl || p > loc_rsl) { printf("Bad r1 = %ld, r2 = %ld, medp = %ld\n", rootsi->r1, rootsi->r2, p); } #endif sieve_rsl[root1] += loc_logp; sieve_rsl[root2] += loc_logp; root1 += p; if (root1 < 0) { sieve_rsl[root1] += loc_logp; root1 += p; /* Advance smaller root by p */ root2 += p; if (root2 < 0) { sieve_rsl[root2] += loc_logp; /* Difference unchanged */ } else { rootsi->r2 = root1 - root2; /* New difference */ root1 = root2; /* New smaller root */ } } rootsi->r1 = root1; rootsi++; } /* for i */ } #endif /* ASM_medp_siever */ long prim_atmost(p) long p; { /* Return the largest index i such that primbase[i] <= p. Use a binary search. We treat primbase[0] as -infinity and primbase[sizefacbase+1] as +infinity. */ register long mn = 0, mx = sizefacbase; while (mx != mn) { /* primbase[mn] <= p < primbase[mx+1] */ register long mid = (mn + mx + 1) >> 1; if (primbase[mid] > p) { mx = mid - 1; } else { mn = mid; } } /* while */ return mn; } #if !(ASM_sieveroot_adjust) sieveroot_adjust(amt) long amt; /* assumed nonnegative */ { /* findnextpolynomial stored the roots of a*x^2 + b*x + c mod p into r1 and r2, in arbitrary order. That is also the convention used in checkroot. Add (amt % p) to these roots, keeping them in [0, p-1). Then store the smaller root in r1 and the difference in r2 (a convention used in bigp_siever, lowp_siever, medp_siever, and sievget). */ register row pi = &primbase[nbsmallpsp]; register roots_t *rootsi = &roots[nbsmallpsp]; register roots_t *rootsi_bound = &roots[sizefacbase]; for (; rootsi <= rootsi_bound; rootsi++) { register long p = *pi++; register long srem = amt % p; register long root1 = rootsi->r1 + srem; register long root2 = rootsi->r2 + srem; if (root1 >= p) root1 -= p; /* Ensure in [0, p-1) */ if (root2 >= p) root2 -= p; if (root1 < root2) { rootsi->r1 = root1; /* smaller root */ rootsi->r2 = root2-root1; /* difference */ } else { rootsi->r1 = root2; rootsi->r2 = root1-root2; } } /* for rootsi */ } #endif /* ASM_sieveroot_adjust */ siever(a,b,d2inv) verylong a, b, d2inv; { register long i; long value = (primbase[1] == 2 ? BYTSOF1 : 0); long ishift; long rsl_prev = 0, rsl_crossing, rsl2_crossing; /* primbase[rsl_crossing ] <= rsl <= primbase[rsl_crossing+1 ] */ /* primbase[rsl2_crossing] <= rsl/2 <= primbase[rsl2_crossing+1] */ sieveroot_adjust(sievelength); /* Add sievelength %p to roots */ /* (negative of initial ishift) */ for (ishift = -sievelength; ishift < sievelength; ishift += sievesize) { long rshift = ( ishift + sievesize < sievelength ? ishift + sievesize : sievelength); /* right end of interval, plus 1 */ long rsl = rshift - ishift; /* width of current sieve interval */ register uchar logp; #ifdef RISC for (i=0;i> 1); } for (logp = firstlog; logp <= lastlog; logp++) { register long ibeg = logpstart[logp]; long iendbig = logpstart[logp+1] - 1; long iendlow = (iendbig <= rsl2_crossing ? iendbig : rsl2_crossing); long iendmed = (iendbig <= rsl_crossing ? iendbig : rsl_crossing ); /* Sieve over primes p <= rsl/2 with this logarithm */ if (ibeg <= iendlow) { lowp_siever(ibeg, iendlow, rsl, logp); ibeg = iendlow + 1; } /* Sieve over primes p with rsl/2 < p <= rsl and this logarithm */ if (ibeg <= iendmed) { medp_siever(ibeg, iendmed, rsl, logp); ibeg = iendmed + 1; } /* Sieve over primes p > rsl with this logarithm */ if (ibeg <= iendbig) { bigp_siever(ibeg, iendbig, rsl, logp); } } /* for logp */ sievget_check(1, rshift); } /* for ishift */ } report() { /* This report is produced at intervals specified by REPORT. The output has one long line, resembling: S40; F1=137; P2(126,113); p4(105,102.25,100); B148(118,102.82,99); L13.86;12765946021730812111390339 S40 - 40 polynomials have been sieved. F=137 - 1 full found by sievget, with ival = 137 (ival is total logarithm from sieving). P2(126,113) - Two partials found, with ival = 126 and 113. p4(105,102.25,100) - Four double partials have been found, with largest ival = 105, smallest = 100, average = 102.25. B148(118,102.82,99) - 148 false reports to sievget occurred, with largest ival = 118, smallest = 99, average = 102.82. L13.86 - Average value of index on dumpfac call (divisors from factor base per relation) = 13.86. 12765946021730812111390339 - Latest startd for findnextpolynomial. */ printf("S%ld",pol); if (ful > 0) { if (ful == 1) printf("; F1=%ld",maxful); else if (ful == 2) printf("; F2(%ld,%ld)",maxful,minful); else printf("; F%ld(%ld,%.2f,%ld)",ful,maxful, (double)sumful/ful,minful); } if (par > 0) { if (par == 1) printf("; P1=%ld",maxpar); else if (par == 2) printf("; P2(%ld,%ld)",maxpar,minpar); else printf("; P%ld(%ld,%.2f,%ld)",par,maxpar, (double)sumpar/par,minpar); } if (ppr > 0) { if (ppr == 1) printf("; p1=%ld",maxppr); else if (ppr == 2) printf("; p2(%ld,%ld)",maxppr,minppr); else printf("; p%ld(%ld,%.2f,%ld)",ppr,maxppr, (double)sumppr/ppr,minppr); } if (bad > 0) { if (bad == 1) printf("; B1=%ld",maxbad); else if (bad == 2) printf("; B2(%ld,%ld)",maxbad,minbad); else printf("; B%ld(%ld,%.2f,%ld)",bad,maxbad, (double)sumbad/bad,minbad); } if ((ful+par+ppr) > 0) printf("; L%.2f",(double)numfac/(ful+par+ppr)); if (started > 0) { printf(";"); if (interrupted > 0) zwrite(altstartd); else zwrite(startd); } timestamp(); printf(" prime %ld; ecm %ld(+%ld); mis %ld; bigfac %ld; large %ld\n", prc,ecc,ecs,mic,woc,lac); fflush(stdout); fsync(fileno(stdout)); /* make sure bits get out */ } void csaver() { int status; if (outputcounter &&( how == SLAVE )) { if (tobemailed) { printf("mailing...\n"); status = system(callsendit); #ifdef VMS if ((status & 01) != 1) #else if (status) #endif printf("mail failed '%s' (status = %d)\n", callsendit, status); else unlink(nameoutput); /* delete only if mail worked */ } else printf("please mail %s to the factor account\n",nameoutput); } printf("Exit: "); report(); exit(0); } void sigtermexit() { printf("caught TERM signal, exit"); timestamp(); csaver(); exit(0); } void sigintexit() { printf("caught INT signal, exit"); timestamp(); csaver(); exit(0); } void sigstoppause() { if (!signalinterrupt) { printf("caught TSTP signal, pause"); timestamp(); } kill(getpid(), SIGSTOP); signalinterrupt = 1; } void goagain() { if (signalinterrupt) { printf("caught CONT signal, continue"); timestamp(); signalinterrupt = 0; } } void setupsignals () { signal(SIGTERM, sigtermexit); signal(SIGINT, sigintexit); signal(SIGTSTP, sigstoppause); signal(SIGCONT, goagain); } iterate() { static verylong a = 0, b = 0, d2inv = 0; sievesize = MAXSIEVE; if (sievesize > sievelength) sievesize = sievelength; if (sievesize < 400) { if ((sieve = (uchar*)malloc(luchar*(400+FIRST_GE_SENTINEL)))==NULL) { printf("unable to get memory for sieve, exit\n"); exit(0); } } else { while ((sieve=(uchar*)malloc( luchar*(sievesize+FIRST_GE_SENTINEL)))==NULL) { if ((sievesize -= MINSIEVE) < 0) { printf("unable to get memory for sieve, exit\n"); exit(0); } } } fsync(fileno(stdout)); /* make sure bits get out */ for (pol=1; pol<=MAXPOL; pol++) { findnextpolynomial(&a,&b,&d2inv); if (how == SLAVE) if (((pol%CHECK) == 0) /* || (pol == 1 ) */ ) checkroot(a,b); siever(a,b,d2inv); if ((pol%REPORT) == 0) report(); } csaver(); } #ifdef VMS unlink(filename) char* filename; { int status; char unlink_command[100]; sprintf(unlink_command, "DELETE %s.;", filename); status = system(unlink_command); if ((status & 01) != 1) { printf("VMS delete failed '%s' (status = %d)\n", unlink_command, status); fflush(stdout); } } #endif /*VMS*/ main (argc, argv) char *argv[]; int argc; { #ifdef SYSV setpgrp2(0,0); #endif if (argc >= 4) { restartfn = argv[1]; if (strcmp(argv[2], "-") && !freopen(argv[2], "r", stdin)) { fprintf(stderr, "Can't open input file '%s'\n", argv[2]); exit(1); } if (!freopen(argv[3], "a", stdout)) { fprintf(stderr, "Can't open output file '%s'\n", argv[3]); exit(1); } } if (argc == 1) restartfn = default_restart; else { restartfn = argv[1]; printf("Restart file is %s\n", restartfn); } if (argc >= 5) file_prefix = argv[4]; if (argc >= 6) mail_cmd = argv[5]; fflush(stdout); setupsignals(); timestart = getutime(); scanf("%ld",&how); setup(); /* Find where the logarithm changes. */ { register uchar logp = firstlog; logpstart[logp] = nbsmallpsp; while (1) { long logbr = nextlog[logp-firstlog]; if (logbr > primbase[sizefacbase]) break; logpstart[++logp] = prim_atmost(logbr-1) + 1; } lastlog = logp; logpstart[lastlog+1] = sizefacbase + 1; } iterate(); }