Data update
This commit is contained in:
Ingy döt Net
parent
72eb4943cb
commit
4d5544505c
2347 changed files with 62432 additions and 16731 deletions
|
|
@ -0,0 +1,110 @@
|
|||
from math import sqrt, gcd
|
||||
|
||||
M = [1, 3, 5, 7, 11]
|
||||
UINT64_MAX = 18446744073709551615
|
||||
|
||||
def isqrt(n: int):
|
||||
return int(sqrt(n))
|
||||
|
||||
def squfof(n: int):
|
||||
if n % 2 == 0:
|
||||
return 2
|
||||
|
||||
h = int(sqrt(n) + 0.5)
|
||||
|
||||
if h ** 2 == n:
|
||||
return h
|
||||
|
||||
for m in M:
|
||||
if m > 1 and n % m == 0:
|
||||
return m
|
||||
|
||||
if n > UINT64_MAX / m:
|
||||
break
|
||||
|
||||
mn = m * n
|
||||
r = isqrt(mn)
|
||||
|
||||
if r ** 2 > mn:
|
||||
r -= 1
|
||||
|
||||
rn = r
|
||||
|
||||
b = r
|
||||
a = 1
|
||||
h = ((rn + b) // a) * a - b
|
||||
c = (mn - h * h) // a
|
||||
|
||||
for i in range(2, 4 * isqrt(2 * r)):
|
||||
a, c = c, a
|
||||
q = (rn + b) // a
|
||||
t = b
|
||||
b = q * a - b
|
||||
c += q * (t - b)
|
||||
|
||||
if i % 2 == 0:
|
||||
r = int(sqrt(c) + 0.5)
|
||||
|
||||
if r ** 2 == c:
|
||||
q = (rn - b) // r
|
||||
v = q * r + b
|
||||
w = (mn - v * v) // r
|
||||
|
||||
u = r
|
||||
|
||||
while True:
|
||||
w, u = u, w
|
||||
r = v
|
||||
q = (rn + v) // u
|
||||
v = q * u - v
|
||||
|
||||
if v == r:
|
||||
break
|
||||
|
||||
w += q * (r - v)
|
||||
|
||||
h = gcd(n, u)
|
||||
|
||||
if h != 1:
|
||||
return h
|
||||
|
||||
return 1
|
||||
|
||||
data = [
|
||||
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
|
||||
]
|
||||
|
||||
print("N f N/f")
|
||||
print("======================================")
|
||||
|
||||
for n in data:
|
||||
f = squfof(n)
|
||||
res = 'fail' if f == 1 else f'{f:<10} {n // f}'
|
||||
print(f'{n:<22} {res}')
|
||||
|
|
@ -1,54 +1,94 @@
|
|||
/*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
|
||||
-- 8 May 2025
|
||||
include Settings
|
||||
numeric digits 40
|
||||
|
||||
say 'SQUARE FORM FACTORIZATION'
|
||||
say version
|
||||
say
|
||||
call TestNumbers
|
||||
call Selected
|
||||
call Randomized
|
||||
exit
|
||||
|
||||
TestNumbers:
|
||||
procedure expose test.
|
||||
t = '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'
|
||||
do i = 1 to Words(t)
|
||||
test.i = Word(t,i)
|
||||
end
|
||||
test.0 = i-1
|
||||
return
|
||||
|
||||
Selected:
|
||||
procedure expose mult. test. glob.
|
||||
call Time('r')
|
||||
say 'Find a factor for 28 selected numbers...'
|
||||
do t = 1 for test.0
|
||||
n = test.t
|
||||
call Time('r')
|
||||
f = Squfof(n)
|
||||
if f = 0 then
|
||||
u = 'failed'
|
||||
else
|
||||
u = f 'x' n/f
|
||||
say Format(Time('e'),3,3)'s Squfof ' n '('Xpon(n)+1 'digits) =' u
|
||||
call Time('r')
|
||||
f = Pollardrho(n)
|
||||
if f = 0 then
|
||||
u = 'failed'
|
||||
else
|
||||
u = f 'x' n/f
|
||||
say Format(Time('e'),3,3)'s Pollard' n '('Xpon(n)+1 'digits) =' u
|
||||
if n < 1e17 then do
|
||||
call Time('r')
|
||||
f = Trialdiv(n); u = f 'x' n/f
|
||||
say Format(Time('e'),3,3)'s Trial ' n '('Xpon(n)+1 'digits) =' u
|
||||
end
|
||||
say
|
||||
end
|
||||
say
|
||||
return
|
||||
|
||||
Randomized:
|
||||
procedure expose mult. glob.
|
||||
say 'Find a factor for 28 random numbers...'
|
||||
x = 0
|
||||
do until x = 28
|
||||
n = (Right(Randu(),Randu(1,10))||Right(Randu(),Randu(1,10))||Right(Randu(),Randu(1,10)))/1
|
||||
if Pos(Right(n,1),05) > 0 | Even(n) | Digitsum(n)//3 = 0 then
|
||||
iterate
|
||||
call Time('r')
|
||||
f = Squfof(n)
|
||||
if f = 0 then
|
||||
u = 'failed'
|
||||
else
|
||||
u = f 'x' n/f
|
||||
say Format(Time('e'),3,3)'s Squfof ' n '('Xpon(n)+1 'digits) =' u
|
||||
if n < 1e24 then do
|
||||
call Time('r')
|
||||
f = Pollardrho(n)
|
||||
if f = 0 then
|
||||
u = 'failed'
|
||||
else
|
||||
u = f 'x' n/f
|
||||
say Format(Time('e'),3,3)'s Pollard' n '('Xpon(n)+1 'digits) =' u
|
||||
end
|
||||
if n < 1e17 then do
|
||||
call Time('r')
|
||||
f = Trialdiv(n); u = f 'x' n/f
|
||||
say Format(Time('e'),3,3)'s Trial ' n '('Xpon(n)+1 'digits) =' u
|
||||
end
|
||||
say
|
||||
x = x+1
|
||||
end
|
||||
say
|
||||
return
|
||||
|
||||
include Numbers
|
||||
include Functions
|
||||
include Abend
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue