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

Content added Content deleted
(→‎Version 2: Added a note and changes needed to use a BitArray for the sieve.)
(→‎resursive alternative: changed to search for count of cubefree numbers)
Line 553: Line 553:
</pre>
</pre>
===resursive alternative===
===resursive alternative===
only counting til limit.Using Apéry's Constant, which is a quite good estimate.
Using Apéry's Constant, which is a quite good estimate.<br>Only checking powers of 10.Not willing to test prime factors > 6,98e12-> 0
<syntaxhighlight lang="pascal">
<syntaxhighlight lang="pascal">
program CubeFree3;
program CubeFree3;
Line 568: Line 568:
const
const
//Apéry's Constant
//Apéry's Constant
Z3 = 1.20205690315959428539973816151144999;
Z3 = 1.20205690315959428539973816151144999;
RezZ3 = 0.831907372580707468683126278821530734417;
RezZ3 = 0.831907372580707468683126278821530734417;
{
{
Line 580: Line 580:
tPrimes = array[tPrimeIdx] of Uint32;
tPrimes = array[tPrimeIdx] of Uint32;
tDl3 = UInt64;
tDl3 = UInt64;
tDelCube = array[tPrimeIdx] of tDl3;
tPrmCubed = array[tPrimeIdx] of tDl3;
var
var
{$ALIGN 8}
{$ALIGN 8}
SmallPrimes: tPrimes;
SmallPrimes: tPrimes;
DelCube : tDelCube;
PrmCubed : tPrmCubed;


procedure InitSmallPrimes;
procedure InitSmallPrimes;
Line 629: Line 629:
end;
end;


procedure InitDelCube(var DC:tDelCube);
procedure InitPrmCubed(var DC:tPrmCubed);
var
var
i,q : Uint64;
i,q : Uint64;
Line 664: Line 664:
dec(j,4);
dec(j,4);
end;
end;
end;
end;

function highestDiv(n: uint64):Uint64;
//can test til 2642246^2 ~ 6,98E12
var
pr : Uint64;
i : integer;
begin
⚫
result := n;
for i in tPrimeIdx do
begin
pr := Smallprimes[i];
if pr*pr>result then
BREAK;
while (result > pr) AND (result MOD pr = 0) do
result := result DIV pr;
end;
end;
end;
end;


procedure OutNum(lmt,n,CntDivs:Uint64);
procedure OutNum(lmt,n,CntDivs:Uint64);
var
MaxPrimeFac : Uint64;
begin
begin
MaxPrimeFac := highestDiv(lmt);
⚫
writeln(Numb2Usa(lmt):26,'|',Numb2Usa(n):26,'|',Numb2Usa(CntDivs):10);
if MaxPrimeFac > sqr(SmallPrimes[high(tPrimeIdx)]) then
MaxPrimeFac := 0;
⚫
writeln(Numb2Usa(lmt):26,'|',Numb2Usa(n):26,'|',Numb2Usa(MaxPrimeFac):15,'|',Numb2Usa(CntDivs):7);
end;
end;


Line 682: Line 704:
For i := i to high(tPrimeIdx) do
For i := i to high(tPrimeIdx) do
begin
begin
p := DelCube[i];
p := PrmCubed[i];
if lmt < p then
if lmt < p then
BREAK;
BREAK;
Line 691: Line 713:
else
else
cnt += p;
cnt += p;
if p >= DelCube[i+1] then
if p >= PrmCubed[i+1] then
check(p,i+1,not(flip));
check(p,i+1,not(flip));
end;
end;
Line 699: Line 721:
var
var
limit : extended;
limit : extended;
lmt: Uint64;
lmt,dezLmt: Int64;
lastCntDivs : Uint64;
begin
begin
limit := Z3;
limit := Z3;
⚫
dezLmt := 1;
while LmtIdx>0 do
while LmtIdx>0 do
Begin
Begin
limit *= 10;
limit *= 10;
dezLmt *=10;
dec(LmtIdx);
dec(LmtIdx);
end;
end;
lmt := trunc(limit);
lmt := trunc(limit);

⚫
cnt := lmt;
repeat
⚫
CntDivs := 0;
check(lmt,0,true);
cnt := lmt;
CntDivs := 0;
check(lmt,0,true);
//new apromitaion
inc(lmt,trunc(Z3*(dezlmt-cnt)));
until cnt = dezLmt;
// OutNum(lmt,cnt,CntDivs);

//maybe lmt is not cubefree, like 1200 for cnt 1000
//1200 was the only one ....
//faster than checking for cubefree of lmt
repeat
lastCntDivs := CntDivs;
dec(lmt);
cnt := lmt;
CntDivs := 0;
check(lmt,0,true);
until cnt < dezLmt;

inc(lmt);
cnt := dezLmt;
CntDivs := lastCntDivs;

OutNum(lmt,cnt,CntDivs);
OutNum(lmt,cnt,CntDivs);
end;
end;
Line 719: Line 766:
Begin
Begin
InitSmallPrimes;
InitSmallPrimes;
InitPrmCubed(PrmCubed);
InitDelCube(DelCube);
Writeln('Tested with Apéry´s Constant approximation of ',Z3:19:15);
write(' ');
write(' ');
writeln('Limit | cube free numbers |count of divisions');
writeln('Limit | cube free numbers |max prim factor| divs');
T0 := GetTickCount64;
T0 := GetTickCount64;
For i := 0 to 18 do
For i := 0 to 18 do
Line 731: Line 779:
{{out|@home}}
{{out|@home}}
<pre>
<pre>

Limit | cube free numbers |count of divisions
Tested with Apéry´s Constant approximation of 1.202056903159594
1| 1| 0
12| 11| 1
Limit | cube free numbers |max prim factor| divs
120| 101| 2
1| 1| 1| 0
1,202| 1,002| 6
11| 10| 11| 1
12,020| 10,001| 14
118| 100| 59| 2
120,205| 100,001| 30
1,199| 1,000| 109| 6
1,202,056| 999,999| 65
12,019| 10,000| 101| 14
12,020,569| 9,999,999| 141
120,203| 100,000| 1,693| 30
120,205,690| 100,000,004| 301
1,202,057| 1,000,000| 1,202,057| 65
1,202,056,903| 999,999,986| 645
12,020,570| 10,000,000| 1,202,057| 141
12,020,569,031| 10,000,000,008| 1,392
120,205,685| 100,000,000| 20,743| 301
120,205,690,315| 100,000,000,014| 3,003
1,202,056,919| 1,000,000,000| 215,461| 645
1,202,056,903,159| 1,000,000,000,020| 6,465
12,020,569,022| 10,000,000,000| 1,322,977| 1,392
12,020,569,031,595| 9,999,999,999,963| 13,924
120,205,690,298| 100,000,000,000| 145,823| 3,003
120,205,690,315,959| 100,000,000,000,027| 30,006
1,202,056,903,137| 1,000,000,000,000|400,685,634,379| 6,465
1,202,056,903,159,594| 1,000,000,000,000,087| 64,643
12,020,569,031,641| 10,000,000,000,000| 1,498,751| 13,924
12,020,569,031,595,942| 9,999,999,999,999,948| 139,261
120,205,690,315,927| 100,000,000,000,000| 57,349| 30,006
120,205,690,315,959,428| 100,000,000,000,000,094| 300,023
1,202,056,903,159,489| 1,000,000,000,000,000| 74,509,198,733| 64,643
1,202,056,903,159,594,285| 1,000,000,000,000,000,317| 646,394
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
⚫
runtime 0.002 s
1,202,056,903,159,593,905| 1,000,000,000,000,000,000| 89,387|646,394
⚫
runtime 0.013 s real 0m0,019
// as an experiment with factor of 1 for changing approximation.
Tested with Apéry´s Constant approximation of 1.000000000000000
runtime 0.065 s real 0m0,071s
</pre>
</pre>


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