#
#       Dixon's random squares factoring
#       Important in theory.
#       find values x, y with x^2=y^2 #mod N
#       The real problem is finding so called smooth quadratic
#       residues (mod N). Continued fraction and variants of the
#       quadratic sieve are common.
#       Matrix algebra is a bottleneck in large problems because of
#       the vast amount of memory required.
#

Z<-RANDOM_SQUARES STR;B;E;FB;FW;I;J;M;NR;NC;QV;SZ;X;Y;XI;XSQ
|Random squares method
B<-100; QV<-0; SZ<-#screen
FW<- -1+rM<-"-b B<-#fi SHIFT|-v QV<-1" GETOPTS STR
Z<-0 r 5 #sed":d"; I<-J<-0; MA<-i0; XLIST<-""
#print "M= ",M
#print "digits= ###  digit sum= ####" z (rM), +/(#av M)-#av "0"
NR<-16+NC<-rFB<-2, SIEVE (B+1)%2
#print "b= ######  [####] " z B, rFB
#print "ufb=", z FB
#print "#","matrix #### ####" z NR, NC
WHILE J<NR
XSQ<-(X*X<-YRAND M) #mod M; I++
Z<-FB FB_DIV XSQ
IF QV;(60 t "[######] ######" zI, 1tZ) #wput 1 12 1 60
REPEATIF 1 #ne 1tZ
XI<-1dZ
E<-XI%2; XI<-XI-2*E
IF 0<c/E; X<-(X*((PRODUCT E/FB) YMOD_INV M)) #mod M
XSQ<-(X*X) #mod M
#print V<-(FWtX),(FWtXSQ)," ", ".o"[XI]
(SZ[1]-6) t ("####/#### "z J,NR),V
J++
WEND
#print "#end"
IF QV; #sed":u rsq.log"
"Matrix phase"
Z<-RSQ_MATRIX_PHASE "-w 5 -v"
IF QV; #sed":u rsq.log"

Z<-RSQ_MATRIX_PHASE STR;DM;E;FB;FW;I;M;QR;QV;TAB;V;WP;XS;X
|Analyse a file constructed in phase one
Z<-""; 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
FW<-1+rM
"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; #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

#
# Relations amongst rows
# Find a set of rows which sum to zero (mod 2)
#

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<-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<-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<-CFRAC STR;A;AC;E;FB;FW;I;J;M;NR;NC;P;Q;QV;SQRTM;SZ;T;V;X;XI;XSQ
|Continued fraction method to factorize M
QV<-0; SZ<- #screen; B<-1200
FW<- -1+rM<-"-b B<-#fi SHIFT|-v QV<-1" GETOPTS STR
Z<-10001 #use "X.TAB"
Z<-"" & SQRTM<-Y_SQRT M
IF "0" #equiv M-SQRTM*SQRTM; Z<-SQRTM & ->0
Z<-0 r 5 #sed":d"; I<-J<-0; MA<-i0; XLIST<-""
NR<-16+NC<-rFB<-2, SIEVE (B+1)%2
#print "M= ",M
#print "digits= ###  digit sum= ####" z (rM), +/(#av M)-#av "0"
#print "b= ######  [####] " z B, NC
#print "ufb=", zFB
#print "#","matrix #### ####" z NR, NC
P<-"0" & Q<-"1" & I<-0 & AC<-"0 1"
WHILE I<NR
A<-(P+SQRTM)%Q
T<-(A*Q)-P; P<-T; Q<-(M-T*T)%Q
AC<-(AC YDX 1)," ",((AC YDX 0)+(A*AC YDX 1)) #mod M
J++
REPEATIF (J=1) c ~ J #and 1
X<-AC YDX 0
Z<-FB FB_DIV XSQ<-(X*X) #mod M
IF QV;(60 t "[######] ######" zJ, 1tZ) #wput 1 12 1 60
REPEATIF 1 #ne 1tZ
REPEATIF 0<X.TAB[X]; X.TAB[X]<-1
XI<-1dZ; E<-XI%2; XI<-XI-2*E
IF 0<c/E; X<-(X*((PRODUCT E/FB) YMOD_INV M)) #mod M
REPEATIF X #equiv "1"   | Should try immediate factorisation
REPEATIF 0<X.TAB[X]
XSQ<-(X*X) #mod M
IF (XSQ #equiv "1") f ~M #equiv E<-X+"1";Z<-(M%Z)," ", Z<-M EA E;->0
#print V<-(FWtX),(FWtXSQ)," ", ".o"[XI]; I++
(SZ[1]-8) t ("####/#### "z I,NR ), V
WEND
#print "#end"
IF QV; #sed":u cfrac.log"
"Matrix phase"
Z<-RSQ_MATRIX_PHASE"-w 5 -v"
IF QV; #sed":u cfrac.log"

Z<-QUAD_ZERO STR;A;E;FB;FW;I;J;M;NR;NC;QV;SQRTM;SZ;T;V;X;XI;XSQ
|Use values X*X-M. Try to factorize at once, without sieve
B<-100; QV<-0; SZ<-#screen
FW<- -1+rM<-"-b B<-#fi SHIFT|-v QV<-1" GETOPTS STR
Z<-0 r 5 #sed":d"; I<-J<-0; MA<-i0; XLIST<-""
SQRTM<-Y_SQRT M
#print "M= ",M
#print "digits= ###  digit sum= ####" z (rM), +/(#av M)-#av "0"
NR<-16+NC<-rFB<-2, SIEVE (B+1)%2
#print "b= ######  [####] " z B, rFB
#print "ufb=", z FB
#print "#","matrix #### ####" z NR, NC
I<-J<-0
WHILE I<NR
X<-SQRTM + zJ<-J+1
XSQ<-(X*X)-M
Z<-FB FB_DIV XSQ
IF QV;(60 t "[######] ######" zJ, 1tZ) #wput 1 12 1 60
REPEATIF 1 #ne 1tZ
XI<-1dZ
E<-XI%2; XI<-XI-2*E
IF 0<c/E; X<-(X*((PRODUCT E/FB) YMOD_INV M)) #mod M
XSQ<-(X*X) #mod M
#print V<-(FWtX),(FWtXSQ)," ", ".o"[XI]
(SZ[1]-6) t ("####/#### "z I,NR),V
I++
WEND
IF QV; #sed":u qs.log"
"Matrix phase"
Z<-RSQ_MATRIX_PHASE "-w 5 -v -r"
IF QV; #sed":u qs.log"