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
Dart
import 'dart:math';
void main() {
List<int> res = [1];
int count = 1;
int i = 2;
int lim1 = 100;
int lim2 = 1000;
double max = 1e7;
var t0 = DateTime.now();
while (count < max) {
bool cubeFree = false;
List<int> factors = primeFactors(i);
if (factors.length < 3) {
cubeFree = true;
} else {
cubeFree = true;
for (int i = 2; i < factors.length; i++) {
if (factors[i - 2] == factors[i - 1] && factors[i - 1] == factors[i]) {
cubeFree = false;
break;
}
}
}
if (cubeFree) {
if (count < lim1) res.add(factors.last);
count += 1;
if (count == lim1) {
print("First $lim1 terms of a[n]:");
print(res.take(lim1).join(', '));
print("");
} else if (count == lim2) {
print("The $count term of a[n] is ${factors.last}");
lim2 *= 10;
}
}
i += 1;
}
print("${DateTime.now().difference(t0).inSeconds} sec.");
}
List<int> primeFactors(int n) {
List<int> factors = [];
while (n % 2 == 0) {
factors.add(2);
n ~/= 2;
}
for (int i = 3; i <= sqrt(n); i += 2) {
while (n % i == 0) {
factors.add(i);
n ~/= i;
}
}
if (n > 2) {
factors.add(n);
}
return factors;
}
- 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 1000 term of a[n] is 109 The 10000 term of a[n] is 101 The 100000 term of a[n] is 1693 The 1000000 term of a[n] is 1202057 The 10000000 term of a[n] is 1202057
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
Pascal
Free Pascal
Uses factors of integer
program CubeFree;
// gets factors of consecutive integers fast
// limited to 1.2e11
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
;
//######################################################################
//prime decomposition
const
//HCN(86) > 1.2E11 = 128,501,493,120 count of divs = 4096 7 3 1 1 1 1 1 1 1
HCN_DivCnt = 4096;
type
tItem = Uint64;
tDivisors = array [0..HCN_DivCnt] of tItem;
tpDivisor = pUint64;
const
//used odd size for test only
SizePrDeFe = 32768;//*72 <= 64kb level I or 2 Mb ~ level 2 cache
type
tdigits = array [0..31] of Uint32;
//the first number with 11 different prime factors =
//2*3*5*7*11*13*17*19*23*29*31 = 2E11
//56 byte
tprimeFac = packed record
pfSumOfDivs,
pfRemain : Uint64;
pfDivCnt : Uint32;
pfMaxIdx : Uint32;
pfpotPrimIdx : array[0..9] of word;
pfpotMax : array[0..11] of byte;
end;
tpPrimeFac = ^tprimeFac;
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
tPrimes = array[0..65535] of Uint32;
var
{$ALIGN 8}
SmallPrimes: tPrimes;
{$ALIGN 32}
PrimeDecompField :tPrimeDecompField;
pdfIDX,pdfOfs: NativeInt;
function Numb2USA(n:Uint64):Ansistring;
const
//extend s by the count of comma to be inserted
deltaLength : array[0..24] of byte =
(0,0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7);
var
pI :pChar;
i,j : NativeInt;
Begin
str(n,result);
i := length(result);
//extend s by the count of comma to be inserted
// j := i+ (i-1) div 3;
j := i+deltaLength[i];
if i<> j then
Begin
setlength(result,j);
pI := @result[1];
dec(pI);
while i > 3 do
Begin
//copy 3 digits
pI[j] := pI[i];
pI[j-1] := pI[i-1];
pI[j-2] := pI[i-2];
// insert comma
pI[j-3] := ',';
dec(i,3);
dec(j,4);
end;
end;
end;
procedure InitSmallPrimes;
//get primes. #0..65535.Sieving only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
var
pr : array[0..MAXLIMIT] of byte;
p,j,d,flipflop :NativeUInt;
Begin
SmallPrimes[0] := 2;
fillchar(pr[0],SizeOf(pr),#0);
p := 0;
repeat
repeat
p +=1
until pr[p]= 0;
j := (p+1)*p*2;
if j>MAXLIMIT then
BREAK;
d := 2*p+1;
repeat
pr[j] := 1;
j += d;
until j>MAXLIMIT;
until false;
SmallPrimes[1] := 3;
SmallPrimes[2] := 5;
j := 3;
d := 7;
flipflop := (2+1)-1;//7+2*2,11+2*1,13,17,19,23
p := 3;
repeat
if pr[p] = 0 then
begin
SmallPrimes[j] := d;
inc(j);
end;
d += 2*flipflop;
p+=flipflop;
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
end;
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
//n must be multiple of base aka n mod base must be 0
var
q,r: Uint64;
i : NativeInt;
Begin
fillchar(dgt,SizeOf(dgt),#0);
i := 0;
n := n div base;
result := 0;
repeat
r := n;
q := n div base;
r -= q*base;
n := q;
dgt[i] := r;
inc(i);
until (q = 0);
//searching lowest pot in base
result := 0;
while (result<i) AND (dgt[result] = 0) do
inc(result);
inc(result);
end;
function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt;
var
q :NativeInt;
Begin
result := 0;
q := dgt[result]+1;
if q = base then
repeat
dgt[result] := 0;
inc(result);
q := dgt[result]+1;
until q <> base;
dgt[result] := q;
result +=1;
end;
function SieveOneSieve(var pdf:tPrimeDecompField):boolean;
var
dgt:tDigits;
i,j,k,pr,fac,n,MaxP : Uint64;
begin
n := pdfOfs;
if n+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then
EXIT(FALSE);
//init
for i := 0 to SizePrDeFe-1 do
begin
with pdf[i] do
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfRemain := n+i;
pfMaxIdx := 0;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := 0;
end;
end;
//first factor 2. Make n+i even
i := (pdfIdx+n) AND 1;
IF (n = 0) AND (pdfIdx<2) then
i := 2;
repeat
with pdf[i] do
begin
j := BsfQWord(n+i);
pfMaxIdx := 1;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := j;
pfRemain := (n+i) shr j;
pfSumOfDivs := (Uint64(1) shl (j+1))-1;
pfDivCnt := j+1;
end;
i += 2;
until i >=SizePrDeFe;
//i now index in SmallPrimes
i := 0;
maxP := trunc(sqrt(n+SizePrDeFe))+1;
repeat
//search next prime that is in bounds of sieve
if n = 0 then
begin
repeat
inc(i);
pr := SmallPrimes[i];
k := pr-n MOD pr;
if k < SizePrDeFe then
break;
until pr > MaxP;
end
else
begin
repeat
inc(i);
pr := SmallPrimes[i];
k := pr-n MOD pr;
if (k = pr) AND (n>0) then
k:= 0;
if k < SizePrDeFe then
break;
until pr > MaxP;
end;
//no need to use higher primes
if pr*pr > n+SizePrDeFe then
BREAK;
//j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
repeat
with pdf[k] do
Begin
pfpotPrimIdx[pfMaxIdx] := i;
pfpotMax[pfMaxIdx] := j;
pfDivCnt *= j+1;
fac := pr;
repeat
pfRemain := pfRemain DIV pr;
dec(j);
fac *= pr;
until j<= 0;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfMaxIdx);
k += pr;
j := IncByBaseInBase(dgt,pr);
end;
until k >= SizePrDeFe;
until false;
//correct sum of & count of divisors
for i := 0 to High(pdf) do
Begin
with pdf[i] do
begin
j := pfRemain;
if j <> 1 then
begin
pfSumOFDivs *= (j+1);
pfDivCnt *=2;
end;
end;
end;
result := true;
end;
function NextSieve:boolean;
begin
dec(pdfIDX,SizePrDeFe);
inc(pdfOfs,SizePrDeFe);
result := SieveOneSieve(PrimeDecompField);
end;
function GetNextPrimeDecomp:tpPrimeFac;
begin
if pdfIDX >= SizePrDeFe then
if Not(NextSieve) then
EXIT(NIL);
result := @PrimeDecompField[pdfIDX];
inc(pdfIDX);
end;
function Init_Sieve(n:NativeUint):boolean;
//Init Sieve pdfIdx,pdfOfs are Global
begin
pdfIdx := n MOD SizePrDeFe;
pdfOfs := n-pdfIdx;
result := SieveOneSieve(PrimeDecompField);
end;
function CheckCubeFree(pPrimeDecomp :tpPrimeFac):boolean;
var
i : NativeInt;
begin
with pPrimeDecomp^ do
begin
For i := 0 to pfMaxIdx-1 do
if pfpotMax[i]>2 then
EXIT(false);
EXIT(True)
end;
end;
var
pPrimeDecomp :tpPrimeFac;
T0:Int64;
n,cnt,lmt : Uint64;
Begin
InitSmallPrimes;
T0 := GetTickCount64;
cnt := 0;
n := 1;
Init_Sieve(n);
writeln('First 100 terms of a[n]:');
repeat
pPrimeDecomp:= GetNextPrimeDecomp;
if CheckCubeFree(pPrimeDecomp) then
begin
with pPrimeDecomp^ do
begin
if pfMaxIdx=0 then
write(pfRemain:4)
else
write(SmallPrimes[pfpotPrimIdx[pfMaxIdx-1]]:4);
end;
inc(cnt);
if cnt mod 10 = 0 then
writeln;
end;
inc(n);
until cnt >= 100;
writeln;
writeln(' Limit Number highest divisor');
lmt := 1000;
repeat
pPrimeDecomp:= GetNextPrimeDecomp;
if CheckCubeFree(pPrimeDecomp)then
begin
inc(cnt);
if cnt = lmt then
begin
with pPrimeDecomp^ do
begin
write(Numb2USA(lmt):17,Numb2USA(n):16);
if pfRemain <>1 then
write(Numb2USA(pFRemain):16)
else
write(Numb2USA(SmallPrimes[pfpotPrimIdx[pfMaxIdx-1]]):16);
writeln;
end;
lmt :=lmt*10;
end;
end;
inc(n);
until cnt >= 10*1000*1000;
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
end.
- @home:
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
Limit Number highest divisor
1,000 1,199 109
10,000 12,019 101
100,000 120,203 1,693
1,000,000 1,202,057 1,202,057
10,000,000 12,020,570 1,202,057
runtime 0.273 s
100,000,000 120,205,685 20,743
1,000,000,000 1,202,056,919 215,461
10,000,000,000 12,020,569,022 1,322,977
100,000,000,000 120,205,690,298 145,823
runtime 3637.753 s
real 60m37,756s
Phix
Quite possibly flawed, but it does at least complete to 1e9 in the blink of an eye, and 93s for 1e11.
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.
-- (it has o/c passed every test I can think of, but the logic of combinations/permutes behind it all eludes me)
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.
atom r = 0
bool xpm = true -- extend prior multiples
sequence pm = {}
for i=1 to floor(power(n,1/3)) do
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, bc = count_bits(mask)
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 bc=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 bc=1 then
if odd(bc) then
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...
atom t1 = time()+1
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
if time()>t1 then
progress("bisecting %,d..%,d...\r",{lo,hi})
t1 = time()+1
end if
end while
return mid
end function
function A370833(atom nth)
if nth=1 then return {1,1,1} end if
atom n = cube_free(nth)
sequence f = prime_powers(n)
return {nth,n,f[$][1]}
end function
atom t0 = time()
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(11,3)) do
printf(1,"The %,dth term of a[n] is %,d with highest divisor %,d.\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 1,199 with highest divisor 109. The 10,000th term of a[n] is 12,019 with highest divisor 101. 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 with highest divisor 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. 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"
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
Version 1
Simple brute force approach though skipping numbers which are multiples of 8 or 27 as these can't be cubefree.
Runtime about 3 minutes 36 seconds on my system (Core i7).
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
if (i % 8 == 0 || i % 27 == 0) 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.
Version 2
Astonishing speed-up, now taking only 0.04 seconds to reach 10 million and 38 seconds to reach 10 billion.
import "./math" for Int
import "./fmt" for Conv, Fmt
import "./iterate" for Stepped
var primes = Int.primeSieve(15 * 1e6)
var countBits = Fn.new { |n| Conv.bin(n).count { |c| c == "1" } }
var cubesBefore = Fn.new { |n|
var r = 0
var xpm = true
var pm = []
for (i in 1..n) {
var p3 = primes[i-1].pow(3)
if (p3 > n) break
var k = (n/p3).floor
for (mask in Stepped.ascend(1..2.pow(pm.count)-1)) {
var m = mask
var mi = 0
var kpm = p3
while (m != 0) {
mi = mi + 1
if (m % 2 == 1) kpm = kpm * pm[mi-1]
m = (m/2).floor
}
if (kpm > n) {
if (countBits.call(mask) == 1) {
xpm = false
pm = pm[0...mi-1]
break
}
} else {
var l = (n/kpm).floor
if (countBits.call(mask) % 2 == 1) {
k = k - l
} else {
k = k + l
}
}
}
r = r + k
if (xpm) pm.add(p3)
}
return r
}
var cubeFree = Fn.new { |nth|
var lo = nth
var hi = lo * 2
var mid
var cb
var k
while (hi - cubesBefore.call(hi) < nth) {
lo = hi
hi = lo * 2
}
while (true) {
mid = ((lo + hi)/2).floor
cb = cubesBefore.call(mid)
k = mid - cb
if (k == nth) {
while (cubesBefore.call(mid-1) != cb) {
mid = mid - 1
cb = cb - 1
}
break
} else if (k < nth) {
lo = mid
} else {
hi = mid
}
}
return mid
}
var a370833 = Fn.new { |n|
if (n == 1) return 1
var nth = cubeFree.call(n)
return Int.primeFactors(nth)[-1]
}
var start = System.clock
System.print("First 100 terms of a[n]:")
Fmt.tprint("$3d", (1..100).map { |i| a370833.call(i) }, 10)
System.print()
var n = 1000
while (n <= 1e10) {
Fmt.print("The $,r term of a[n] is $,d.", n, a370833.call(n))
n = n * 10
}
System.print("\nTook %(System.clock - start) seconds.")>
- 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. The 100,000,000th term of a[n] is 20,743. The 1,000,000,000th term of a[n] is 215,461. The 10,000,000,000th term of a[n] is 1,322,977. Took 37.888876 seconds.