#
#       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<rFB; R<- #fi M #mod z P<-FB[I++]
QR<-Pr0; QR[(X*X<-iP) #mod P]<-1        | Mark squares
IF QV; V_SLIDER 0 8, I,rFB
Z<-Z,:QR[(R*LST) #mod P]
WEND
Z<-((#ln FB)%FB-1) #inner Z
Z<-Z-0.5* #ln LST
Z<-1 t (Z=c/Z)/LST

Z<-V_SLIDER CX
|Progress bar. position, fraction<1
Z<-"X."[((50*CX[2])%CX[3])<=i50]
Z<-0 r Z #wput CX[0 1],1, rZ

Z<-X MOD_INV M;B;R;Q;S;T;U
|Find Z with X*Z=1 (mod M) and zero if gcd (X,M) > 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 I<rEFB
Z<-Z,"/ /\r/*" #rxsubs T<-BP.TAB[zEFB[I++]]
WEND
#print "ufb=", z UFB<-FB,EFB
I<-0; NR<-1tr:Z<-#sstomat Z
#print "#", "matrix #### ####" z NR, rUFB
WHILE I<NR; X<-#deb Z[I++]; XSQ<-(X*X)-MK
REPEATIF 1 #ne 1tT<-UFB FB_DIV XSQ
T<-1dT;E<-T%2;T<-T-2*E
IF 0<c/E; X<-(X*((PRODUCT E/UFB) YMOD_INV M)) #mod M; XSQ<-(X*X) #mod M
#print V<-(FWtX),(FWtXSQ)," ", ".o"[T]
Q_PRINT V
WEND
#print "#end"
IF 0 < r Z<-OFILE; 0 r #sed":u ", OFILE

Z<-QS Z2_SCAN_ROWS X;C;DX;DZ;I;J;K;NR;NC;R
|Brute force attempt to return a set of rows which sum to zero
|It is slow for large sets
QS<-0 #ifndef "QS"      | Slider
Z<-(iNR<-1tr:X) #outer = iNC<-1dr:X; I<-0
WHILE I<NR; R<-X[I]
IF QS; V_SLIDER 1 8, I, NR
IF 0=rJ<-(R>0)/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 0<rSRC; Z<-WP #sed":r ", SRC
IF 0 = r Z<- #nfind "M="; ->0
M<-(2d$T[1tZ]) YDX 0            |Number to factor
"M=", M
IF 0 = r Z<- #nfind "ufb="; ->0
NC<-rFB<- #fi 4d $T[1tZ]        |Factor base
IF 0<r Z<-#nfind"k="; KA<-#deb 2d $T[1tZ]
FW<-1+r M*KA
"FB =", 50t zFB
IF 0 = r Z<-#nfind "#matrix"; ->0
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 0<Z<-+/E=0; #print "crunch #### columns" z Z
#print V<-"matrix size = ##### ##### " z NR, NC; IF QV; V
IF QR; XS<-XS[I<-NR?NR]; MA<-MA[I]      | Permute rows
Z<-M;DM<-0; MA<-MA #ne "."
WHILE DM+NC<NR
#print V<-"z2_scan_rows (####)" z DM; IF QV; V
IF ~1=1tr:I<-QV Z2_SCAN_ROWS MA[DM+iNC]; ->0
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 I<NR; Z<-(Z * #deb TAB[I++ :]) #mod M; WEND

Z<-QUAD_PNPP STR;B;BS;DM;FB;I;LN2FB;LP;M;MAXPP;NP;MP;P;PN;QV;X;Y
|Prepare factor base with some prime powers
BS<-64000; B<-100; QV<-0
M<-"-b B<-#fi SHIFT|-bs BS<-#fi SHIFT|-v QV<-1" GETOPTS STR
FB<-2, SIEVE (1+B)%2
LN2FB<- (#ln FB)% #ln 2
MAXPP<- (#ln B c #sqrt BS) % #ln 2
Z<-iI<-0
WHILE I<rFB; LP<-LN2FB[I]; MP<- #fi M #mod z P<-FB[I++]
IF QV; V_SLIDER 1 8, I, rFB
REPEATIF MP=0
REPEATIF 0=+/Y<-0=((X*X<-iP)-MP) #mod P
DM<-2 r Y/iP
NP<- "N" #av MAXPP % LP
IF NP=1; ->AA
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 I<rIJ
F<-LST[J<-IJ[I]]
IF 0 #ne M #mod F;I++;REPEATIF 1
M<-M% zF; R[J]++;
BREAKIF "1" #equiv M
WEND
IF ~"1" #equiv M; H<-_1; IF "-"=1t M-zMAX; H<-#fi M
R<-H,R

Z<-QUAD_MPQS STR;AC;B;BP;BS;FC;FB;FLIST;IB;IX;J;JP;K;KA;M;MK;QV;SQRTM;SZ;TAB;THETA;X;XP;XSQ
|Quadratic sieve on (AX+B)^2*X-K*M
|multiplier K
|Large primes added as required to factor base
|Global BP.TAB, X.TAB, XP.TAB
FC<-"4"; QV<-0; BS<-64000; B<-1600
M<-"-fc FC<-SHIFT|-bs BS<-#fi SHIFT|-b B<-#fi SHIFT|-v QV<-1" GETOPTS STR
IF 0<rZ<-QUAD_TESTM "-b ",(zB), " ", M; ->0
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]<SZ[1]
IF JP=1tr:FLIST; IB++; JP<-0;
XP<-FLIST[JP++ :]; AC<-#fi XP YDX 0; XC<-(XP YDX 1)+z IB*AC*BS
#print Z<-"#", "A=#### block ####/#### (##%)" z AC, SZ,(100*SZ[0])%SZ[1] & Z
THETA<-TAB READ_SIEVE XP," ",zIB
Q_PRINT "select phase:"
IX<-(THETA>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 0<c/K; Q_PRINT ("####/#### X=" z SZ),X," Q=",zZ
->AA
WEND
Z<-QS_PRINT_MATRIX "-k ", z KA
Z<-RSQ_MATRIX_PHASE "-v -r"
0 r #sed":u mpqs.log"

Z<-TAB READ_SIEVE STR;A;AJ;DI;DX;I;IA;IB;J;N;NB;NC;MP;NP;P;PN;PT;Q;RP;T;X;Y;XZ
|Sieve a block (ax+b)^2-M using a table
|Import AI.LST, AI.TAB  -- inverse of A table
|Import BI.TAB          -- offsets XZ mod P
I<-0;N<-1tr:TAB
A<-#fi SHIFT; XZ<-SHIFT; BS<-#fi SHIFT; RP<-#fi SHIFT
IA<-AI.TAB[AI.LST i A :]
IB<-BI.TAB[AI.LST i A :]
PT<-"N" #av #sqrt BS    | Large primes
Z<-BSr0
XZ<-XZ+ z A*BS*RP
WHILE I<N; T<-TAB[I]
BREAKIF PT< P<-T[0]
NP<-T[1];PN<-T[2];NB<-T[3];MP<-T[5]
DX<-XZ #mod PN
X<-(DX+A*iPN) #mod PN; Y<-((X*X)-MP) #mod PN
IF NP=1; Y<-Y=0; ->AA
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<N
J<-I+iDI<-50 f N-I
T<-TAB[J:]
P<-,T[:0]
NC<-(BS+P[0]-1)%P[0]
NB<-(T[0:])[3]
IF QV;(70 t "sieving ", 6 z P[0], NC, NB) #wput 1 0 1 70
Q<- yT[:6 7]
AJ<-IA[J]
DX<-IB[J]
X<-(BS*RP) #mod P; DX<-(DX+A*X) #mod P
PN<-P #outer * iNC
XI<-BSr0
X<-(AJ*Q[0]-DX)#mod P; Y<-,PN+NC/yX; XI[Y]++
X<-(AJ*Q[1]-DX)#mod P; Y<-,PN+NC/yX; XI[Y]++
Z<-Z+NB*XI
I<-I+DI
WEND