# # Quadratic sieve with extra features # Needs li.d4f # Extended factor base # Sieve forms (ax+b)^2-M # "mpqs.doc" >> EOF The idea is to sieve values y=(ax+b)^2-M and pick out those with small factors. Sieving takes much time with large factor bases. The file system interface works as DATA<-#read N, START, SIZE, RP which asks for SIZE bytes from START+SIZE*RP. The same could be done for the sieve function: VALUES<-FBTAB SIEVE_READ A," ", B, " ",(zBS), " ", z RP The big problem is the number of different interfaces involved. The paramaters for SIEVE_READ are a numeric table, and the start, block-size, index combination. The values A are assigned by random selection within a certain range (32-80 for a start) and corresponding values B are computed to make sure (AX+B)^2-M is divisible by A. Sieving takes place over two sets of primes. Sieve for prime powers amongst small primes, and sieve for single divisors for most of the larger primes. This can be speeded up by making sure it is easy to compute the two offsets corresponding to solutions (AX+B)^2=M mod q. If r1 and r2 are roots of X^2=M mod q, then the required values are A'(r1-B) mod q where A' is the inverse of A mod q (AA'=1 mod q). EOF Z<-QUAD_TESTM STR;B;M;X;P |Test whether number is prime, or has easy factors B<-1200 M<-"-b B<-#fi SHIFT" GETOPTS STR IF (rM)<9; Z<-z FACTORIZE32 #fi M; ->0 P<-2, SIEVE (B+1)%2 IF 0< r X<-(0=M #mod P)/P; Z<-X," ", M% X<-PRODUCT X; ->0 IF "0" #equiv M-X*X<-Y_SQRT M; Z<-X," ",X ; ->0 IF M #equiv X<-SCAN_FTEST M; Z<-X,"p"; ->0 Z<-"" | Nothing known Z<-LST SELECT_MULTIPLIER STR;B;FB;M;QR;R |Select multiplier k to sieve values X^2-k*M LST<-1 3 5 7 11 13 17 19 23 29 31 37 41 43 47 #ifndef "LST" B<-1200; QV<-0 M<-"-b B<-#fi SHIFT|-v QV<-1" GETOPTS STR B<-4200 f B Z<-iI<-0;FB<-SIEVE (1+B)%2 |Create a (rFB) x rLST table. WHILE I 1 |Supercedes brute force Z<-(1=(A*Z) #mod M)/Z<-iM B<-M & S<-0 & Z<-1 & U<-0 WHILE 0< R<-B-X*Q<-B%X B<-X & X<-R & S<-1-S & T<-Z & Z<-U+Z*Q & U<-T WEND IF ~1= X; Z<-X & S<-0 IF S;Z<-M-Z Z<-A MOD_SQRT M;X |Square root of A mod M |Brute force for integers in a first level factor base Z<-(X*X<-iM) #mod M Z<-(Z=A)/X Z<-Q_PRINT X;SZ SZ<-#screen; Z<-(SZ[1]-8)tX Z<-QS_PRINT_MATRIX STR;E;EFB;FW;I;KA;MK;NN;NR;OFILE;Q;T;UFB;V;X;XSQ |Print rows of matrix in full. |Import FB, M, X,TAB, BP.TAB KA<-1; OFILE<-Z<-"" STR<-"-o OFILE<-SHIFT|-k KA<-#fi SHIFT" GETOPTS "" #ifndef "STR" #print "M=", M #print "k=", zKA #print "fb= ", zFB #print "start columns ", zrFB FW<--1+rMK<-M*zKA; I<-0; EFB<-""; NN<-#sstomat #use "BP.TAB" #print "bp.tab(######)" z 1tr:NN #print "x.tab (######)" z 1tr: #sstomat #use "X.TAB" WHILE I<1tr:NN Q<-#deb NN[I++];T<-BP.TAB[Q]; IF 1<+/T=" "; EFB<-EFB,"/",Q WEND EFB<-#fi , " ", #sstomat EFB EFB<-EFB[#sort EFB] #print "xfb= ", z EFB |Create rows Z<-#use "X.TAB" | string of numbers I<-0 WHILE I0)/iNC; Z<-Z[I]; ->0 | A row sum is zero C<-,X[:J<-1tJ]; K<-(C[K] * I #ne K)/K<-iNR; DX<-(rK)/:R; DZ<-(rK)/:Z[I] X[K]<-X[K] #ne DX; Z[K]<-Z[K] #ne DZ I++ WEND Z<-RSQ_MATRIX_PHASE STR;DM;E;FB;FW;I;KA;M;QR;QV;TAB;V;WP;XS;X |Analyse a file constructed in phase one Z<-""; KA<-"1"; QV<-0; QR<-0; WP<-5 | Default transcript SRC<-"-w WP<-#fi SHIFT|-v QV<-1|-r QR<-1" GETOPTS STR<-"" #ifndef "STR" IF 00 NR<-#fi (\$T[Z]) YDX 1 XS<-(NR,FW) t TAB<-\$T[Z+1+iNR]; MA<-(NR,-NC) t TAB E<-~f/MA="."; NC<-rFB<-E/FB;MA<-E/MA | Crunch blank columns IF 00 I<-I/irI X<-M C_PRODUCT XS[DM+I] Y<-PRODUCT ((+/MA[DM+I])%2)/FB Y<-Y #mod M #print V<-"(x,y) ",X," ",Y; IF QV;V I<-M EA X+Y |Here we have X^2=Y^2 #mod M BREAKIF ~ (I #equiv "1") c I #equiv M #print V<-"trivial factor"; IF QV;V DM++; WEND #print Z<-I," ",M%I; IF QV; Z Z<-M C_PRODUCT TAB;I;NR |Product of long table mod M Z<-"1"; I<-0; NR<-1tr:TAB WHILE IAA X<-iPN<-*/NP r P | Deal with powers of primes MP<-#fi M #mod z PN Y<-((X*X)-MP) #mod PN Y<-+/y 0=Y #outer #mod *\NPrP AA: NZ<-+/Y Z<-Z,: P, NP, (*/NPrP), ("N" #av LP), NZ, MP, DM WEND Z<-BP_ADD STR;N;NC;X;T |Large prime stuff. Update number of relations |Global BP.TAB, XP.TAB Z<-0 0 IF 0 < XP.TAB[X<-STR YDX 1]; ->0 |Reject duplicates XP.TAB[X]<-1 NC<-+/" "=T<-BP.TAB[N<-STR YDX 0] BP.TAB[N]<-T," ",X IF NC=1; Z<-2 1 IF NC>1; Z<-1 0 Z<-TAB MK_FLIST STR;A;AMAX;BS;DX;FB;FC;J;M;MA;X;XC;XI;XS;Y |Make list of pairs a, b to sieve (ax+b)^2-M BS<-64000; NF<-4 FC<-SHIFT; BS<-#fi SHIFT; M<-SHIFT AMAX<-64 80 80 128 256 512[+/(rM)>10 20 24 30 36] "AMAX=", zAMAX XS<-"1" + Y_SQRT M IF FC="1"; Z<-"1 ",XS," ",zBS; ->0 IF FC e "23456789"; NF<-#fi FC XI<-AMAXr0 |A Less than AMAX FB<-(FB<80)/FB<-,TAB[:0] | Tiny factor base XI[FB]<-1 |Any useful prime XI[,FB #outer * FB]<-1 |Product of pairs XI[FB*FB]<-0 |Exclude squares Z<-(Z>40)/Z<-XI/irXI FB<-1, Z[(NF-1)?rZ] Z<-"" WHILE 0 < r FB; A<-1tFB;FB<-1dFB |construct b so that b>x; (ax+b)^2-M = 0 #mod a IF A=1; Z<-Z,"/1 ", XS, " ",zBS; REPEATIF 1 MA<-M #mod A DX<-XS #mod A Y<-((X*X<-iA)-MA) #mod A J<-1 t (Y=0)/iA XC<-XS + z (J-DX) #mod A Z<-Z,"/",(zA)," ", XC, " ", zBS WEND Z<-#sstomat Z Z<-A V_MODINV LST;I;J;K;N;NC |Calulate inverse values A' so that AA'=1 (mod LST[i]) N<-(rA)*NC<-rLST; Z<-Nr0 K<-I<-J<-0 WHILE N>K Z[K++]<-A[I] MOD_INV LST[J++] IF J=NC; J<-0; I++ WEND Z<-((rA), NC) r Z Z<-B V_MOD LST;I;N;NR |Compute integer matrix #fi B[i] #mod FB[j] I<-0; NR<-1tr:B; Z<-(NR, rLST)r0 WHILE I< NR; Z[I:]<- (#deb B[I:]) #mod LST; I++; WEND R<-LST FB_DIV M;F;H;I;J;IJ;MAX |Factorise M with a given factor base H<-1; I<-0; R<-(rLST)r0; MAX<-MAX*MAX<-c/LST IJ<-(0 = M #mod LST)/irLST |Primes dividing M WHILE I0 Z<-0 r 5 #sed":d" BP<-"N" #av #exp 1.5 * #ln B | Test size for large primes Z<-2 16013 #use"BP.TAB";Z<-10001 #use "X.TAB"; Z<-10001 #use "XP.TAB" KA<-SELECT_MULTIPLIER "-v -b ",(zB)," ", M MK<-M*z KA SIEVE_T<-"N" #av (0.5 * Y_LN MK)% #ln 2 TAB<-QUAD_PNPP "-v -b ", (zB)," ", MK FB<-,TAB[:0] SZ<-0, 16+rFB FLIST<-TAB MK_FLIST FC," ", (zBS), " ", MK #print "M=",M #print "digits= ###" zrM #print "k=", zKA #print "bs=########" zBS #print "b=###### [####] " z B, rFB #print "P=", zBP AI.LST<- #fi ," ",0 Y_COLUMN FLIST AI.TAB<- AI.LST V_MODINV FB BI.TAB<- (1 Y_COLUMN FLIST) V_MOD FB #print "a=", z AI.LST IB<-JP<-0 WHILE SZ[0]SIEVE_T-3)/irTHETA Q_PRINT "selection completed (####)" z rIX J<-0 AA: | treat selection from block BREAKIF SZ[0]>=SZ[1] REPEATIF J=rIX X<-XC + z AC*IX[J++]; XSQ<-(X*X)-MK IF KA>1; IF "0" #equiv XSQ #mod z KA; ->AA Z<-1 t FB FB_DIV XSQ IF QV;0 r (70 t ("####/#### [####] ######## " zSZ, J, Z)) #wput 1 0 1 70 IF (Z<0) c Z>BP; ->AA |Check whether X value is already used IF (1 = Z) f X.TAB[X]>0; ->AA IF 1=Z; X.TAB[X]++; SZ[0]++; Q_PRINT("####/#### X="z SZ),X; ->AA SZ<-SZ+K<-BP_ADD (zZ)," ",X; IF 0AA Y<-+/y 0=Y #outer #mod *\NPrP | Prime powers AA: IF QV;(70 t "sieving ", 6 z P, NP, PN, MP, DX) #wput 1 0 1 70 Z<-Z + BS r NB * Y |Accumulate by dense addition I++; WEND |Sieve tail of list WHILE I