Greatest prime dividing the n-th cubefree number
A cubefree number is a positive integer whose prime factorization does not contain any third (or higher) power factors. If follows that all primes are trivially cubefree and the first cubefree number is 1 because it has no prime factors.
Greatest prime dividing the n-th cubefree number is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
- Definitions
Let a[n] be the greatest prime dividing the n-th cubefree number for n >= 2. By convention, let a[1] = 1 even though the first cubefree number, 1, has no prime factors.
- Examples
a[2] is clearly 2 because it is the second cubefree number and also prime. The fourth cubefree number is 4 and it's highest prime factor is 2, hence a[4] = 2.
- Task
Compute and show on this page the first 100 terms of a[n]. Also compute and show the 1,000th, 10,000th and 100,000th members of the sequence.
- Stretch
Compute and show the 1 millionth and 10 millionth terms of the sequence.
This may take a while for interpreted languages.
- Reference
FreeBASIC
Without using external libraries, it takes about 68 seconds to run on my system (Core i5).
Dim Shared As Integer factors()
Dim As Integer res(101)
res(0) = 1
Dim As Integer count = 1
Dim As Integer j, i = 2
Dim As Integer lim1 = 100
Dim As Integer lim2 = 1000
Dim As Integer max = 1e7
Dim As Integer cubeFree = 0
Sub primeFactors(n As Integer, factors() As Integer)
Dim As Integer i = 2, cont = 2
While (i * i <= n)
While (n Mod i = 0)
n /= i
cont += 1
Redim Preserve factors(1 To cont)
factors(cont) = i
Wend
i += 1
Wend
If (n > 1) Then
cont += 1
Redim Preserve factors(1 To cont)
factors(cont) = n
End If
End Sub
While (count < max)
primeFactors(i, factors())
If (Ubound(factors) < 3) Then
cubeFree = 1
Else
cubeFree = 1
For j = 2 To Ubound(factors)
If (factors(j-2) = factors(j-1) And factors(j-1) = factors(j)) Then
cubeFree = 0
Exit For
End If
Next j
End If
If (cubeFree = 1) Then
If (count < lim1) Then
res(count) = factors(Ubound(factors))
End If
count += 1
If (count = lim1) Then
Print "First "; lim1; " terms of a[n]:"
For j = 1 To lim1
Print Using "####"; res(j-1);
If (j Mod 10 = 0) Then Print
Next j
Print
Elseif (count = lim2) Then
Print "The"; count; " term of a[n] is"; factors(Ubound(factors))
lim2 *= 10
End If
End If
i += 1
Wend
Sleep
- Output:
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 1000th term of a[n] is 109 The 10000th term of a[n] is 101 The 100000th term of a[n] is 1693 The 1000000th term of a[n] is 1202057 The 10000000th term of a[n] is 1202057
Phix
Quite possibly flawed, but it does finish in the blink of an eye.
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)
- Output:
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"
Python
This takes under 12 minutes to run on my system (Core i5).
#!/usr/bin/python
from sympy import primefactors
res = [1]
count = 1
i = 2
lim1 = 100
lim2 = 1000
max = 1e7
while count < max:
cubeFree = False
factors = primefactors(i)
if len(factors) < 3:
cubeFree = True
else:
cubeFree = True
for j in range(2, len(factors)):
if factors[j-2] == factors[j-1] and factors[j-1] == factors[j]:
cubeFree = False
break
if cubeFree:
if count < lim1:
res.append(factors[-1])
count += 1
if count == lim1:
print("First {} terms of a[n]:".format(lim1))
for k in range(0, len(res), 10):
print(" ".join(map(str, res[k:k+10])))
print("")
elif count == lim2:
print("The {} term of a[n] is {}".format(count, factors[-1]))
lim2 *= 10
i += 1
- Output:
Similar to FreeBASIC entry.
RPL
I confirm it does take a while for an interpreted language like RPL. Getting the 100,000th term is likely to be a question of hours, even on an emulator.
« 1 { } DUP2 + → n res2 res1
« 2
WHILE 'n' INCR 10000 ≤ REPEAT
WHILE DUP FACTORS DUP 1 « 3 < NSUB 2 MOD NOT OR » DOSUBS ΠLIST NOT
REPEAT DROP 1 + END
HEAD
CASE
n 100 ≤ THEN 'res1' OVER STO+ END
n LOG FP NOT THEN 'res2' OVER STO+ END
END
DROP 1 +
END
DROP res1 res2
» » 'TASK' STO
- Output:
2: { 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 }
1: { 109 101 }
Wren
This takes just under 4 minutes to run on my system (Core i7).
Although it looks odd that the 1 millionth and 10 millionth terms are the same (and I have no independent source to check against), it would appear that the 1 millionth cubefree number is 1,202,057 which is prime and the 10 millionth such number is coincidentally 10 times this.
import "./math" for Int
import "./fmt" for Fmt
var res = [1]
var count = 1
var i = 2
var lim1 = 100
var lim2 = 1000
var max = 1e7
while (count < max) {
var cubeFree = false
var factors = Int.primeFactors(i)
if (factors.count < 3) {
cubeFree = true
} else {
cubeFree = true
for (i in 2...factors.count) {
if (factors[i-2] == factors[i-1] && factors[i-1] == factors[i]) {
cubeFree = false
break
}
}
}
if (cubeFree) {
if (count < lim1) res.add(factors[-1])
count = count + 1
if (count == lim1) {
System.print("First %(lim1) terms of a[n]:")
Fmt.tprint("$3d", res, 10)
} else if (count == lim2) {
Fmt.print("\nThe $,r term of a[n] is $,d.", count, factors[-1])
lim2 = lim2 * 10
}
}
i = i + 1
}
- Output:
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.