54 lines
4.3 KiB
Rexx
54 lines
4.3 KiB
Rexx
/*REXX pgm factors an integer using Daniel Shanks' (1917-1996) square form factorization*/
|
|
numeric digits 100 /*ensure enough decimal digits.*/
|
|
call dMults 1,3,5,7,11,3*5,3*7,3*11,5*7,5*11,7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11
|
|
call dTests 2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, ,
|
|
42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021, ,
|
|
60012462237239, 287129523414791, 9007199254740931, 11111111111111111, ,
|
|
314159265358979323, 384307168202281507, 419244183493398773, ,
|
|
658812288346769681, 922337203685477563, 1000000000000000127, ,
|
|
1152921505680588799, 1537228672809128917, 4611686018427387877
|
|
w= length( commas(!.$) ) /*the max width of test numbers*/
|
|
do tests=1 for !.0; n= !.tests; nc= commas(n)
|
|
f= ssff(n); fc= commas(f); wf= length(fc); if f\==0 then nf= commas(n%f)
|
|
if f\==0 then do; nfc= commas(n%f); wnfc= length(nfc); end
|
|
if f ==0 then _= " (Shank's square form factor failed.)"
|
|
else _= ' factors are: ' right( fc, max(w%2 , wf ) ) " and " ,
|
|
right(nfc, max(w%2+4, wnfc) )
|
|
say right(nc, w+5) _
|
|
end /*tests*/
|
|
exit 0 /*stick a fork in it, we're all done. */
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ?
|
|
dMults: @.$= 0; do j=1 for arg(); @.j= arg(j); @.$=max(@.$, @.j); end; @.0=j-1; return
|
|
dTests: !.$= 0; do j=1 for arg(); !.j= arg(j); !.$=max(!.$, !.j); end; !.0=j-1; return
|
|
gcd: procedure; parse arg x,y; do until _==0; _= x // y; x= y; y= _; end; return x
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
iSqrt: procedure; parse arg x; r=0; q=1; do while q<=x; q=q*4; end
|
|
do while q>1; q=q%4; _=x-r-q; r=r%2; if _>=0 then do;x=_;r=r+q; end; end
|
|
return r
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
ssff: procedure expose @.; parse arg n; n= abs(n); er= '***error***'
|
|
s= iSqrt(n); if s**2==n then return s; big= 2**digits()
|
|
do #=1 for @.0; k= @.# /*get a # from the list of low factors*/
|
|
if n>big/k then do; say er 'number is too large: ' commas(k); exit 8; end
|
|
d= n*k; po= iSqrt(d); p= po
|
|
pprev= po; QQ= d - po*po
|
|
qprev= 1; BB= iSqrt(s+s)*6
|
|
do i=2 while i<BB; b= (po+p)%QQ
|
|
p= b*QQ - p; q= QQ
|
|
QQ= qprev + b*(pprev-p); r= iSqrt(QQ)
|
|
if i//2==0 then if r*r==QQ then leave
|
|
qprev= q; pprev= p
|
|
end /*i*/
|
|
if i>=BB then iterate
|
|
b= (po-p)%r; p= b*r + p
|
|
pprev= p; qprev= r
|
|
QQ= (d - pprev*pprev)%qprev
|
|
do until p==pprev; pprev= p
|
|
b= (po+p)%QQ; q= QQ; p= b*QQ - p
|
|
QQ= qprev + b*(pprev-p); qprev= q
|
|
end /*until*/
|
|
r= gcd(n, qprev)
|
|
if r\==1 then if r\==n then return r
|
|
end /*#*/
|
|
return 0
|