Greatest prime dividing the n-th cubefree number: Difference between revisions
Greatest prime dividing the n-th cubefree number (view source)
Revision as of 16:39, 9 April 2024
, 2 months agoPromoted to ‘full’ task.
m (→{{header|Phix}}: admittedly remark) |
m (Promoted to ‘full’ task.) |
||
| (8 intermediate revisions by 4 users not shown) | |||
Line 1:
{{
;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 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}}==
Line 470 ⟶ 666:
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..192724.Sieving only odd numbers
const
MAXLIMIT = (2642246-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;
procedure InitPrmCubed(var DC:tPrmCubed);
var
i,q : Uint64;
begin
for i in tPrimeIdx do
begin
q := SmallPrimes[i];
q *= sqr(q);
DC[i] := q;
end;
end;
function Numb2USA(n:Uint64):Ansistring;
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;
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;
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;
procedure OutNum(lmt,n:Uint64);
var
MaxPrimeFac : Uint64;
begin
MaxPrimeFac := highestDiv(lmt);
if MaxPrimeFac > sqr(SmallPrimes[high(tPrimeIdx)]) then
MaxPrimeFac := 0;
writeln(Numb2Usa(lmt):26,'|',Numb2Usa(n):26,'|',Numb2Usa(MaxPrimeFac):15);
end;
//##########################################################################
var
cnt : Int64;
procedure check(lmt:Uint64;i:integer;flip :Boolean);
var
p : Uint64;
begin
For i := i to high(tPrimeIdx) do
begin
p := PrmCubed[i];
if lmt < p then
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 GetLmtfromCnt(inCnt:Uint64):Uint64;
begin
result := trunc(Z3*inCnt);
repeat
cnt := result;
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;
//##########################################################################
var
T0,lmt:Int64;
i : integer;
Begin
InitSmallPrimes;
InitPrmCubed(PrmCubed);
For i := 1 to 100 do
Begin
lmt := GetLmtfromCnt(i);
write(highestDiv(lmt):4);
if i mod 10 = 0 then
Writeln;
end;
Writeln;
Writeln('Tested with Apéry´s Constant approximation of ',Z3:17:15);
write(' ');
writeln('Limit | cube free numbers |max prim factor');
T0 := GetTickCount64;
lmt := 1;
For i := 0 to 18 do
Begin
OutNum(GetLmtfromCnt(lmt),lmt);
lmt *= 10;
end;
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
Tested with Apéry´s Constant approximation of 1.202056903159594
Limit | cube free numbers |max prim factor| divs
1| 1| 1| 0
11| 10| 11| 1
118| 100| 59| 2
1,199| 1,000| 109| 6<
12,019| 10,000| 101| 14
120,203| 100,000| 1,693| 30
1,202,057| 1,000,000| 1,202,057| 65<
12,020,570| 10,000,000| 1,202,057| 141
120,205,685| 100,000,000| 20,743| 301
1,202,056,919| 1,000,000,000| 215,461| 645<
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>
=={{header|Phix}}==
Line 780 ⟶ 1,225:
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) {
| |||