Powerful numbers: Difference between revisions

Content added Content deleted
Line 412: Line 412:
9-powerful numbers <= 10ⁿ (where n == 0 through 18): 1, 1, 1, 2, 6, 11, 16, 24, 38, 59, 94, 145, 217, 317, 453, 644, 919, 1308, 1868
9-powerful numbers <= 10ⁿ (where n == 0 through 18): 1, 1, 1, 2, 6, 11, 16, 24, 38, 59, 94, 145, 217, 317, 453, 644, 919, 1308, 1868
10-powerful numbers <= 10ⁿ (where n == 0 through 19): 1, 1, 1, 1, 5, 9, 14, 21, 28, 43, 68, 104, 155, 227, 322, 447, 621, 858, 1192, 1651</pre>
10-powerful numbers <= 10ⁿ (where n == 0 through 19): 1, 1, 1, 1, 5, 9, 14, 21, 28, 43, 68, 104, 155, 227, 322, 447, 621, 858, 1192, 1651</pre>

=={{header|Phix}}==
{{libheader|mpfr}}
===priority queue===
Tried a slightly different approach, but it gets 7/10 at best on performance.
<lang Phix>include mpfr.e

integer pq = NULL -- data is the prime powers, eg {2,2} for 2^2*3^2
-- priority is the mpz of that, eg/ie 36

procedure padd(mpz n, np, sequence s)
if pq=NULL then
pq = pq_new(MIN_HEAP,routine_id("mpz_cmp"))
end if
if mpz_cmp(np,n)<=0 then
pq_add({s,mpz_init_set(np)},pq)
end if
end procedure

procedure add_next(integer k, mpz n, p, sequence s)
mpz np = mpz_init()
integer l = length(s)
for i=1 to l do
integer ki = iff(s[i]=0?k:1)
mpz_ui_pow_ui(np,get_prime(i),ki)
mpz_mul(np,np,p)
s[i] += ki
padd(n,np,s)
s[i] -= ki
end for
mpz_ui_pow_ui(np,get_prime(l+1),k)
if l>0 and s[l]=k and s[1..l-1]==repeat(0,l-1) then
padd(n,np,0&s)
elsif s={} then
mpz_mul(np,np,p)
padd(n,np,s&k)
end if
end procedure

function powerful(integer k, mpz n)
mpz prevp = mpz_init(1)
sequence res = {prevp}
add_next(k,n,prevp,{})
while pq_size(pq) do
{sequence s, mpz p} = pq_pop(pq)
if mpz_cmp(p,prevp)!=0 then
res &= p
add_next(k,n,p,s)
prevp = p
end if
end while
return res
end function
mpz p = mpz_init(10)
for k=2 to 10 do
mpz_mul_si(p,p,10)
sequence s = powerful(k,p)
integer l = length(s)
for i=1 to 5 do
s[i] = mpz_get_str(s[i])
s[-i] = mpz_get_str(s[-i])
end for
s[6..-6] = {"..."}
s = join(s,",")
printf(1,"%d %2d-powerful numbers <= 10^%-2d : %s\n",{l,k,k,s})
end for

puts(1,"\n")
--Not really fast enough, so only doing to k+5..k+8 for the first 4, but k+9 for last 5
for k=2 to 10 do
sequence counts = {}
integer lim = k+min(k+3,9)
mpz_set_si(p,1)
for j=0 to lim do
printf(1,"counting (%d/%d)...\r",{j,lim})
counts &= length(powerful(k,p))
mpz_mul_si(p,p,10)
end for
printf(1,"Counts of %2d-powerful numbers <= 10^(0..%d) : %v\n",{k,lim,counts})
end for</lang>
{{out}}
<pre>
14 2-powerful numbers <= 10^2 : 1,4,8,9,16,...,49,64,72,81,100
20 3-powerful numbers <= 10^3 : 1,8,16,27,32,...,625,648,729,864,1000
25 4-powerful numbers <= 10^4 : 1,16,32,64,81,...,5184,6561,7776,8192,10000
32 5-powerful numbers <= 10^5 : 1,32,64,128,243,...,65536,69984,78125,93312,100000
38 6-powerful numbers <= 10^6 : 1,64,128,256,512,...,559872,746496,823543,839808,1000000
46 7-powerful numbers <= 10^7 : 1,128,256,512,1024,...,7558272,8388608,8957952,9765625,10000000
52 8-powerful numbers <= 10^8 : 1,256,512,1024,2048,...,60466176,67108864,80621568,90699264,100000000
59 9-powerful numbers <= 10^9 : 1,512,1024,2048,4096,...,644972544,725594112,816293376,967458816,1000000000
68 10-powerful numbers <= 10^10 : 1,1024,2048,4096,8192,...,7739670528,8589934592,8707129344,9795520512,10000000000

Counts of 2-powerful numbers <= 10^(0..7) : {1,4,14,54,185,619,2027,6553}
Counts of 3-powerful numbers <= 10^(0..9) : {1,2,7,20,51,129,307,713,1645,3721}
Counts of 4-powerful numbers <= 10^(0..11) : {1,1,5,11,25,57,117,235,464,906,1741,3312}
Counts of 5-powerful numbers <= 10^(0..13) : {1,1,3,8,16,32,63,117,211,375,659,1153,2000,3402}
Counts of 6-powerful numbers <= 10^(0..15) : {1,1,2,6,12,21,38,70,121,206,335,551,900,1451,2326,3706}
Counts of 7-powerful numbers <= 10^(0..16) : {1,1,1,4,10,16,26,46,77,129,204,318,495,761,1172,1799,2740}
Counts of 8-powerful numbers <= 10^(0..17) : {1,1,1,3,8,13,19,32,52,85,135,211,315,467,689,1016,1496,2191}
Counts of 9-powerful numbers <= 10^(0..18) : {1,1,1,2,6,11,16,24,38,59,94,145,217,317,453,644,919,1308,1868}
Counts of 10-powerful numbers <= 10^(0..19) : {1,1,1,1,5,9,14,21,28,43,68,104,155,227,322,447,621,858,1192,1651}
</pre>
===translation===
{{trans|Go}} {{trans|Sidef}}
Significantly faster. The last few counts were wrong when using native ints/atoms, so I went gmp, which means it also works on 32-bit.
<lang Phix>include mpfr.e
function powerful(mpz n, integer k, sequence res={}, mpz m=NULL, integer r=0)
if m=NULL then
m = mpz_init(1)
r = 2*k-1
end if
if r<k then
res = append(res,iff(mpz_fits_atom(m)?mpz_get_atom(m):mpz_get_str(m)))
else
mpz t = mpz_init()
mpz_fdiv_q(t, n, m) -- t := floor(n/m)
{} = mpz_root(t,t,r) -- t = floor(rth root of t)
if not mpz_fits_integer(t) then ?9/0 end if
integer lim = mpz_get_integer(t)
for v=1 to lim do
if r<=k
or (square_free(v) and mpz_gcd_ui(NULL,m,v)=1) then
mpz_ui_pow_ui(t,v,r) -- t:= v^r
mpz_mul(t,t,m) -- t*= m
res = powerful(n, k, res, t, r-1)
end if
end for
end if
return res
end function

function powercount(mpz n, integer k, mpz m=NULL, integer r=0)
atom res = 0
if m=NULL then
m = mpz_init(1)
r = 2*k-1
end if
mpz t = mpz_init()
mpz_fdiv_q(t, n, m) -- t := floor(n/m)
{} = mpz_root(t,t,r) -- t = floor(rth root of t)
if not mpz_fits_integer(t) then ?9/0 end if
integer lim = mpz_get_integer(t)
if r<=k then
res += lim
else
for v=1 to lim do
if square_free(v) and mpz_gcd_ui(NULL,m,v)=1 then
mpz_ui_pow_ui(t,v,r) -- t:= v^r
mpz_mul(t,t,m) -- t*= m
res += powercount(n, k, t, r-1)
end if
end for
end if
return res
end function
atom t0 = time()
mpz p = mpz_init(10)
for k=2 to 10 do
mpz_mul_si(p,p,10)
sequence s = sort(powerful(p, k))
integer l = length(s)
s[6..-6] = {"..."}
s = ppf(s,{pp_FltFmt,"%d",pp_StrFmt,-1,pp_IntCh,false,pp_Maxlen,100})
printf(1,"%d %2d-powerful numbers <= 10^%-2d : %s\n",{l,k,k,s})
end for
puts(1,"\n")
for k=2 to 10 do
mpz_set_si(p,1)
sequence counts = {}
for j=0 to k+9 do
counts = append(counts, powercount(p, k))
mpz_mul_si(p,p,10)
end for
printf(1,"Counts of %2d-powerful numbers <= 10^(0..%d) : %v\n",{k,k+9,counts})
end for
?elapsed(time()-t0)</lang>
{{out}}
<pre>
14 2-powerful numbers <= 10^2 : {1,4,8,9,16, ..., 49,64,72,81,100}
20 3-powerful numbers <= 10^3 : {1,8,16,27,32, ..., 625,648,729,864,1000}
25 4-powerful numbers <= 10^4 : {1,16,32,64,81, ..., 5184,6561,7776,8192,10000}
32 5-powerful numbers <= 10^5 : {1,32,64,128,243, ..., 65536,69984,78125,93312,100000}
38 6-powerful numbers <= 10^6 : {1,64,128,256,512, ..., 559872,746496,823543,839808,1000000}
46 7-powerful numbers <= 10^7 : {1,128,256,512,1024, ..., 7558272,8388608,8957952,9765625,10000000}
52 8-powerful numbers <= 10^8 : {1,256,512,1024,2048, ..., 60466176,67108864,80621568,90699264,100000000}
59 9-powerful numbers <= 10^9 : {1,512,1024,2048,4096, ..., 644972544,725594112,816293376,967458816,1000000000}
68 10-powerful numbers <= 10^10 : {1,1024,2048,4096,8192, ..., 7739670528,8589934592,8707129344,9795520512,10000000000}

Counts of 2-powerful numbers <= 10^(0..11) : {1,4,14,54,185,619,2027,6553,21044,67231,214122,680330}
Counts of 3-powerful numbers <= 10^(0..12) : {1,2,7,20,51,129,307,713,1645,3721,8348,18589,41136}
Counts of 4-powerful numbers <= 10^(0..13) : {1,1,5,11,25,57,117,235,464,906,1741,3312,6236,11654}
Counts of 5-powerful numbers <= 10^(0..14) : {1,1,3,8,16,32,63,117,211,375,659,1153,2000,3402,5770}
Counts of 6-powerful numbers <= 10^(0..15) : {1,1,2,6,12,21,38,70,121,206,335,551,900,1451,2326,3706}
Counts of 7-powerful numbers <= 10^(0..16) : {1,1,1,4,10,16,26,46,77,129,204,318,495,761,1172,1799,2740}
Counts of 8-powerful numbers <= 10^(0..17) : {1,1,1,3,8,13,19,32,52,85,135,211,315,467,689,1016,1496,2191}
Counts of 9-powerful numbers <= 10^(0..18) : {1,1,1,2,6,11,16,24,38,59,94,145,217,317,453,644,919,1308,1868}
Counts of 10-powerful numbers <= 10^(0..19) : {1,1,1,1,5,9,14,21,28,43,68,104,155,227,322,447,621,858,1192,1651}
"0.8s"
</pre>


=={{header|Sidef}}==
=={{header|Sidef}}==
Retrieved from "https://rosettacode.org/wiki/Powerful_numbers"