Rosetta Code

Square form factorization: Difference between revisions

← Older editNewer edit →
Content deleted Content added
mNo edit summary
Petelomax (talk | contribs)
7,878 edits
→‎{{header|Phix}}: replaced with translation of C
Line 1,052:
 
=={{header|Phix}}==
{{trans|C|<small>(fixes the two incorrectly failing cases of the previous version)</small>}}
{{trans|Julia|and wikipedia}}
An extra couple of the high 19-digit tests fail compared to other entries, but I can live with that. 0 moved up to mark the 32-bit limit.
<lang Phix>--requires(64) -- (decided to limit 32-bit explicitly instead)
constant MxN = power(2,63),
constant multiplier = {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},
INT_MAXm = power(2{1,iff(machine_bits()=32?53:64)) 3, 5, 7, 11}
function squfof(atom N)
-- square form factorization
 
integer h, a=0, b, c, u=0, v, w, rN, q, r, t
if remainder(N,2)==0 then return 2 end if
h = floor(sqrt(N) + 0.5)
if h*h==N then return h end if
for k=1 to length(m) do
integer mk = m[k]
if mk>1 and remainder(N,mk)==0 then return mk end if
//check overflow m * N
if N>MxN/mk then exit end if
atom mN = N*mk
r = floor(sqrt(mN))
if r*r>mN then r -= 1 end if
rN = r
//principal form
{b,a,c} = {r,1,a}
h = floor((rN+b)/a)*a-b
c = floor((mN-h*h)/a)
for i=2 to floor(sqrt(2*r)) * 4-1 do
//search principal cycle
{a,c,t} = {c,a,b}
q = floor((rN+b)/a)
b = q*a-b
c += q*(t-b)
 
if remainder(i,2)==0 then
r = floor(sqrt(c)+0.5)
if r*r==c then
//square form found
//inverse square root
q = floor((rN-b)/r)
v = q*r+b
w = floor((mN-v*v)/r)
//search ambiguous cycle
u = r
while true do
{u,w,r} = {w,u,v}
q = floor((rN+v)/u)
v = q*u-v
if v==r then exit end if
w += q*(r-v)
end while
 
//symmetry point
function square_form_factor(atom n)
atom s h = round(sqrtgcd(n)N,u)
if s * s == n if h!=1 then return sh end if
end if
for mi=1 to length(multiplier) do
atom k = multiplier[mi] end if
if n > INT_MAX/k then exit end iffor
end for
atom d = k * n,
return 1
p0 = floor(sqrt(d)),
pprev = p0, p = p0,
qprev = 1,
Q = d - p0 * p0,
l = floor(2 * sqrt(2 * s)),
B = 3*l,
i = 2, r, b, q
while i < B do
b = floor((p0 + p) / Q)
p = b * Q - p
q = Q
Q = qprev + b * (pprev - p)
r = round(sqrt(Q))
if remainder(i,2)!=1 and r * r == Q then exit end if
qprev = q
pprev = p
i += 1
end while
if i < B then
b = floor((p0 - p) / r)
p = b * r + p
pprev = p
qprev = r
Q = floor((d - pprev * pprev) / qprev)
i = 0
while true do
b = floor((p0 + p) / Q)
pprev = p
p = b * Q - p
q = Q
Q = qprev + b * (pprev - p)
qprev = q
i += 1
if p == pprev then exit end if
end while
r = gcd(n, qprev)
if r != 1 and r != n then return r end if
end if
end for
return 0
end function
constant tests = {2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 132900595214317, 4285444720834839, 223553581,
202765128123515517, 1111111111133409583, 10089559816944524219, 100274262802113290059, 60012462237239223553581, 28712952341479142854447, 223553581, 2027651281,
011111111111, 100895598169, 1002742628021, -- (limitprime/expected ofto 32-bitfail)
900719925474093160012462237239, 11111111111111111287129523414791, 3141592653589793239007199254740931, 38430716820228150732, 419244183493398773,-- (limit of 32-bit)
11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773,
658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799,
1537228672809128917, 4611686018427387877}
 
printf(1,"IntegerN Factor Quotient f N/f\n")
printf(1,"======================================\n")
for i=1 to length(tests) do
atom nN = tests[i],
if N=32 then
factr = square_form_factor(n)
if machine_bits()=32 then exit end if
string f = iff(factr=0?"fail":sprintf("%d",n/factr))
else
printf(1,"%-22d %-10d %s\n",{n,factr,f})
atom f = squfof(N)
if machine_bits()=32 and n=0 then exit end if
printf(1,"%-22d %s\n",{N,iff(f=1?"fail":sprintf("%-10d %d",{f,N/f}))})
end if
end for</lang>
{{out}}
<small>(on 64-bit, whereas the last 10 entries, ie 11111111111111111 on, are simply omitted on 32-bit)</small>
<pre>
IntegerN Factor Quotient f N/f
======================================
2501 61 41
12851 71 181
13289 13797 97 137
75301 293 257
120787 43 2809
Line 1,135 ⟶ 1,152:
997417 257 3881
7091569 2663 2663
5214317 73 71429
20834839 3361 6199
23515517 53 443689
33409583 991 33713
44524219 593 75083
13290059 3119 4261
42854447 9689 4423
223553581 11213 19937
202765128142854447 46061 4423 44021 9689
223553581 11213 19937
2027651281 44021 46061
11111111111 21649 513239
100895598169 112303 898423
1002742628021 0 fail
60012462237239 6862753 8744663
287129523414791 6059887 47381993
0 0 fail
9007199254740931 10624181 847801751
11111111111111111 2071723 5363222357
Line 1,152 ⟶ 1,174:
658812288346769681 62222119 10588072199
922337203685477563 110075821 8379108103
1000000000000000127 0111756107 fail8948056861
1152921505680588799 139001459 8294312261
1537228672809128917 026675843 fail57626245319
4611686018427387877 343242169 13435662733
</pre>
Retrieved from "https://rosettacode.org/wiki/Square_form_factorization"