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 13:50, 18 May 2024

15,446 bytes added ,  1 month ago
Added XPL0 example.
(Added a ref to "The number of integers in a given interval which are multiple of at least one of the given numbers")
(Added XPL0 example.)
 
(12 intermediate revisions by 6 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 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}}==
Quite possibly flawed, but it does at least completeCompletes 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.
<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 508 ⟶ 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 odd(bc=1) then
k -= l
if odd(bc) then
k -= lelse
k += l
else
kend += lif
end if
end if
end for
Line 553 ⟶ 993:
end if
if time()>t1 then
progress("bisecting %,d..%,d (diff %,d)...\r",{lo,hi,hi-lo})
t1 = time()+1
end if
Line 643 ⟶ 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 743 ⟶ 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) {
Line 939 ⟶ 1,424:
 
Took 36.458647 seconds.
</pre>
 
=={{header|XPL0}}==
Simple brute force takes 43 seconds on Raspberry Pi 4.
<syntaxhighlight lang "XPL0">func MaxFactor(N); \Return maximum factor of N that's cube-free
int N, Div, Count, Max;
[Max:= N; Div:= 2;
while N >= Div*Div do
[Count:= 0;
while rem(N/Div) = 0 do
[Count:= Count+1;
if Count = 3 then return 0;
Max:= Div;
N:= N/Div;
];
Div:= Div+1;
];
if N > 1 then Max:= N; \prime
return Max;
];
 
int I, N, Pow10, Fac;
[Format(4, 0);
I:= 1; N:= 0; Pow10:= 1000;
loop [Fac:= MaxFactor(I);
if Fac # 0 then
[N:= N+1;
if N <= 100 then
[RlOut(0, float(Fac));
if rem(N/10) = 0 then CrLf(0);
]
else if N = Pow10 then
[IntOut(0, Pow10); Text(0, "th term of a[n] is ");
IntOut(0, Fac); CrLf(0);
if Pow10 = 10_000_000 then quit;
Pow10:= Pow10*10;
];
];
I:= I+1;
];
]</syntaxhighlight>
{{out}}
<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
1000th term of a[n] is 109
10000th term of a[n] is 101
100000th term of a[n] is 1693
1000000th term of a[n] is 1202057
10000000th term of a[n] is 1202057
</pre>
Boreal
297

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/361335...364545"