Rosetta Code

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

Page Discussion

Greatest prime dividing the n-th cubefree number (view source)

Revision as of 16:39, 9 April 2024

17,527 bytes added ,  1 month ago
m
Promoted to ‘full’ task.
(→‎{{header|Wren}}: Added a simple optimization & rewrote preamble.)
m (Promoted to ‘full’ task.)
(20 intermediate revisions by 7 users not shown)
Line 1:
{{draft task}}
;Definitions
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.
Line 24:
;Reference
* [https://oeis.org/A370833 OEIS sequence: A370833: a(n) is the greatest prime dividing the n-th cubefree number, for n >= 2; a(1)=1.]
* [https://cp-algorithms.com/combinatorics/inclusion-exclusion.html#:~:text=The%20inclusion%2Dexclusion%20principle%20is,individual%20sets%20with%20their%20union The number of integers in a given interval which are multiple of at least one of the given numbers]
 
=={{header|Dart}}==
{{trans|Wren}}
Line 180:
The 1000000th term of a[n] is 1202057
The 10000000th term of a[n] is 1202057</pre>
 
=={{header|jq}}==
'''Works with jq, the C implementation of jq'''
 
'''Works with gojq, the Go implementation of jq'''
<syntaxhighlight lang="jq">
# The following may be omitted if using the C implementation of jq
def _nwise($n):
def n: if length <= $n then . else .[0:$n] , (.[$n:] | n) end;
n;
 
### Generic functions
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
 
# tabular print
def tprint($columns; $width):
reduce _nwise($columns) as $row ("";
. + ($row|map(lpad($width)) | join(" ")) + "\n" );
 
# like while/2 but emit the final term rather than the first one
def whilst(cond; update):
def _whilst:
if cond then update | (., _whilst) else empty end;
_whilst;
 
## Prime factors
 
# Emit an array of the prime factors of 'n' in order using a wheel with basis [2, 3, 5]
# e.g. 44 | primeFactors => [2,2,11]
def primeFactors:
def out($i): until (.n % $i != 0; .factors += [$i] | .n = ((.n/$i)|floor) );
if . < 2 then []
else [4, 2, 4, 2, 4, 6, 2, 6] as $inc
| { n: .,
factors: [] }
| out(2)
| out(3)
| out(5)
| .k = 7
| .i = 0
| until(.k * .k > .n;
if .n % .k == 0
then .factors += [.k]
| .n = ((.n/.k)|floor)
else .k += $inc[.i]
| .i = ((.i + 1) % 8)
end)
| if .n > 1 then .factors += [ .n ] else . end
| .factors
end;
 
### Cube-free numbers
# If cubefree then emit the largest prime factor, else emit null
def cubefree:
if . % 8 == 0 or . % 27 == 0 then false
else primeFactors as $factors
| ($factors|length) as $n
| {i: 2, cubeFree: true}
| until (.cubeFree == false or .i >= $n;
$factors[.i-2] as $f
| if $f == $factors[.i-1] and $f == $factors[.i]
then .cubeFree = false
else .i += 1
end)
| if .cubeFree then $factors[-1] else null end
end;
 
## The tasks
{ res: [1], # by convention
count: 1, # see the previous line
i: 2,
lim1: 100,
lim2: 1000,
max: 10000 }
| whilst (.count <= .max;
.emit = null
| (.i|cubefree) as $result
| if $result
then .count += 1
| if .count <= .lim1 then .res += [$result] end
| if .count == .lim1
then .emit = ["First \(.lim1) terms of a[n]:"]
| .emit += [.res | tprint(10; 3)]
elif .count == .lim2
then .lim2 *= 10
| .emit = ["The \(.count) term of a[n] is \($result)"]
end
end
| .i += 1
| if .i % 8 == 0 or .i % 27 == 0
then .i += 1
end
)
| select(.emit) | .emit[]
</syntaxhighlight>
{{output}}
<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 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
</pre>
 
=={{header|Julia}}==
<syntaxhighlight lang="julia">using Formatting
using Primes
using ResumableFunctions
 
const MAXINMASK = 10_000_000_000 # memory on test machine could not spare a bitmask much larger than this
 
""" return a bitmask containing at least max_wanted cubefreenumbers """
function cubefreemask(max_wanted)
size_wanted = Int(round(max_wanted * 1.21))
mask = trues(size_wanted)
p = primes(Int(floor(size_wanted^(1/3))))
for i in p
interval = i^3
for j in interval:interval:size_wanted
mask[j] = false
end
end
return mask
end
 
""" generator for cubefree numbers up to max_wanted in number """
@resumable function nextcubefree(max_wanted = MAXINMASK)
cfmask = cubefreemask(max_wanted)
@yield 1
for i in firstindex(cfmask)+1:lastindex(cfmask)
if cfmask[i]
@yield i
end
end
@warn "past end of allowable size of sequence A370833"
end
 
""" various task output with OEIS sequence A370833 """
function testA370833(toprint)
println("First 100 terms of a[n]:")
for (i, a) in enumerate(nextcubefree())
if i < 101
f = factor(a).pe # only factor the ones we want to print
highestprimefactor = isempty(f) ? 1 : f[end][begin]
print(rpad(highestprimefactor, 4), i % 10 == 0 ? "\n" : "")
elseif i ∈ toprint
highestprimefactor = (factor(a).pe)[end][begin]
println("\n The ", format(i, commas = true), "th term of a[n] is ",
format(highestprimefactor, commas = true))
end
i >= toprint[end] && break
end
end
 
testA370833(map(j -> 10^j, 3:Int(round(log10(MAXINMASK)))))
</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
 
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
</pre>
 
=={{header|Pascal}}==
==={{header|Free Pascal}}===
Uses factorssieving with cube powers of integerprimes.<BR>
Nearly linear runtime. 0.8s per Billion. Highest Prime for 1E9 is 997
<syntaxhighlight lang="pascal">
program CubeFreeCubeFree2;
// gets factors of consecutive integers fast
// limited to 1.2e11
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$OPTIMIZATION ON,ALL}
//{$CODEALIGN proc=16,loop=8} //TIO.RUN (Intel Xeon 2.3 Ghz) loves it
{$COPERATORS ON}
{$ELSE}
{$APPTYPE CONSOLE}
Line 197 ⟶ 395:
{$IFDEF WINDOWS},Windows{$ENDIF}
;
//######################################################################
//prime decomposition
const
SizeCube235 =4* (2*2*2* 3*3*3 *5*5*5);//2*27000 <= 64kb level I
//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 tPrimeIdx = Uint640..65535;
tDivisorstPrimes = array [0..HCN_DivCnttPrimeIdx] of tItemUint32;
tSv235IDx = 0..SizeCube235-1;
tpDivisor = pUint64;
tSieve235 = array[tSv235IDx] of boolean;
const
tDl3 = record
//used odd size for test only
dlPr3 : UInt64;
SizePrDeFe = 32768;//*72 <= 64kb level I or 2 Mb ~ level 2 cache
dlSivMod,
type
tdigits = array [0..31] of dlSivNum : Uint32;
end;
//the first number with 11 different prime factors =
tDelCube = array[tPrimeIdx] of tDl3;
//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;
Sieve235,
{$ALIGN 32}
Sieve : tSieve235;
PrimeDecompField :tPrimeDecompField;
DelCube : tDelCube;
pdfIDX,pdfOfs: NativeInt;
 
procedure InitSmallPrimes;
function Numb2USA(n:Uint64):Ansistring;
//get primes. #0..65535.Sieving only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
//extend s by the count of comma to be inserted
var
deltaLength : array[0..24] of byte =
pr : array[0..MAXLIMIT] 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);
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;
 
procedure Init235(var Sieve235:tSieve235);
var
i,j,k : NativeInt;
begin
fillchar(Sieve235,SizeOf(Sieve235),Ord(true));
Sieve235[0] := false;
for k in [2,3,5] do
Begin
j := k*k*k;
i := j;
while i < SizeCube235 do
begin
Sieve235[i] := false;
inc(i,j);
end;
end;
end;
 
procedure InitDelCube(var DC:tDelCube);
var
i,q,r : Uint64;
begin
for i in tPrimeIdx do
begin
q := SmallPrimes[i];
q *= sqr(q);
with DC[i] do
begin
dlPr3 := q;
r := q div SizeCube235;
dlSivNum := r;
dlSivMod := q-r*SizeCube235;
end;
end;
end;
function Numb2USA(n:Uint64):Ansistring;
var
pI :pChar;
Line 246 ⟶ 502:
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
Line 256 ⟶ 511:
Begin
//copy 3 digits
pI[j] := pI[i];pI[j-1] := pI[i-1];pI[j-2] := pI[i-2];
pI[j-1] := pI[i-1];
pI[j-2] := pI[i-2];
// insert comma
pI[j-3] := ',';
Line 266 ⟶ 519:
end;
end;
 
function highestDiv(n: uint64):Uint64;
//can test til 821641^2 ~ 6,75E11
var
pr : Uint64;
i : integer;
begin
result := n;
i := 0;
repeat
pr := Smallprimes[i];
if pr*pr>result then
BREAK;
while (result > pr) AND (result MOD pr = 0) do
result := result DIV pr;
inc(i);
until i > High(SmallPrimes);
end;
 
procedure OutNum(lmt,n:Uint64);
begin
writeln(Numb2Usa(lmt):18,Numb2Usa(n):19,Numb2Usa(highestDiv(n)):18);
end;
 
 
var
sieveNr,minIdx,maxIdx : Uint32;
procedure SieveOneSieve;
var
j : Uint64;
i : Uint32;
begin
// sieve with previous primes
Sieve := Sieve235;
For i := minIdx to MaxIdx do
with DelCube[i] do
if dlSivNum = sievenr then
begin
j := dlSivMod;
repeat
sieve[j] := false;
inc(j,dlPr3);
until j >= SizeCube235;
dlSivMod := j Mod SizeCube235;
dlSivNum += j div SizeCube235;
end;
//sieve with new primes
while DelCube[maxIdx+1].dlSivNum = sieveNr do
begin
inc(maxIdx);
with DelCube[maxIdx] do
begin
j := dlSivMod;
repeat
sieve[j] := false;
inc(j,dlPr3);
until j >= SizeCube235;
dlSivMod := j Mod SizeCube235;
dlSivNum := sieveNr + j div SizeCube235;
end;
end;
end;
 
var
T0:Int64;
cnt,lmt : Uint64;
i : integer;
 
Begin
T0 := GetTickCount64;
InitSmallPrimes;
Init235(Sieve235);
InitDelCube(DelCube);
 
sieveNr := 0;
minIdx := low(tPrimeIdx);
while Smallprimes[minIdx]<=5 do
inc(minIdx);
MaxIdx := minIdx;
while DelCube[maxIdx].dlSivNum <= sieveNr do
inc(maxIdx);
 
SieveOneSieve;
i := 1;
cnt := 0;
lmt := 100;
repeat
if sieve[i] then
begin
inc(cnt);
write(highestDiv(i):4);
if cnt mod 10 = 0 then
Writeln;
end;
inc(i);
until cnt = lmt;
Writeln;
 
cnt := 0;
lmt *=10;
repeat
For i in tSv235IDx do
begin
if sieve[i] then
begin
inc(cnt);
if cnt = lmt then
Begin
OutNum(lmt,i+sieveNr*SizeCube235);
lmt*= 10;
end;
end;
end;
inc(sieveNr);
SieveOneSieve;
until lmt > 1*1000*1000*1000;
 
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
end.
</syntaxhighlight>
{{out|@home}}
<pre>
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,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
100,000,000 120,205,685 20,743
1,000,000,000 1,202,056,919 215,461 // TIO.RUN 2.4 s
10,000,000,000 12,020,569,022 1,322,977
100,000,000,000 120,205,690,298 145,823
1,000,000,000,000 1,202,056,903,137 400,685,634,379
runtime 805.139 s
real 13m25,140s
</pre>
===resursive alternative===
Using Apéry's Constant, which is a quite good estimate.<br>Only checking powers of 10.Not willing to test prime factors > 2,642,246-> 0
<syntaxhighlight lang="pascal">
program CubeFree3;
{$IFDEF FPC}
{$MODE DELPHI}{$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$CODEALIGN proc=16,loop=8} //TIO.RUN loves loop=8
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
;
const
//Apéry's Constant
Z3 : extended = 1.20205690315959428539973816151144999;
RezZ3 = 0.831907372580707468683126278821530734417;
 
type
tPrimeIdx = 0..192724;// primes til 2,642,246 ^3 ~ 2^64-1= High(Uint64)
tPrimes = array[tPrimeIdx] of Uint32;
tDl3 = UInt64;
tPrmCubed = array[tPrimeIdx] of tDl3;
 
var
SmallPrimes: tPrimes;
{$ALIGN 32}
PrmCubed : tPrmCubed;
 
procedure InitSmallPrimes;
//get primes. #0..65535192724.Sieving only odd numbers
const
MAXLIMIT = (8216412642246-1) shr 1;
var
pr : array[0..MAXLIMIT] of byte;
Line 310 ⟶ 739:
end;
 
procedure InitPrmCubed(var DC:tPrmCubed);
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
//n must be multiple of base aka n mod base must be 0
var
qi,rq : Uint64;
begin
i : NativeInt;
for i in tPrimeIdx do
Begin
begin
fillchar(dgt,SizeOf(dgt),#0);
i q := 0SmallPrimes[i];
n := nq div*= basesqr(q);
result DC[i] := 0q;
repeatend;
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 IncByBaseInBaseNumb2USA(var dgtn:tDigits;base:NativeIntUint64):NativeIntAnsistring;
var
qpI :NativeIntpChar;
i,j : NativeInt;
Begin
str(n,result := 0);
qi := dgt[length(result]+1);
//extend s by the count of comma to be inserted
if q = base then
j := i+ (i-1) div 3;
repeat
if i<> j then
dgt[result] := 0;
Begin
inc(result);
q := dgt[setlength(result]+1,j);
untilpI q:= <> base@result[1];
dec(pI);
dgt[result] := q;
while i > 3 do
result +=1;
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;
 
function SieveOneSievehighestDiv(var pdfn:tPrimeDecompField uint64):booleanUint64;
//can test til 2642246^2 ~ 6,98E12
var
dgtpr :tDigits Uint64;
i,j,k,pr,fac,n,MaxP : Uint64integer;
begin
nresult := pdfOfsn;
for i in tPrimeIdx do
if n+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then
EXIT(FALSE);
//init
for i := 0 to SizePrDeFe-1 do
begin
withpr pdf:= Smallprimes[i] do;
if pr*pr>result then
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;
while (result > pr) AND (result MOD pr = 0) do
 
//j is powerresult of:= primeresult DIV pr;
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;
 
procedure OutNum(lmt,n:Uint64);
function NextSieve:boolean;
var
MaxPrimeFac : Uint64;
begin
MaxPrimeFac := highestDiv(lmt);
dec(pdfIDX,SizePrDeFe);
if MaxPrimeFac > sqr(SmallPrimes[high(tPrimeIdx)]) then
inc(pdfOfs,SizePrDeFe);
MaxPrimeFac := 0;
result := SieveOneSieve(PrimeDecompField);
writeln(Numb2Usa(lmt):26,'|',Numb2Usa(n):26,'|',Numb2Usa(MaxPrimeFac):15);
end;
//##########################################################################
var
cnt : Int64;
 
procedure check(lmt:Uint64;i:integer;flip :Boolean);
function GetNextPrimeDecomp:tpPrimeFac;
var
p : Uint64;
begin
For i := i to high(tPrimeIdx) do
if pdfIDX >= SizePrDeFe then
begin
if Not(NextSieve) then
p := EXIT(NIL)PrmCubed[i];
if lmt < p then
result := @PrimeDecompField[pdfIDX];
inc(pdfIDX) BREAK;
p := lmt DIV p;
if flip then
cnt -= p
else
cnt += p;
if p >= PrmCubed[i+1] then
check(p,i+1,not(flip));
end;
end;
 
function Init_SieveGetLmtfromCnt(ninCnt:NativeUintUint64):booleanUint64;
//Init Sieve pdfIdx,pdfOfs are Global
begin
pdfIdxresult := n MOD SizePrDeFetrunc(Z3*inCnt);
repeat
pdfOfs := n-pdfIdx;
cnt := result;
result := SieveOneSieve(PrimeDecompField);
check(result,0,true);
//new approximation
inc(result,trunc(Z3*(inCnt-Cnt)));
until cnt = inCnt;
//maybe lmt is not cubefree, like 1200 for cnt 1000
//faster than checking for cubefree of lmt for big lmt
repeat
dec(result);
cnt := result;
check(result,0,true);
until cnt < inCnt;
inc(result);
end;
//##########################################################################
 
function CheckCubeFree(pPrimeDecomp :tpPrimeFac):boolean;
var
T0,lmt:Int64;
i : NativeInt;
i : integer;
begin
Begin
with pPrimeDecomp^ do
InitSmallPrimes;
begin
InitPrmCubed(PrmCubed);
For i := 0 to pfMaxIdx-1 do
 
if pfpotMax[i]>2 then
For i := 1 to 100 EXIT(false);do
Begin
EXIT(True)
lmt := GetLmtfromCnt(i);
write(highestDiv(lmt):4);
if i mod 10 = 0 then
Writeln;
end;
Writeln;
end;
 
Writeln('Tested with Apéry´s Constant approximation of ',Z3:17:15);
var
write(' ');
pPrimeDecomp :tpPrimeFac;
writeln('Limit | cube free numbers |max prim factor');
T0:Int64;
n,cnt,lmt : Uint64;
Begin
InitSmallPrimes;
T0 := GetTickCount64;
cntlmt := 01;
nFor i := 1;0 to 18 do
Begin
Init_Sieve(n);
OutNum(GetLmtfromCnt(lmt),lmt);
writeln('First 100 terms of a[n]:');
lmt *= 10;
repeat
end;
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.
end.</syntaxhighlight>
{{out|@home}}
<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
Line 567 ⟶ 889:
107 109 11 37 113 19 23 29 13 59
 
Tested with Apéry´s Constant approximation of 1.202056903159594
Limit Number highest divisor
1,000 Limit | 1,199 cube free numbers |max prim factor| 109 divs
10,000 12,019 1| 101 1| 1| 0
100,000 120,203 11| 10| 11| 1,693
1,000,000 1,202,057 1,202,057 118| 100| 59| 2
10,000,000 12 1,020,570199| 1,202,057000| 109| 6<
12,019| 10,000| 101| 14
runtime 0.273 s
100,000,000 120,205203| 100,685000| 201,743693| 30
1,000202,000,000057| 1,202000,056,919000| 1,202,057| 215,461 65<
10,000,000,000 12,020,569,022570| 10,000,000| 1,322202,977057| 141
100,000,000,000 120,205,690685| 100,298000,000| 14520,823743| 301
1,202,056,919| 1,000,000,000| 215,461| 645<
runtime 3637.753 s
12,020,569,022| 10,000,000,000| 1,322,977| 1,392
120,205,690,298| 100,000,000,000| 145,823| 3,003
1,202,056,903,137| 1,000,000,000,000|400,685,634,379| 6,465<
12,020,569,031,641| 10,000,000,000,000| 1,498,751| 13,924
120,205,690,315,927| 100,000,000,000,000| 57,349| 30,006
1,202,056,903,159,489| 1,000,000,000,000,000| 74,509,198,733| 64,643<
12,020,569,031,596,003| 10,000,000,000,000,000| 0|139,261
120,205,690,315,959,316| 100,000,000,000,000,000| 0|300,023
1,202,056,903,159,593,905| 1,000,000,000,000,000,000| 89,387|646,394<
runtime 0.008 s //best value.
real 0m0,013s
Tested with Apéry´s Constant approximation of 1.000000000000000
runtime 0.065 s
real 0m0,071s</pre>
 
real 60m37,756s
</pre>
=={{header|Phix}}==
Completes to 1e9 in the blink of an eye, and 93s for 1e11. Admittedly 1e12 is well out of reach, and this is far from the best way to get the full series.
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.
-- (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.
Line 595 ⟶ 925:
-- but only if cubes_before(n-1)==cubes_before(n),
-- otherwise cubes_before(cubicate) isn't really very meaningful.
integeratom r = 0
bool xpm = true -- extend prior multiples
sequence pm = {}
for i=1 to n do -- (or eg floor(cbrtpower(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
Line 613 ⟶ 943:
end while
if kpm>n then
if count_bits(mask)bc=1 then
xpm = false
pm = pm[1..mi-1]
Line 619 ⟶ 949:
end if
else
-- subtract?account for already counted multiples..
integer l = floor(n/kpm)
-- -pairs +triples -quads ... as per link above
-- DEV insufficient kludge... (probably)
-- if count_bitsodd(maskbc)=1 then
k -= l
if odd(count_bits(mask)) then -- match Wren
k -= lelse
k += l
else
kend += lif
end if
end if
end for
Line 646 ⟶ 975:
end while
-- bisect until we have a possible...
atom t1 = time()+1
while true do
mid = floor((lo+hi)/2)
Line 661 ⟶ 991:
else
hi = mid
end if
if time()>t1 then
progress("bisecting %,d..%,d (diff %,d)...\r",{lo,hi,hi-lo})
t1 = time()+1
end if
end while
Line 666 ⟶ 1,000:
end function
 
function A370833(atom nnth)
if nnth=1 then return {1,1,1} end if
atom nthn = cube_free(nnth)
sequence f = prime_powers(nthn)
return {nth,n,f[$][1]}
end function
 
atom t0 = time()
printf(1,"Firstsequence 100f100 terms= of a[n]:\n%s\n",join_byvslice(apply(tagset(100),A370833),1,10," ",fmt:="%3d")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
?elapsed(time()-t0)
Line 694 ⟶ 1,029:
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.
"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>
 
Line 744 ⟶ 1,083:
{{out}}
<pre>Similar to FreeBASIC entry.</pre>
 
=={{header|Raku}}==
 
<syntaxhighlight lang="raku">
use Prime::Factor;
 
sub max_factor_if_cubefree ($i) {
my @f = prime-factors($i);
return @f.tail if @f.elems < 3
or @f.Bag.values.all < 3;
}
 
constant @Aₙ = lazy flat 1, map &max_factor_if_cubefree, 2..*;
 
say 'The first terms of A370833 are:';
say .fmt('%3d') for @Aₙ.head(100).batch(10);
 
say '';
 
for 10 X** (3..6) -> $k {
printf "The %8dth term of A370833 is %7d\n", $k, @Aₙ[$k-1];
}
</syntaxhighlight>
 
{{out}}
<pre>
The first terms of A370833 are:
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 A370833 is 109
The 10000th term of A370833 is 101
The 100000th term of A370833 is 1693
The 1000000th term of A370833 is 1202057</pre>
 
=={{header|RPL}}==
Line 771 ⟶ 1,152:
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
 
===Version 1===
Simple brute force approach though skipping numbers which are multiples of 8 or 27 as these can't be cubefree.
 
Line 835 ⟶ 1,218:
 
The 10,000,000th term of a[n] is 1,202,057.
</pre>
 
===Version 2===
This uses a simple sieve to find cubefree numbers up to a given limit which means we only now need to factorize the numbers of interest to find the greatest prime factor.
 
The 10 millionth term is now found in 0.85 seconds and the 1 billionth in about 94 seconds.
 
However, a lot of memory is needed for the sieve since all values in Wren (including bools) need 8 bytes of storage each.
 
We could use only 1/32nd as much memory by importing the BitArray class from [[:Category:Wren-array|Wren-array]] (see program comments for changes needed). However, unfortunately this is much slower to index than a normal List of booleans and the 10 millionth term would now take just over 2 seconds to find and the 1 billionth just under 4 minutes.
<syntaxhighlight lang="wren">import "./math" for Int
//import "./array" for BitArray
import "./fmt" for Fmt
 
var cubeFreeSieve = Fn.new { |n|
var cubeFree = List.filled(n+1, true) // or BitArray.new(n+1, true)
var primes = Int.primeSieve(n.cbrt.ceil)
for (p in primes) {
var p3 = p * p * p
var k = 1
while (p3 * k <= n) {
cubeFree[p3 * k] = false
k = k + 1
}
}
return cubeFree
}
 
var al = [1]
var count = 1
var i = 2
var lim1 = 100
var lim2 = 1000
var max = 1e9
var cubeFree = cubeFreeSieve.call(max * 1.25)
while (count < max) {
if (cubeFree[i]) {
count = count + 1
if (count <= lim1) {
var factors = Int.primeFactors(i)
al.add(factors[-1])
if (count == lim1) {
System.print("First %(lim1) terms of a[n]:")
Fmt.tprint("$3d", al, 10)
}
} else if (count == lim2) {
var factors = Int.primeFactors(i)
Fmt.print("\nThe $,r term of a[n] is $,d.", count, factors[-1])
lim2 = lim2 * 10
}
}
i = i + 1
}</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.
 
The 100,000,000th term of a[n] is 20,743.
 
The 1,000,000,000th term of a[n] is 215,461.
</pre>
 
===Version 3===
{{trans|Phix}}
{{libheader|Wren-iterate}}
Even greater speed-up, now taking only 0.04 seconds to reach 10 million and 36.5 seconds to reach 10 billion.
<syntaxhighlight lang="wren">import "./math" for Int
import "./fmt" for Conv, Fmt
import "./iterate" for Stepped
 
var start = System.clock
 
var primes = Int.primeSieve(3000)
 
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.cbrt.floor) {
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 bc = countBits.call(mask)
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 (bc == 1) {
xpm = false
pm = pm[0...mi-1]
break
}
} else {
var l = (n/kpm).floor
if (bc % 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, 1]
var nth = cubeFree.call(n)
return [Int.primeFactors(nth)[-1], nth]
}
 
System.print("First 100 terms of a[n]:")
Fmt.tprint("$3d", (1..100).map { |i| a370833.call(i)[0] }, 10)
System.print()
var n = 1000
while (n <= 1e10) {
var res = a370833.call(n)
Fmt.print("The $,r term of a[n] is $,d for cubefree number $,d.", n, res[0], res[1])
n = n * 10
}
System.print("\nTook %(System.clock - start) seconds.")</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 for cubefree number 1,199.
The 10,000th term of a[n] is 101 for cubefree number 12,019.
The 100,000th term of a[n] is 1,693 for cubefree number 120,203.
The 1,000,000th term of a[n] is 1,202,057 for cubefree number 1,202,057.
The 10,000,000th term of a[n] is 1,202,057 for cubefree number 12,020,570.
The 100,000,000th term of a[n] is 20,743 for cubefree number 120,205,685.
The 1,000,000,000th term of a[n] is 215,461 for cubefree number 1,202,056,919.
The 10,000,000,000th term of a[n] is 1,322,977 for cubefree number 12,020,569,022.
 
Took 36.458647 seconds.
</pre>
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/361294...363410"