Radical of an integer
Radical of an integer 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.
- Definition
The radical of a positive integer is defined as the product of its distinct prime factors.
Although the integer 1 has no prime factors, its radical is regarded as 1 by convention.
- Example
The radical of 504 = 2³ x 3² x 7 is: 2 x 3 x 7 = 42.
- Task
1. Find and show on this page the radicals of the first 50 positive integers.
2. Find and show the radicals of the integers: 99999, 499999 and 999999.
3. By considering their radicals, show the distribution of the first one million positive integers by numbers of distinct prime factors (hint: the maximum number of distinct factors is 7).
- Bonus
By (preferably) using an independent method, calculate the number of primes and the number of powers of primes less than or equal to one million and hence check that your answer in 3. above for numbers with one distinct prime is correct.
- Related tasks
- References
- Wikipedia article Radical of an integer
- OEIS sequence A007947: Largest square free number dividing n
FreeBASIC[edit]
'#include "isprime.bas"
Function Radical(n As Integer) As Integer 'Return radical of n
Dim As Integer d = 2, d0 = 0, p = 1
While n >= d*d
While n Mod d = 0
If d <> d0 Then
p *= d
d0 = d
End If
n /= d
Wend
d += 1
Wend
If d <> d0 Then p *= n
Return p
End Function
Function DistinctFactors(n As Integer) As Integer 'Return count of distinct factors of n
Dim As Integer d = 2, d0 = 0, c = 0
While n >= d*d
While n Mod d = 0
If d <> d0 Then
c += 1
d0 = d
End If
n = n/d
Wend
d += 1
Wend
If d <> d0 And n <> 1 Then c += 1
Return c
End Function
Dim As Integer n, c, count(0 To 9), pcount, ppcount, p2, p
Print "The radicals for the first 50 positive integers are:"
For n = 1 To 50
Print Using "####"; Radical(n);
If n Mod 10 = 0 Then Print
Next
Print Using !"\nRadical for ###,###: ###,###"; 99999; Radical(99999)
Print Using "Radical for ###,###: ###,###"; 499999; Radical(499999)
Print Using "Radical for ###,###: ###,###"; 999999; Radical(999999)
For n = 0 To 9
count(n) = 0
Next
For n = 1 To 1e6
c = DistinctFactors(n)
count(c) += 1
Next
Print !"\nBreakdown of numbers of distinct prime factors \nfor positive integers from 1 to 1,000,000:"
c = 0
For n = 0 To 9
If count(n) > 0 Then Print Using " #: ###,###"; n; count(n)
c += count(n)
Next
Print Using !" ---------\n #,###,###"; c
Print !"\nor graphically:\n"; String(50, "-")
For n = 0 To 9
If count(n) > 0 Then Print n; " "; String(count(n)/1e4, Chr(177)); " "; count(n)
c += count(n)
Next
Print String(50, "-")
'Bonus (algorithm from Wren):
pcount = 0
For n = 1 To 1e6
If isPrime(n) Then pcount += 1
Next
Print !"\nFor primes or powers (>1) thereof <= 1,000,000:"
Print Using " Number of primes: ##,###"; pcount
ppcount = 0
For p = 1 To Sqr(1e6)
If isPrime(p) Then
p2 = p
Do
p2 *= p
If p2 > 1e6 Then Exit Do
ppcount += 1
Loop
End If
Next
Print Using " Number of powers: ##,###"; ppcount
Print Using "Add 1 for number 1: ##,###"; 1
If isPrime(n) Then pcount += 1
Print Spc(19); "-------"
Print Spc(13); Using !"Total: ##,###"; pcount + ppcount + 1
Sleep
- Output:
The radicals for the first 50 positive integers are:
1 2 3 2 5 6 7 2 3 10
11 6 13 14 15 2 17 6 19 10
21 22 23 6 5 26 3 14 29 30
31 2 33 34 35 6 37 38 39 10
41 42 43 22 15 46 47 6 7 10
Radical for 99,999: 33,333
Radical for 499,999: 3,937
Radical for 999,999: 111,111
Breakdown of numbers of distinct prime factors
for positive integers from 1 to 1,000,000:
0: 1
1: 78,734
2: 288,726
3: 379,720
4: 208,034
5: 42,492
6: 2,285
7: 8
---------
1,000,000
or graphically:
--------------------------------------------------
0 1
1 ▒▒▒▒▒▒▒▒ 78734
2 ▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ 288726
3 ▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ 379720
4 ▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒ 208034
5 ▒▒▒▒ 42492
6 2285
7 8
--------------------------------------------------
For primes or powers (>1) thereof <= 1,000,000:
Number of primes: 78,498
Number of powers: 236
Add 1 for number 1: 1
-------
Total: 78,735
J[edit]
~.&.q: 1+i.5 10 NB. radicals of first 50 positive integers
1 2 3 2 5 6 7 2 3 10
11 6 13 14 15 2 17 6 19 10
21 22 23 6 5 26 3 14 29 30
31 2 33 34 35 6 37 38 39 10
41 42 43 22 15 46 47 6 7 10
~.&.q: 99999 499999 999999 NB. radicals of these three...
33333 3937 111111
(~.,.#/.~) 1>.#@~.@q: 1+i.1e6 NB. distribution of number of prime factors of first million positive integers
1 78735
2 288726
3 379720
4 208034
5 42492
6 2285
7 8
p:inv 1e6 NB. number of primes not exceeding 1 million
78498
+/_1+<.(i.&.(p:inv) 1000)^.1e6 NB. count of prime powers (square or above) up to 1 million
236
78498+236+1 NB. and we "claimed" that 1 had a prime factor
78735
jq[edit]
Adapted from Wren
Also works with gojq, the Go implementations of jq, except that gojq is likely to run out of memory before completing part2.
# Utility functions
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
def prod(s): reduce s as $x (1; . * $x);
def sum(s): reduce s as $x (0; . + $x);
def uniq(s):
foreach s as $x (null;
if . and $x == .[0] then .[1] = false
else [$x, true]
end;
if .[1] then .[0] else empty end);
# Prime number functions
# Returns the prime factors of . in order using a wheel with basis [2, 3, 5].
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;
# Input: a positive integer
# Output: an array, $a, of length .+1 such that
# $a[$i] is $i if $i is prime, and false otherwise.
def primeSieve:
# erase(i) sets .[i*j] to false for integral j > 1
def erase($i):
if .[$i] then
reduce (range(2*$i; length; $i)) as $j (.; .[$j] = false)
else .
end;
(. + 1) as $n
| (($n|sqrt) / 2) as $s
| [null, null, range(2; $n)]
| reduce (2, 1 + (2 * range(1; $s))) as $i (.; erase($i)) ;
# Number of primes up to and including .
def primeCount:
sum(primeSieve[] | select(.) | 1);
## Radicals
def task1:
{ radicals: [0],
counts: [range(0;8)|0] }
| .radicals[1] = 1
| .counts[1] = 1
| foreach range(2; 1+1e6) as $i (.;
.factors = [uniq($i|primeFactors[])]
| (.factors|length) as $fc
| .counts[$fc] += 1
| if $i <= 50 then .radicals[$i] = prod(.factors[]) else . end ;
if $i == 50
then "The radicals for the first 50 positive integers are:",
(.radicals[1:] | _nwise(10) | map(lpad(4)) | join(" ")),
""
elif $i | IN( 99999, 499999, 999999)
then "Radical for \($i|lpad(8)): \(prod(.factors[])|lpad(8))"
elif $i == 1e6
then "\nBreakdown of numbers of distinct prime factors",
"for positive integers from 1 to 1,000,000:",
(range(1; 8) as $i
| " \($i): \(.counts[$i]|lpad(8))"),
" ---------",
" \(sum(.counts[]))"
else empty
end);
def task2:
def pad: lpad(6);
(1000|primeSieve|map(select(.))) as $primes1k
| { pcount: (1e6|primeCount),
ppcount: 0 }
| reduce $primes1k[] as $p (.;
.p2 = $p
| .done = false
| until(.done;
.p2 *= $p
| if .p2 > 1e6 then .done = true
else .ppcount += 1
end ) )
| "\nFor primes or powers (>1) thereof <= 1,000,000:",
" Number of primes = \(.pcount|pad)",
" Number of powers = \(.ppcount|pad)",
" Add 1 for number 1 = \(1|pad)",
" ------",
" \( (.pcount + .ppcount + 1)|pad)" ;
task1, task2The radicals for the first 50 positive integers are:
1 2 3 2 5 6 7 2 3 10
11 6 13 14 15 2 17 6 19 10
21 22 23 6 5 26 3 14 29 30
31 2 33 34 35 6 37 38 39 10
41 42 43 22 15 46 47 6 7 10
Radical for 99999: 33333
Radical for 499999: 3937
Radical for 999999: 111111
Breakdown of numbers of distinct prime factors
for positive integers from 1 to 1,000,000:
1: 78735
2: 288726
3: 379720
4: 208034
5: 42492
6: 2285
7: 8
---------
1000000
For primes or powers (>1) thereof <= 1,000,000:
Number of primes = 78498
Number of powers = 236
Add 1 for number 1 = 1
------
78735
Julia[edit]
using Formatting, Primes, UnicodePlots
radical(n) = prod(map(first, factor(n).pe))
radicallength(n) = length(factor(n).pe)
println("The radicals for the first 50 positive integers are:")
foreach(p -> print(rpad(p[2], 4), p[1] % 10 == 0 ? "\n" : ""), enumerate(map(radical, 1:50)))
for i in [99999, 499999, 999999]
println("\nRadical for ", format(i, commas=true), ": ", format(radical(i), commas=true))
end
println("\nBreakdown of numbers of distinct prime factors for positive integers from 1 to 1,000,000:")
histogram(map(radicallength, 1:1_000_000), nbins = 8)
println("\nCheck on breakdown:")
primecount = length(primes(1_000_000)) # count of primes to 1 million
powerscount = mapreduce(p -> Int(floor(6 / log10(p)) - 1), +, primes(1000))
println("Prime count to 1 million: ", lpad(primecount, 6))
println("Prime powers less than 1 million: ", lpad(powerscount, 6))
println("Subtotal:", lpad(primecount + powerscount, 32))
println("The integer 1 has 0 prime factors: ", lpad(1, 6))
println("-"^41, "\n", "Overall total:", lpad(primecount + powerscount + 1, 27))
- Output:
The radicals for the first 50 positive integers are:
1 2 3 2 5 6 7 2 3 10
11 6 13 14 15 2 17 6 19 10
21 22 23 6 5 26 3 14 29 30
31 2 33 34 35 6 37 38 39 10
41 42 43 22 15 46 47 6 7 10
Radical for 99,999: 33,333
Radical for 499,999: 3,937
Radical for 999,999: 111,111
Breakdown of numbers of distinct prime factors for positive integers from 1 to 1,000,000:
┌ ┐
[0.0, 1.0) ┤▏ 1
[1.0, 2.0) ┤██████▌ 78 734
[2.0, 3.0) ┤███████████████████████▌ 288 726
[3.0, 4.0) ┤███████████████████████████████ 379 720
[4.0, 5.0) ┤████████████████▉ 208 034
[5.0, 6.0) ┤███▌ 42 492
[6.0, 7.0) ┤▎ 2 285
[7.0, 8.0) ┤▏ 8
└ ┘
Frequency
Check on breakdown:
Prime count to 1 million: 78498
Prime powers less than 1 million: 236
Subtotal: 78734
The integer 1 has 0 prime factors: 1
-----------------------------------------
Overall total: 78735
Nim[edit]
import std/[math, strformat, strutils]
const N = 1_000_000
### Build list of primes.
func isPrime(n: Natural): bool =
## Return true if "n" is prime.
## "n" should not be a mutiple of 2 or 3.
var k = 5
var delta = 2
while k * k <= n:
if n mod k == 0: return false
inc k, delta
delta = 6 - delta
result = true
var primes = @[2, 3]
var n = 5
var step = 2
while n <= N:
if n.isPrime:
primes.add n
inc n, step
step = 6 - step
### Build list of distinct prime factors to
### compute radical and distinct factor count.
var primeFactors: array[1..N, seq[int]]
for p in primes:
for n in countup(p, N, p):
primeFactors[n].add p
template radical(n: int): int = prod(primeFactors[n])
template factorCount(n: int): int = primeFactors[n].len
### Task ###
echo "Radical of first 50 positive integers:"
for n in 1..50:
stdout.write &"{radical(n):2}"
stdout.write if n mod 10 == 0: '\n' else: ' '
echo()
for n in [99_999, 499_999, 999_999]:
echo &"Radical for {insertSep($n):>7}: {insertSep($radical(n)):>7}"
echo()
echo "Distribution of the first one million positive"
echo "integers by numbers of distinct prime factors:"
var counts: array[0..7, int]
for n in 1..1_000_000:
inc counts[factorCount(n)]
for n, count in counts:
echo &"{n}: {insertSep($count):>7}"
echo()
### Bonus ###
echo "Number of primes and powers of primes"
echo "less than or equal to one million:"
var count = 0
const LogN = ln(N.toFloat)
for p in primes:
inc count, int(LogN / ln(p.toFloat))
echo insertSep($count)
- Output:
Radical of first 50 positive integers: 1 2 3 2 5 6 7 2 3 10 11 6 13 14 15 2 17 6 19 10 21 22 23 6 5 26 3 14 29 30 31 2 33 34 35 6 37 38 39 10 41 42 43 22 15 46 47 6 7 10 Radical for 99_999: 33_333 Radical for 499_999: 3_937 Radical for 999_999: 111_111 Distribution of the first one million positive integers by numbers of distinct prime factors: 0: 1 1: 78_734 2: 288_726 3: 379_720 4: 208_034 5: 42_492 6: 2_285 7: 8 Number of primes and powers of primes less than or equal to one million: 78_734
Pascal[edit]
Free Pascal[edit]
program Radical;
{$IFDEF FPC} {$MODE DELPHI}{$Optimization ON,ALL} {$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
//much faster would be
//https://rosettacode.org/wiki/Factors_of_an_integer#using_Prime_decomposition
const
LIMIT = 1000*1000;
DeltaMod235 : array[0..7] of Uint32 = (4, 2, 4, 2, 4, 6, 2, 6);
type
tRadical = record
number,radical,PrFacCnt: Uint64;
isPrime : boolean;
end;
function GetRadical(n: UInt32):tRadical;forward;
function CommaUint(n : Uint64):AnsiString;
//commatize only positive Integers
var
fromIdx,toIdx :Int32;
pRes : pChar;
Begin
str(n,result);
fromIdx := length(result);
toIdx := fromIdx-1;
if toIdx < 3 then
exit;
toIdx := 4*(toIdx DIV 3)+toIdx MOD 3 +1 ;
setlength(result,toIdx);
pRes := @result[1];
dec(pRes);
repeat
pRes[toIdx] := pRes[FromIdx];
pRes[toIdx-1] := pRes[FromIdx-1];
pRes[toIdx-2] := pRes[FromIdx-2];
pRes[toIdx-3] := ',';
dec(toIdx,4);
dec(FromIdx,3);
until FromIdx<=3;
while fromIdx>=1 do
Begin
pRes[toIdx] := pRes[FromIdx];
dec(toIdx);
dec(fromIdx);
end;
end;
procedure OutRadical(n: Uint32);
Begin
writeln('Radical for ',CommaUint(n):8,':',CommaUint(GetRadical(n).radical):8);
end;
function GetRadical(n: UInt32):tRadical;
var
q,divisor, rest: UInt32;
nxt : Uint32;
begin
with result do
Begin
number := n;
radical := n;
PrFacCnt := 1;
isPrime := false;
end;
if n <= 1 then
EXIT;
if n in [2,3,5,7,11,13,17,19,23,29,31] then
Begin
with result do
Begin
isprime := true;
PrFacCnt := 1;
end;
EXIT;
end;
with result do
Begin
radical := 1;
PrFacCnt := 0;
end;
rest := n;
if rest AND 1 = 0 then
begin
with result do begin radical := 2; PrFacCnt:= 1;end;
repeat
rest := rest shr 1;
until rest AND 1 <> 0;
end;
if rest < 3 then
EXIT;
q := rest DIV 3;
if rest-q*3= 0 then
begin
with result do begin radical *= 3; inc(PrFacCnt);end;
repeat
rest := q;
q := rest DIV 3;
until rest-q*3 <> 0;
end;
if rest < 5 then
EXIT;
q := rest DIV 5;
if rest-q*5= 0 then
begin
with result do begin radical *= 5;inc(PrFacCnt);end;
repeat
rest := q;
q := rest DIV 5;
until rest-q*5 <> 0;
end;
divisor := 7;
nxt := 0;
repeat;
if rest < sqr(divisor) then
BREAK;
q := rest DIV divisor;
if rest-q*divisor= 0 then
begin
with result do begin radical *= divisor; inc(PrFacCnt);end;
repeat
rest := q;
q := rest DIV divisor;
until rest-q*divisor <> 0;
end;
divisor += DeltaMod235[nxt];
nxt := (nxt+1) AND 7;
until false;
//prime ?
if rest = n then
with result do begin radical := n;PrFacCnt:=1;isPrime := true; end
else
if rest >1 then
with result do begin radical *= rest;inc(PrFacCnt);end;
end;
var
Rad:tRadical;
CntOfPrFac : array[0..9] of Uint32;
j,sum,countOfPrimes,CountOfPrimePowers: integer;
begin
writeln('The radicals for the first 50 positive integers are:');
for j := 1 to 50 do
Begin
write (GetRadical(j).radical:4);
if j mod 10 = 0 then
Writeln;
end;
writeln;
OutRadical( 99999);
OutRadical(499999);
OutRadical(999999);
writeln;
writeln('Breakdown of numbers of distinct prime factors');
writeln('for positive integers from 1 to ',CommaUint(LIMIT));
countOfPrimes:=0;
CountOfPrimePowers :=0;
For j := Low(CntOfPrFac) to High(CntOfPrFac) do
CntOfPrFac[j] := 0;
For j := 1 to LIMIT do
Begin
Rad := GetRadical(j);
with rad do
Begin
IF isPrime then
inc(countOfPrimes)
else
if (j>1)AND(PrFacCnt= 1) then
inc(CountOfPrimePowers);
end;
inc(CntOfPrFac[Rad.PrFacCnt]);
end;
sum := 0;
For j := Low(CntOfPrFac) to High(CntOfPrFac) do
if CntOfPrFac[j] > 0 then
Begin
writeln(j:3,': ',CommaUint(CntOfPrFac[j]):10);
inc(sum,CntOfPrFac[j]);
end;
writeln('sum: ',CommaUint(sum):10);
writeln;
sum := countOfPrimes+CountOfPrimePowers+1;
writeln('For primes or powers (>1) there of <= ',CommaUint(LIMIT));
Writeln(' Number of primes =',CommaUint(countOfPrimes):8);
Writeln(' Number of prime powers =',CommaUint(CountOfPrimePowers):8);
Writeln(' Add 1 for number = 1');
Writeln(' sums to =',CommaUint(sum):8);
{$IFDEF WINDOWS}readln;{$ENDIF}
end.
- @TIO.RUN:
The radicals for the first 50 positive integers are: 1 2 3 2 5 6 7 2 3 10 11 6 13 14 15 2 17 6 19 10 21 22 23 6 5 26 3 14 29 30 31 2 33 34 35 6 37 38 39 10 41 42 43 22 15 46 47 6 7 10 Radical for 99,999: 33,333 Radical for 499,999: 3,937 Radical for 999,999: 111,111 Breakdown of numbers of distinct prime factors for positive integers from 1 to 1,000,000 1: 78,735 2: 288,726 3: 379,720 4: 208,034 5: 42,492 6: 2,285 7: 8 sum: 1,000,000 For primes or powers (>1) there of <= 1,000,000 Number of primes = 78,498 Number of prime powers = 236 Add 1 for number = 1 sums to = 78,735 Real time: 0.310 s User time: 0.285 s Sys. time: 0.022 s CPU share: 99.16 %
alternative[edit]
Using modified factors of integer inserted radical.
program Radical_2;
// gets factors of consecutive integers fast limited to 1.2e11
{$IFDEF FPC}{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$ENDIF}
{$IFDEF Windows}{$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;//*+64 <= 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
//64 byte
tprimeFac = packed record
pfSumOfDivs,
pfRemain,
pfRadical : Uint64;
pfDivCnt,
pfPrimeCnt : 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;
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 OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring;
var
s: String[31];
chk,p,i: NativeInt;
Begin
str(n,s);
result := s+' :';
with pd^ do
begin
str(pfDivCnt:3,s);
result += s+' : ';
chk := 1;
For n := 0 to pfPrimeCnt-1 do
Begin
if n>0 then
result += '*';
p := SmallPrimes[pfpotPrimIdx[n]];
chk *= p;
str(p,s);
result += s;
i := pfpotMax[n];
if i >1 then
Begin
str(pfpotMax[n],s);
result += '^'+s;
repeat
chk *= p;
dec(i);
until i <= 1;
end;
end;
p := pfRemain;
If p >1 then
Begin
str(p,s);
chk *= p;
result += '*'+s;
end;
str(chk,s);
result += '_chk_'+s+'<';
str(pfSumOfDivs,s);
result += '_SoD_'+s+'<';
end;
end;
function smplPrimeDecomp(n:Uint64):tprimeFac;
var
pr,i,pot,fac,q :NativeUInt;
Begin
with result do
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfRemain := n;
pfRadical := 1;
pfPrimeCnt := 0;
pfpotPrimIdx[0] := 1;
pfpotMax[0] := 0;
i := 0;
while i < High(SmallPrimes) do
begin
pr := SmallPrimes[i];
q := n DIV pr;
//if n < pr*pr
if pr > q then
BREAK;
if n = pr*q then
Begin
pfpotPrimIdx[pfPrimeCnt] := i;
pot := 0;
fac := pr;
pfRadical *= pr;
repeat
n := q;
q := n div pr;
pot+=1;
fac *= pr;
until n <> pr*q;
pfpotMax[pfPrimeCnt] := pot;
pfDivCnt *= pot+1;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfPrimeCnt);
end;
inc(i);
end;
pfRemain := n;
if n > 1 then
Begin
pfRadical *= pfRemain;
pfDivCnt *= 2;
pfSumOfDivs *= n+1
end;
end;
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;
pfRadical := 1;
pfRemain := n+i;
pfPrimeCnt := 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);
pfPrimeCnt := 1;
pfRadical := 2;
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[pfPrimeCnt] := i;
pfpotMax[pfPrimeCnt] := j;
pfDivCnt *= j+1;
fac := pr;
pfRadical *= pr;
repeat
pfRemain := pfRemain DIV pr;
dec(j);
fac *= pr;
until j<= 0;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfPrimeCnt);
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
pfRadical *= j;
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;
procedure InsertSort(pDiv:tpDivisor; Left, Right : NativeInt );
var
I, J: NativeInt;
Pivot : tItem;
begin
for i:= 1 + Left to Right do
begin
Pivot:= pDiv[i];
j:= i - 1;
while (j >= Left) and (pDiv[j] > Pivot) do
begin
pDiv[j+1]:=pDiv[j];
Dec(j);
end;
pDiv[j+1]:= pivot;
end;
end;
procedure GetDivisors(pD:tpPrimeFac;var Divs:tDivisors);
var
pDivs : tpDivisor;
pPot : UInt64;
i,len,j,l,p,k: Int32;
Begin
pDivs := @Divs[0];
pDivs[0] := 1;
len := 1;
l := 1;
with pD^ do
Begin
For i := 0 to pfPrimeCnt-1 do
begin
//Multiply every divisor before with the new primefactors
//and append them to the list
k := pfpotMax[i];
p := SmallPrimes[pfpotPrimIdx[i]];
pPot :=1;
repeat
pPot *= p;
For j := 0 to len-1 do
Begin
pDivs[l]:= pPot*pDivs[j];
inc(l);
end;
dec(k);
until k<=0;
len := l;
end;
p := pfRemain;
If p >1 then
begin
For j := 0 to len-1 do
Begin
pDivs[l]:= p*pDivs[j];
inc(l);
end;
len := l;
end;
end;
//Sort. Insertsort much faster than QuickSort in this special case
InsertSort(pDivs,0,len-1);
//end marker
pDivs[len] :=0;
end;
procedure AllFacsOut(var Divs:tdivisors;proper:boolean=true);
var
k,j: Int32;
Begin
k := 0;
j := 1;
if Proper then
j:= 2;
repeat
IF Divs[j] = 0 then
BREAK;
write(Divs[k],',');
inc(j);
inc(k);
until false;
writeln(Divs[k]);
end;
const
//LIMIT =2*3*5*7*11*13*17*19*23;
LIMIT =1000*1000;
var
pPrimeDecomp :tpPrimeFac;
Mypd :tPrimeFac;
//Divs:tDivisors;
CntOfPrFac : array[0..9] of Uint32;
T0:Int64;
n : NativeUInt;
Begin
InitSmallPrimes;
T0 := GetTickCount64;
n := 0;
Init_Sieve(0);
pPrimeDecomp := @Mypd;
repeat
// Mypd := smplPrimeDecomp(n);
pPrimeDecomp:= GetNextPrimeDecomp;
if pPrimeDecomp^.pfRemain <> 1 then
inc(CntOfPrFac[pPrimeDecomp^.pfPrimeCnt+1])
else
inc(CntOfPrFac[pPrimeDecomp^.pfPrimeCnt]);
inc(n);
until n > Limit;
T0 := GetTickCount64-T0;
writeln(' Limit = ',OutPots(pPrimeDecomp,LIMIT));
writeln(' runtime ',T0/1000:0:3,' s');
For n := Low(CntOfPrFac) to High(CntOfPrFac) do
writeln(n:2,' : ',CntOfPrFac[n]:8);
end.
- @TIO.RUN fpc -O3 -XX:
Limit = 1000000 : 49 : 2^6*5^6_chk_1000000<_SoD_2480437< runtime 0.058 s //@home runtime 0.017 s 0 : 1 1 : 78735 2 : 288726 3 : 379720 4 : 208034 5 : 42492 6 : 2285 7 : 8 8 : 0 9 : 0 Real time: 0.183 s User time: 0.155 s Sys. time: 0.026 s CPU share: 99.09 % Limit = 223092870 :512 : 2*3*5*7*11*13*17*19*23_chk_223092870<_SoD_836075520< runtime 15.066 s//@home runtime 4.141 s //vs runtime 101.005 s =smplPrimeDecomp 0 : 1 1 : 12285486 2 : 48959467 3 : 76410058 4 : 58585602 5 : 22577473 6 : 4008166 7 : 262905 8 : 3712 9 : 1 Real time: 15.218 s User time: 15.068 s Sys. time: 0.033 s CPU share: 99.23 %
Phix[edit]
with javascript_semantics sequence radicals = reinstate(repeat(0,50),1,1), counts = reinstate(repeat(0,8),1,1) for i=2 to 1e6 do sequence f = vslice(prime_powers(i),1) counts[length(f)] += 1 if i<=50 then radicals[i] = product(f) end if if i=50 then printf(1,"The radicals for the first 50 positive integers are:\n%s\n", {join_by(radicals,1,10," ",fmt:="%3d")}) elsif i=99999 or i=499999 or i=999999 then printf(1,"Radical for %,7d: %,7d\n", {i, product(f)}) elsif i=1e6 then printf(1,"\nBreakdown of numbers of distinct prime factors\n") printf(1,"for positive integers from 1 to 1,000,000:\n") for c=1 to 7 do printf(1," %d: %,8d\n", {c, counts[c]}) end for printf(1," ---------\n") printf(1," %,8d\n\n", sum(counts)) end if end for integer pcount = length(get_primes_le(1e6)), ppcount = 0 for p in get_primes_le(1000) do atom p2 = p while true do p2 *= p if p2>1e6 then exit end if ppcount += 1 end while end for printf(1,"For primes or powers (>1) thereof <= 1,000,000:\n") printf(1," Number of primes = %,6d\n", pcount) printf(1," Number of powers = %,6d\n", ppcount) printf(1," Add 1 for number 1 = %,6d\n", 1) printf(1," ------\n") printf(1," %,6d\n", pcount + ppcount + 1)
Output matches other entries, but w/o a charbarchart - should you pine see here, here, here, or here.
Raku[edit]
use Prime::Factor;
use List::Divvy;
use Lingua::EN::Numbers;
sub radical ($_) { [×] unique .&prime-factors }
say "First fifty radicals:\n" ~
(1..50).map({.&radical}).batch(10)».fmt("%2d").join: "\n";
say '';
printf "Radical for %7s => %7s\n", .&comma, comma .&radical
for 99999, 499999, 999999;
my %rad = 1 => 1;
my $limit = 1e6.Int;
%rad.push: $_ for (2..$limit).race(:1000batch).map: {(unique .&prime-factors).elems => $_};
say "\nRadical factor count breakdown, 1 through {comma $limit}:";
say .key ~ " => {comma +.value}" for sort %rad;
my @primes = (2..$limit).grep: &is-prime;
my int $powers;
@primes.&upto($limit.sqrt.floor).map: -> $p {
for (2..*) { ($p ** $_) < $limit ?? ++$powers !! last }
}
say qq:to/RADICAL/;
Up to {comma $limit}:
Primes: {comma +@primes}
Powers: $powers
Plus 1: 1
Total: {comma 1 + $powers + @primes}
RADICAL
- Output:
First fifty radicals: 1 2 3 2 5 6 7 2 3 10 11 6 13 14 15 2 17 6 19 10 21 22 23 6 5 26 3 14 29 30 31 2 33 34 35 6 37 38 39 10 41 42 43 22 15 46 47 6 7 10 Radical for 99,999 => 33,333 Radical for 499,999 => 3,937 Radical for 999,999 => 111,111 Radical factor count breakdown, 1 through 1,000,000: 1 => 78,735 2 => 288,726 3 => 379,720 4 => 208,034 5 => 42,492 6 => 2,285 7 => 8 Up to 1,000,000: Primes: 78,498 Powers: 236 Plus 1: 1 Total: 78,735
RPL[edit]
PDIV is defined at Prime decomposition
NODUP is defined at Remove duplicate elements
Standard version
≪ PDIV NODUP 1 1 3 PICK SIZE FOR j OVER j GET * NEXT SWAP DROP ≫ 'RADIX' STO
HP-48G and later models
≪ PDIV NODUP ΠLIST ≫ 'RADIX' STO
Calls to solve the task:
≪ { 1 } 2 50 FOR n n RADIX + NEXT ≫ EVAL
99999 RADIX
499999 RADIX
999999 RADIX
≪ { 7 } 0 CON 1 1 PUT 2 10000 FOR n n PDIV NODUP SIZE DUP2 GET 1 + PUT NEXT ≫ EVAL
Looping one million times would need days to run on HP calculators and would wake up emulator's watchdog timer.
- Output:
5: { 1 2 3 2 5 6 7 2 3 10 11 6 13 14 15 2 17 6 19 10 21 22 23 6 5 26 3 14 29 30 31 2 33 34 35 6 37 38 39 10 41 42 43 22 15 46 47 6 7 10 }
4: 33333
3: 3937
2: 111111
1: { 1281 4097 3695 894 33 0 0 }
Wren[edit]
import "./math" for Int, Nums
import "./seq" for Lst
import "./fmt" for Fmt
var radicals = List.filled(51, 0)
radicals[1] = 1
var counts = List.filled(8, 0)
counts[1] = 1
for (i in 2..1e6) {
var factors = Lst.prune(Int.primeFactors(i))
var fc = factors.count
counts[fc] = counts[fc] + 1
if (i <= 50) radicals[i] = Nums.prod(factors)
if (i == 50) {
System.print("The radicals for the first 50 positive integers are:")
Fmt.tprint("$2d ", radicals.skip(1), 10)
System.print()
} else if (i == 99999 || i == 499999 || i == 999999) {
Fmt.print("Radical for $,7d: $,7d", i, Nums.prod(factors))
} else if (i == 1e6) {
System.print("\nBreakdown of numbers of distinct prime factors")
System.print("for positive integers from 1 to 1,000,000:")
for (i in 1..7) {
Fmt.print(" $d: $,8d", i, counts[i])
}
Fmt.print(" ---------")
Fmt.print(" $,8d", Nums.sum(counts))
Fmt.print("\nor graphically:")
Nums.barChart("", 50, Nums.toStrings(1..7), counts[1..-1])
}
}
var pcount = Int.primeCount(1e6)
var ppcount = 0
var primes1k = Int.primeSieve(1000)
for (p in primes1k) {
var p2 = p
while (true) {
p2 = p2 * p
if (p2 > 1e6) break
ppcount = ppcount + 1
}
}
Fmt.print("\nFor primes or powers (>1) thereof <= 1,000,000:")
Fmt.print(" Number of primes = $,6d", pcount)
Fmt.print(" Number of powers = $,6d", ppcount)
Fmt.print(" Add 1 for number 1 = $,6d", 1)
Fmt.print(" ------")
Fmt.print(" $,6d", pcount + ppcount + 1)
- Output:
The radicals for the first 50 positive integers are:
1 2 3 2 5 6 7 2 3 10
11 6 13 14 15 2 17 6 19 10
21 22 23 6 5 26 3 14 29 30
31 2 33 34 35 6 37 38 39 10
41 42 43 22 15 46 47 6 7 10
Radical for 99,999: 33,333
Radical for 499,999: 3,937
Radical for 999,999: 111,111
Breakdown of numbers of distinct prime factors
for positive integers from 1 to 1,000,000:
1: 78,735
2: 288,726
3: 379,720
4: 208,034
5: 42,492
6: 2,285
7: 8
---------
1,000,000
or graphically:
--------------------------------------------------
1 ■■■■■■■■ 78735
2 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 288726
3 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 379720
4 ■■■■■■■■■■■■■■■■■■■■■■ 208034
5 ■■■■ 42492
6 ◧ 2285
7 ◧ 8
--------------------------------------------------
For primes or powers (>1) thereof <= 1,000,000:
Number of primes = 78,498
Number of powers = 236
Add 1 for number 1 = 1
------
78,735
XPL0[edit]
include xpllib; \for Print and IsPrime
proc Radical(N); \Return radical of N
int N, D, D0, P;
[D:= 2; D0:= 0; P:= 1;
while N >= D*D do
[while rem(N/D) = 0 do
[if D # D0 then
[P:= P*D;
D0:= D;
];
N:= N/D;
];
D:= D+1;
];
if D # D0 then P:= P*N;
return P;
];
func DistinctFactors(N); \Return count of distinct factors of N
int N, D, D0, C;
[D:= 2; D0:= 0; C:= 0;
while N >= D*D do
[while rem(N/D) = 0 do
[if D # D0 then
[C:= C+1;
D0:= D;
];
N:= N/D;
];
D:= D+1;
];
if D # D0 and N # 1 then C:= C+1;
return C;
];
int N, C, A(10), PC, PPC, P2, P;
[Print("The radicals for the first 50 positive integers are:\n");
for N:= 1 to 50 do
[Print("%4d", Radical(N));
if rem(N/10) = 0 then CrLf(0);
];
Print("\n");
Print("Radical for %6,d: %6,d\n", 99_999, Radical( 99_999));
Print("Radical for %6,d: %6,d\n", 499_999, Radical(499_999));
Print("Radical for %6,d: %6,d\n", 999_999, Radical(999_999));
for N:= 0 to 9 do A(N):= 0;
for N:= 1 to 1_000_000 do
[C:= DistinctFactors(N);
A(C):= A(C)+1;
];
Print("\nBreakdown of numbers of distinct prime factors
for positive integers from 1 to 1,000,000:\n");
C:= 0;
for N:= 0 to 9 do
[if A(N) > 0 then
Print(" %d: %6,d\n", N, A(N));
C:= C + A(N);
];
Print(" ---------\n %,d\n", C);
\Bonus (algorithm from Wren):
PC:= 0;
for N:= 1 to 1_000_000 do
if IsPrime(N) then PC:= PC+1;
Print("\nNumber of primes: %5,d\n", PC);
PPC:= 0;
for P:= 1 to sqrt(1_000_000) do
[if IsPrime(P) then
[P2:= P;
loop [P2:= P2 * P;
if P2 > 1_000_000 then quit;
PPC:= PPC+1;
];
];
];
Print("Number of powers: %5,d\n", PPC);
if IsPrime(N) then PC:= PC+1;
Print("Total : %5,d\n", PC+PPC);
]- Output:
The radicals for the first 50 positive integers are:
1 2 3 2 5 6 7 2 3 10
11 6 13 14 15 2 17 6 19 10
21 22 23 6 5 26 3 14 29 30
31 2 33 34 35 6 37 38 39 10
41 42 43 22 15 46 47 6 7 10
Radical for 99,999: 33,333
Radical for 499,999: 3,937
Radical for 999,999: 111,111
Breakdown of numbers of distinct prime factors
for positive integers from 1 to 1,000,000:
0: 1
1: 78,734
2: 288,726
3: 379,720
4: 208,034
5: 42,492
6: 2,285
7: 8
---------
1,000,000
Number of primes: 78,498
Number of powers: 236
Total : 78,734