Greatest prime dividing the n-th cubefree number: Difference between revisions

Content added Content deleted
m (→‎Version 2: Added a libheader.)
m (→‎{{header|Phix}}: extended to 1e11, matching the Pascal output)
Line 583: Line 583:
</pre>
</pre>
=={{header|Phix}}==
=={{header|Phix}}==
Quite possibly flawed, but it does finish in the blink of an eye.
Quite possibly flawed, but it does at least complete to 1e9 in the blink of an eye, and 93s for 1e11.
<syntaxhighlight lang="phix">
<syntaxhighlight lang="phix">
with javascript_semantics
with javascript_semantics
Line 595: Line 595:
-- but only if cubes_before(n-1)==cubes_before(n),
-- but only if cubes_before(n-1)==cubes_before(n),
-- otherwise cubes_before(cubicate) isn't really very meaningful.
-- otherwise cubes_before(cubicate) isn't really very meaningful.
integer r = 0
atom r = 0
bool xpm = true -- extend prior multiples
bool xpm = true -- extend prior multiples
sequence pm = {}
sequence pm = {}
for i=1 to n do -- (or eg floor(cbrt(n)))
for i=1 to floor(power(n,1/3)) do
atom p3 = power(get_prime(i),3)
atom p3 = power(get_prime(i),3)
if p3>n then exit end if
if p3>n then exit end if
integer k = floor(n/p3)
integer k = floor(n/p3)
for mask=1 to power(2,length(pm))-1 do
for mask=1 to power(2,length(pm))-1 do
integer m = mask, mi = 0
integer m = mask, mi = 0, bc = count_bits(mask)
atom kpm = p3
atom kpm = p3
while m do
while m do
Line 613: Line 613:
end while
end while
if kpm>n then
if kpm>n then
if count_bits(mask)=1 then
if bc=1 then
xpm = false
xpm = false
pm = pm[1..mi-1]
pm = pm[1..mi-1]
Line 622: Line 622:
integer l = floor(n/kpm)
integer l = floor(n/kpm)
-- DEV insufficient kludge... (probably)
-- DEV insufficient kludge... (probably)
--if count_bits(mask)=1 then
--if bc=1 then
if odd(count_bits(mask)) then -- match Wren
if odd(bc) then
k -= l
k -= l
else
else
Line 646: Line 646:
end while
end while
-- bisect until we have a possible...
-- bisect until we have a possible...
atom t1 = time()+1
while true do
while true do
mid = floor((lo+hi)/2)
mid = floor((lo+hi)/2)
Line 661: Line 662:
else
else
hi = mid
hi = mid
end if
if time()>t1 then
progress("bisecting %,d..%,d...\r",{lo,hi})
t1 = time()+1
end if
end if
end while
end while
Line 667: Line 672:


function A370833(atom n)
function A370833(atom n)
if n=1 then return 1 end if
if n=1 then return {1,1,1} end if
atom nth = cube_free(n)
atom nth = cube_free(n)
sequence f = prime_powers(nth)
sequence f = prime_powers(nth)
return f[$][1]
return {n,nth,f[$][1]}
end function
end function


atom t0 = time()
atom t0 = time()
printf(1,"First 100 terms of a[n]:\n%s\n",join_by(apply(tagset(100),A370833),1,10," ",fmt:="%3d"))
sequence f100 = vslice(apply(tagset(100),A370833),3)
printf(1,"First 100 terms of a[n]:\n%s\n",join_by(f100,1,10," ",fmt:="%3d"))
for n in sq_power(10,tagset(7,3)) do
for n in sq_power(10,tagset(11,3)) do
printf(1,"The %,dth term of a[n] is %,d.\n",{n,A370833(n)})
printf(1,"The %,dth term of a[n] is %,d with highest divisor %,d.\n",A370833(n))
end for
end for
?elapsed(time()-t0)
?elapsed(time()-t0)
Line 694: Line 700:
107 109 11 37 113 19 23 29 13 59
107 109 11 37 113 19 23 29 13 59


The 1,000th term of a[n] is 109.
The 1,000th term of a[n] is 1,199 with highest divisor 109.
The 10,000th term of a[n] is 101.
The 10,000th term of a[n] is 12,019 with highest divisor 101.
The 100,000th term of a[n] is 1,693.
The 100,000th term of a[n] is 120,203 with highest divisor 1,693.
The 1,000,000th term of a[n] is 1,202,057.
The 1,000,000th term of a[n] is 1,202,057 with highest divisor 1,202,057.
The 10,000,000th term of a[n] is 1,202,057.
The 10,000,000th term of a[n] is 12,020,570 with highest divisor 1,202,057.
The 100,000,000th term of a[n] is 120,205,685 with highest divisor 20,743.
"0.1s"
The 1,000,000,000th term of a[n] is 1,202,056,919 with highest divisor 215,461.
The 10,000,000,000th term of a[n] is 12,020,569,022 with highest divisor 1,322,977.
The 100,000,000,000th term of a[n] is 120,205,690,298 with highest divisor 145,823.
"1 minute and 33s"
</pre>
</pre>


Retrieved from "https://rosettacode.org/wiki/Greatest_prime_dividing_the_n-th_cubefree_number"