Greatest prime dividing the n-th cubefree number: Difference between revisions
Greatest prime dividing the n-th cubefree number (view source)
Revision as of 00:28, 6 March 2024
, 3 months ago→{{header|Phix}}
(Added FreeBASIC) |
|||
Line 109:
The 1000000th term of a[n] is 1202057
The 10000000th term of a[n] is 1202057</pre>
=={{header|Phix}}==
Quite possibly flawed, but it does finish in the blink of an eye.
<syntaxhighlight lang="phix">
with javascript_semantics
-- Note this routine had a major hiccup at 27000 = 2^3*3^3*5^3 and another was predicted at 9261000 = 27000*7^3.
-- The "insufficient kludge" noted below makes it match Wren over than limit, but quite possibly only by chance.
function cubes_before(atom n)
-- nb: if n is /not/ cube-free it /is/ included in the result.
-- nth := n-cubes_before(n) means n is the nth cube-free integer,
-- but only if cubes_before(n-1)==cubes_before(n),
-- otherwise cubes_before(cubicate) isn't really very meaningful.
integer r = 0
bool xpm = true -- extend prior multiples
sequence pm = {}
for i=1 to n do -- (or eg floor(cbrt(n)))
atom p3 = power(get_prime(i),3)
if p3>n then exit end if
integer k = floor(n/p3)
for mask=1 to power(2,length(pm))-1 do
integer m = mask, mi = 0
atom kpm = p3
while m do
mi += 1
if odd(m) then
kpm *= pm[mi]
end if
m = floor(m/2)
end while
if kpm>n then
if count_bits(mask)=1 then
xpm = false
pm = pm[1..mi-1]
exit
end if
else
-- subtract? already counted multiples..
integer l = floor(n/kpm)
-- DEV insufficient kludge... (probably)
--if count_bits(mask)=1 then
if odd(count_bits(mask)) then -- match Wren
k -= l
else
k += l
end if
end if
end for
r += k
if xpm then
pm &= p3
end if
end for
return r
end function
function cube_free(atom nth)
-- get the nth cube-free integer
atom lo = nth, hi = lo*2, mid, cb, k
while hi-cubes_before(hi)<nth do
lo = hi
hi = lo*2
end while
-- bisect until we have a possible...
while true do
mid = floor((lo+hi)/2)
cb = cubes_before(mid)
k = mid-cb
if k=nth then
-- skip omissions
while cubes_before(mid-1)!=cb do
mid -= 1
cb -= 1
end while
exit
elsif k<nth then
lo = mid
else
hi = mid
end if
end while
return mid
end function
function A370833(atom n)
if n=1 then return 1 end if
atom nth = cube_free(n)
sequence f = prime_powers(nth)
return f[$][1]
end function
atom t0 = time()
printf(1,"First 100 terms of a[n]:\n%s\n",join_by(apply(tagset(100),A370833),1,10," ",fmt:="%3d"))
for n in sq_power(10,tagset(7,3)) do
printf(1,"The %,dth term of a[n] is %,d.\n",{n,A370833(n)})
end for
?elapsed(time()-t0)
</syntaxhighlight>
{{out}}
<pre>
First 100 terms of a[n]:
1 2 3 2 5 3 7 3 5 11
3 13 7 5 17 3 19 5 7 11
23 5 13 7 29 5 31 11 17 7
3 37 19 13 41 7 43 11 5 23
47 7 5 17 13 53 11 19 29 59
5 61 31 7 13 11 67 17 23 7
71 73 37 5 19 11 13 79 41 83
7 17 43 29 89 5 13 23 31 47
19 97 7 11 5 101 17 103 7 53
107 109 11 37 113 19 23 29 13 59
The 1,000th term of a[n] is 109.
The 10,000th term of a[n] is 101.
The 100,000th term of a[n] is 1,693.
The 1,000,000th term of a[n] is 1,202,057.
The 10,000,000th term of a[n] is 1,202,057.
"0.1s"
</pre>
=={{header|RPL}}==
| |||