Z<-LI.D4F
|Long integer routines
Z<-ADD_ALIAS LI.D4S

"LI.D4S"
edli|WW.ED"li.d4f ",ARGS
ldli|#copy"li.d4f" & LI.D4F
li_info|3 #sed"X:d;$ li.inf" & Z<-WW.FSCR 3
factorize|G<-P_FACTORS H<-ARGS & G
fretry|Z<-TEST_FILE ARGS
rsanum|Z<-#sstomat _CR,LOAD"rsac.num" & RSA_NUM<-Z[2+?7]YDX _1
ligo|SEL.FILE "*.num"
et|T0<-24 60 60 b 3 d #ts & 0 r DO ARGS & T1<-24 60 60 b 3 d #ts & "et=",8 z T1-T0
indent|IN<-#fi ARGS; LC<-#sed"." & $T[LC+iIN[0]]<-(INr" "),$T[LC+iIN[0]]
smtest|Z<-SM_TEST ARGS

"li.inf">>EOF
    Some routines for 32-bit integers for complete factorisation.

    See Pollard's rho-method of finding factors. This was discovered
around 1976, and is very simple to program. The method has a striking
similarity to the generation of Mandelbrot set images. The same
dynamic system z<-z^2+c is used in both algorithms.
EOF

"pema.lnk">>EOF
<A href="http://www.faqs.org/rfcs/rfc2313.html">rfc2313 RSA Laboratories</A><BR>
<A href="http://www.rsasecurity.com/rsalabs/challenges/factoring/challengenumbers.txt">The RSA challenge numbers</A>
<A href="jobzap.htm">jobzap.htm</A><BR>
<A href="factor.txt">factor.txt</A><BR>
EOF

Z<-RSA_SELECT STR
|get rsa number by name
Z<-0 #sed":r rsac.num"
->(0=r Z<-#nfind STR)/0
Z<-$T[1tZ] YDX _1

Z<-FACTORS N;X
Z<-(Z>0)/Z<-(0=N #mod X)*X<-1+i "N" #av 1+#sqrt N

Z<-PRIMES N
|Generate prime numbers less than N
Z<-(2=+/y0=N #outer #mod 1+N)/N<-iN

Z<-PRIMES_IN LST
|Extract primes in list
Z<-(2=+/y0=LST #outer #mod 1+ic/LST)/LST

R<-SIEVE N;J;K;Z
|SIEVE OF ERASTOTHENES FOR ODD PRIMES
|Get odd primes up to 2*N
R<-(1,N-1)/"10" & Z<-iJ<-0
WHILE (J*J)<2*N;Z<-Z,K<-Ri"0";R[K+J*i(N+J-K+1)%J<-1+2*K]<-"1";WEND
R[Z]<-"0" & R<-1+2*(R="0")/iN

R<-B_SIEVE BLOC;C;I;K;MAX;MIN;N;P
|STRIKE OUT ARITHMETIC PROGRESSIONS FROM A BLOCK.
|This is useful even for powerful computers
|Example: P B_SIEVE START,SIZE
MIN<-BLOC[0]
MAX<-MIN + N<- BLOC[1]
IF 0 = #nc "P"; P<-2,SIEVE 1+ ("N" #av #sqrt MAX+1)%2
I<-0; R<-Nr"0"
WHILE I<rP
BREAKIF MAX < K*K<-P[I++]
IF 0< C<-MIN #mod K; C<-K-C
R[C+K*i(N+K-C+1)%K]<-"1"
WEND
R<-MIN+(R="0")/iN

Z<-PRODUCT X;N
|Multiple precision product
IF 0=N<-rX; X<-"1"
IF N<=1; Z<-z X & ->0
N<-N%2 & Z<-(PRODUCT NtX) * PRODUCT NdX

Z<-FACTORIZE32 N;FB;T;X
|Factorize ordinary integer (up to about 10^10)
FB<-(Z>0)/Z<-(0=N #mod X)*X<-1+i "N" #av 1+#sqrt N
FB<-(2=+/y0=FB #outer #mod FB)/FB
Z<-i0
WHILE 0<rFB
IF 0=N #mod T<-1tFB; Z<-Z,T & N<-N%T & REPEATIF 1
FB<-1dFB
WEND
Z<-Z,(N>1)/N

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<-ODD N
Z<-1 #and #av _1tN

Z<-EVEN N
Z<-~1 #and #av _1tN

Z<-A EA B;C;D;R
|Euclidean algorithm with A>B
R<-A #mod B & D<-B
WHILE ~R #equiv "0"; C<-D & D<-R & R<-C #mod D & WEND
Z<-D

Z<-G FTEST M;E;A
|Evaluate G^(M-1) mod M
|If M is prime, this should be 1
G<-"2" #ifndef "G"
Z<-"1"  & E<-M-"1" & A<-G
WHILE ~E #equiv "0"
IF ODD E; Z<-(Z*A) #mod M
E<-E%"2" & A<-(A*A) #mod M
WEND

Z<-SCAN_FTEST STR;M;QP;PP
|Scan forward from N for primes
QV<-0
M<-"-v QV<-1" GETOPTS STR
IF (rM)<6; M<-#fi M; Z<-z 1t(M <= PP)/PP<-SIEVE (M+200)%2; ->0
PP<-2, SIEVE 1000
QP<-0 & M<-M-"12"[1 #and #av _1tM]
WHILE QP=0; M<-M+"2"
IF QV;(_72tM) #wput 0 4 1 72
REPEATIF c/0=M #mod PP
REPEATIF ~ "1" #equiv Z<-FTEST M
|test for powers of 3 and 5
IF QV; "N=",M
Z<-"3" FTEST M;   IF QV; "3^(N-1)  =", Z
REPEATIF ~"1" #equiv Z
Z<-"19" FTEST M;  IF QV; "19^(N-1) =", Z
REPEATIF ~"1" #equiv Z
Z<-"163" FTEST M; IF QV; "163^(N-1)=", Z
REPEATIF ~"1" #equiv Z
QP<-1
WEND
Z<-M

Z<-RLONG N
N<-20 #ifndef"N"
Z<-"0123456789"[(Nt1)+?(1,0 c N-1)/9 10]

Z<-Y_PNPP N;XI;J;K;LN;P;Q
|Set of primes with repetitions of prime powers less than N
J<-0 & LN<-#ln N & XI<-(N+1)r_1 & P<-2,SIEVE (N+1)%2
XI[P]<-irP
WHILE N>Q*Q<-P[J++]
K<-"N" #av (LN)% #ln Q
XI[*\KrQ]<-XI[Q]
WEND
Z<-P[(XI>_1)/XI]

Z<-X YMOD_INV M;B;R;S;T;U
|Find Z with X*Z=1 (mod M) and zero if gcd (X,M) > 1
IF "-" = 1tX; X<-M-(1dX) #mod M
B<-M & S<-0 & Z<-"1" & U<-"0"
WHILE ~ "0" #equiv R<-B #mod X
B<-X & X<-R & S<-1-S & T<-Z & Z<-U+Z*$QUOTIENT & U<-T
WEND
IF ~ "1" #equiv X; Z<-"z",X & S<-0
IF S;Z<-M-Z

Z<-XP YMOD_EXP M;X;E
|Z<-X^P mod M
Z<-"1" & X<-XP YDX 0 & E<-XP YDX 1
WHILE ~E #equiv "0"
IF 1 #and #av _1tE; Z<-(Z*X) #mod M
E<-E%"2" & X<-(X*X) #mod M
WEND

Z<-Y_MIN STR;N
|Minimum of long numbers
Z<-SHIFT
WHILE 0<rN<-SHIFT
IF "-" = 1 t N-Z; Z<-N
WEND

Z<-Y_MAX STR;N
|Maximum of long numbers
Z<-SHIFT
WHILE 0<rN<-SHIFT
IF "-" = 1 t Z-N; Z<-N
WEND

Z<-Y_SUM STR;N
|Sum of long numbers
Z<-"0"
WHILE 0<rN<-SHIFT;Z<-Z+N;WEND

Z<-Y_PRODUCT STR;N
|Product of long numbers
Z<-"1"
WHILE 0<rN<-SHIFT;Z<-Z*N;WEND

Z<-Y_XSQ STR;X
|Square elements in list
Z<-""
WHILE 0<rX<-SHIFT; Z<-Z," ",X*X; WEND
Z<-#deb Z

Z<-Y_SQRT NN;X
|Integer square root of long
IF 1=rNN;Z<-"0111222223"[#fi NN] & ->0
X<-"0" & Z<-RLONG 1+(rNN)%2
WHILE ~X #equiv Z
X<-Z
Z<-(X+NN%X)%"2"
IF "1" #equiv Z-X; Z<-X & ->0
WEND

Z<-X YMOD M
|Modulus function including negative X
IF "-" #ne 1tX; Z<-X #mod M & ->0
Z<-M-("0"-X) #mod M

Z<-YRAND X;N
|Long random 0<=Z<X
IF 6>N<-rX; Z<-z ? #fi X & ->0
Z<-(z ? #fi 5tX),"0123456789"[?(N-5)r10]

Z<-VAR YDX J;K;U;V;M
|Get components J from string of long ints.
Z<-"" & V<-+\U<-" "=VAR<-#deb VAR & M<-1++/U
WHILE 0<rJ;K<-1tJ & J<-1dJ
Z<-Z," ",(V=K #mod M)/VAR
WEND
Z<-#deb Z

Z<-Y_LN X;N
|Natural log of long int
IF 8>N<-rX; Z<-#ln #fi X & ->0
Z<-(#ln #fi 8tX) + (N-8) * #ln 10

Z<-N Y_ROOT X;CNT;T;P
|Integer N th root of long X. Not good for N>20 or so
T<-"0" & CNT<-300 & Z<-RLONG 1+(rX)%N
WHILE ~T #equiv Z
T<-Z
P<-(N-1) Y_POWER T
Z<-(X+(zN-1)*T*P)%(zN)*P
BREAKIF 0>=CNT--
WEND

Z<-N Y_POWER X
Z<-"1"
WHILE 0<N--;Z<-Z*X;WEND

Z<-PP M_START M;T
|Divide out small facters up to a limit
Z<-""
PP<-(2, SIEVE 1000) #ifndef "PP"
T<-PP FB_DIV M
Z<-(1dT)/PP
IF T[0]=1;Z<-zZ,1;->0
Z<-(zZ)," ",M % PRODUCT Z

Z<-P_FACTOR STR;CNT;A;MAX;N;P;Q;QV;SEED;T;X;Y
|Pollard's rho-method.
|See Reisel [2] pp 180-181
MAX<-18000 & A<-SEED<-"" & QV<-0 & P<-20 & CNT<-0
N<-"-cnt MAX<-#fi SHIFT|-a A<-SHIFT|-seed SEED<-SHIFT|-v QV<-1" GETOPTS STR
IF 0=rA; A<-YRAND N
IF 0=rSEED; SEED<-YRAND N
Q<-"1" & X<-Y<-SEED & Z<-N
WHILE MAX>CNT++
X<-(A+X*X) #mod N;Y<-(A+Y*Y) #mod N;Y<-(A+Y*Y) #mod N;Q<-(Q*Y-X) #mod N
REPEATIF 0 #ne CNT #mod P
IF QV;(8zCNT),"  Q=", 66 t Q
IF "1" #equiv T<-Q EA N; REPEATIF 0< P<-240 f P+20
Z<-"z", T
IF QV; "cnt =", z CNT
BREAKIF 1
WEND

Z<-P_FACTORS STR;CNT;MAX;N;Q;QV
Z<-"" &  QV<-0 & MAX<-5 & CNT<-"18000"
N<-"-cnt CNT<-SHIFT|-v QV<-1" GETOPTS STR
IF N #equiv "1"; Z<-N; ->0
Z<-M_START N
N<-Z YDX _1 & Z<-#deb (-rN)d Z
IF N #equiv "1"; ->0
IF "1" #equiv "3" FTEST N; CNT<-0
WHILE 0<MAX--
IF QV;"N=",N
BREAKIF "1" #equiv "3" FTEST N
Q<-P_FACTOR "-cnt ", CNT, " ", (QV/"-v "), N
REPEATIF ~ "z" = 1tQ
Z<-Z," ",Q<-1dQ
MAX<-5
BREAKIF "1" #equiv N<-N%Q
WEND
IF ~N #equiv "1"; Z<-#deb Z," ",N
Z<-Z,("1" #equiv "163" FTEST Z YDX _1)/"p"

Z<-TEST_FILE STR;CNT;F;J;HA;LC;NUM;T;QV
|Test old factorizations in a file
CNT<-"50000" & QV<-0
F<-"-cnt CNT<-SHIFT|-v QV<-1" GETOPTS STR
Z<-1 #sed":r ", F & LC<-0
WHILE ~ _1 #equiv T<-#sed "@",zLC++
REPEATIF 0=rJ<-("->"#find T)/irT
HA<-JtT; T<-(J+2)dT; N<-T YDX _1
REPEATIF "p"=_1tN
"testing ", N
Z<-P_FACTORS (QV/"-v "),"-cnt ",CNT, " ", N
"Z=",Z
$T[LC-1] <- $T[LC-1]," .... ",Z
WEND

#
#  Generate a number with 2 prime factors
#

Z<-RND_PRIME N
|Random prime of N digits
IF N=1; Z<-"2357"[?4]; ->0
IF N<5; Z<-B_SIEVE N,9*N<-*/(N-1)/10; Z<-zZ[?rZ]; ->0
Z<-SCAN_FTEST RLONG N

Z<-MK_PQ N;N2;P;Q
|Number of length N with 2 factors
Z<-"6"; IF N<2; ->0
N2<-(N%2)+?2
WHILE N #ne rZ; P<-RND_PRIME N2; Q<-RND_PRIME 1 c N-N2; Z<-P*Q; WEND

#
# circle group mod M
# Zm = integers mod M
# use the multiplicative group of the ring Zm[X]/(1+X^2)
# This is a non trivial group for all M.
#

Z<-SM_POINT M;T
|Random point x,y with x^2+y^2 mod M
WHILE M #equiv Z<-"1"+T*T<-YRAND M; WEND
IF "z" = 1tZ<-Z YMOD_INV M; ->0
Z<-(("2"*Z*T) #mod M)," ",(Z*(T*T)-"1") #mod M

Z<-M SM_PLUS STR;A;B;C;D
|Add two points in circle group mod M
|Each point is an x,y pair
A<-SHIFT;B<-SHIFT;C<-SHIFT;D<-SHIFT
Z<-((A*C)-B*D) #mod M
Z<-Z," ",((A*D)+B*C) #mod M

Z<-M SM_DOUBLE STR;P;Q
|Point doubling routine on circle mod M
P<-SHIFT;Q<-SHIFT
Z<-((P*P)-Q*Q) #mod M
Z<-Z," ",("2"*P*Q) #mod M

Z<-M SM_MULT STR;K;P
|compute [K] P for circle mod M
K<-SHIFT; P<-STR; Z<-"1 0"
WHILE ~K #equiv "0"
IF ODD K; Z<-M SM_PLUS Z," ",P
K<-K%"2"; P<-M SM_DOUBLE P
WEND

Z<-SM_FACTOR STR;B;K;LST;M;P;QV;V
|Factorize by chasing circle orbits (mod M)
B<-1000 & QV<-0
M<-"-b B<-#fi SHIFT|-v QV<-1" GETOPTS STR
IF "z"=1tZ<-P<-SM_POINT M; ->0
LST<-Y_PNPP B
WHILE 0<rLST;K<-z1tLST;LST<-1dLST
IF QV; V<-K," ",P & 72tV
P<-M SM_MULT K," ",P
BREAKIF ~"1" #equiv Z<-M EA Y_PRODUCT P
WEND

Z<-SM_TEST STR;B;F;G;H;I;N;LST;SRC;QV;V
|Test circle group factoring on a file
|The file is normally output of a previous factoring run.
|the format is number->factor list
B<-"1000"; QV<-0
SRC<-"-b B<-SHIFT|-v QV<-1" GETOPTS STR<-"" #ifndef "STR"
IF 0=rSRC; SRC<-"tmp.num"
IF 0= #fstat SRC; Z<-1; ->0
|Get last componenent of hard to factor set
LST<-"";I<-0; N<-1tr:Z<- #sed"X:r",SRC,";v/p/;s/.*->//*;:cm"
WHILE I<N; LST<-LST,_CR,(#deb Z[I]) YDX _1; I++; WEND
Z<-3 #sed":d"; N<-1tr:LST<-#sstomat LST; I<-0
WHILE I<N; H<-#deb LST[I++ :]
REPEATIF "1" #equiv "3" FTEST H
"H= ",H
G<-SM_FACTOR (QV/"-v "),"-b ",B ," ",H
3 #print V<-H,"->", G & 72tV
WEND

Z<-TO_HEX N
|Large integer to hexadecimal
Z<-""
WHILE ~ "0" #equiv N
Z<-Z,"0123456789ABCDEF"[N #mod 16]
N<-N%"16"
WEND
Z<-mZ

Z<-HEX_TO_NUM STR;I
|Turn hexadecimal string to long int
Z<-"0" & STR<-UPALPHA[#av STR] & I<-0
WHILE I<rSTR;Z<-("16"*Z)+z"0123456789ABCDEF" i STR[I]; I++; WEND

Z<-RLONG2 N;K
|Number with N bits
K<-(N+3)%4      | Hex digits
Z<-HEX_TO_NUM "0123456789ABCDEF"[(4r2)b(K,4)r(-4*K) t1, ?(N-1)r2]

Z<-GSS X;Y
|Greedy sum of squares
Z<-""
WHILE ~ X #equiv "0"
Z<-Z," ",Y<-Y_SQRT X
X<-X-Y*Y
WEND
Z<-1dZ

Z<-MK_PQ2 N;K;N2;P;Q;T
|Number of with N bits and 2 factors
|Export $PQ2_FACTORS
IF N<4; Z<-"6"; ->0
N2<-(N%2)+?2
WHILE 1
P<-SCAN_FTEST  Z<-RLONG2 N2
Q<-SCAN_FTEST  Z<-RLONG2 N-N2
Z<-P*Q
K<-rT<-TO_HEX Z
K<-(4*K-1)+0 1 2 2 3 3 3 3 4 4 4 4 4 4 4 4["0123456789ABCDEF"iT[0]]
BREAKIF K=N
WEND
$PQ2_FACTORS<-P," ",Q

Z<-MSIEVE STR
|MSIEVE: A Library for Factoring Large Integers
|Jason Papadopoulos (jasonp@boo.net) 11/27/05
Z<-#sstomat _CR, #cmd"< ../msieve-1.03/msieve -q ", STR

Z<-Y_HEXTON STR;T
|Convert hex longs to decimal longs
Z<-""
WHILE 0<rT<-SHIFT;Z<-Z," ",HEX_TO_NUM T;WEND
Z<-#deb Z

Z<-J Y_COLUMN X;U
|Get Jth column of ' ' delimited matrix
U<-0 _1 d 0,y+\y" "=X<-X," "
Z<- #sstomat _1 m (J=,U)/,X