Radical of an integer

You are encouraged to solve this task according to the task description, using any language you may know.
- 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
ALGOL 68
When running this with ALGOL 68G, a large heap size will need to be specified, e.g.: with -heap 256M
on the command line.
Builds tables of radicals and unique prime factor counts to avoid factorising the numbers.
BEGIN # find the radicals of some integers - the radical of n is the product #
# of thr the distinct prime factors of n, the radical of 1 is 1 #
INT max number = 1 000 000; # maximum number we will consider #
[ 1 : max number ]INT upfc; # unique prime factor counts #
[ 1 : max number ]INT radical; # table of radicals #
FOR i TO UPB upfc DO upfc[ i ] := 0; radical[ i ] := 1 OD;
FOR i FROM 2 TO UPB upfc DO
IF upfc[ i ] = 0 THEN
radical[ i ] := i;
upfc[ i ] := 1;
FOR j FROM i + i BY i TO UPB upfc DO
upfc[ j ] +:= 1;
radical[ j ] *:= i
OD
FI
OD;
# show the radicals of the first 50 positive integers #
print( ( "Radicals of 1 to 50:", newline ) );
FOR i TO 50 DO
print( ( whole( radical[ i ], -5 ) ) );
IF i MOD 10 = 0 THEN print( ( newline ) ) FI
OD;
# radicals of some specific numbers #
print( ( newline ) );
print( ( "Radical of 99999: ", whole( radical[ 99999 ], 0 ), newline ) );
print( ( "Radical of 499999: ", whole( radical[ 499999 ], 0 ), newline ) );
print( ( "Radical of 999999: ", whole( radical[ 999999 ], 0 ), newline ) );
print( ( newline ) );
# find the maximum number of unique prime factors #
INT max upfc := 0;
# show the distribution of the unique prime factor counts #
FOR i TO UPB upfc DO IF upfc[ i ] > max upfc THEN max upfc := upfc[ i ] FI OD;
[ 0 : max upfc ]INT dpfc; FOR i FROM LWB dpfc TO UPB dpfc DO dpfc[ i ] := 0 OD;
FOR i TO UPB upfc DO
dpfc[ upfc[ i ] ] +:= 1
OD;
print( ( "Distribution of radicals:", newline ) );
FOR i FROM LWB dpfc TO UPB dpfc DO
print( ( whole( i, -2 ), ": ", whole( dpfc[ i ], 0 ), newline ) )
OD
END
- Output:
Radicals of 1 to 50: 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 of 99999: 33333 Radical of 499999: 3937 Radical of 999999: 111111 Distribution of radicals: 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
Arturo
radical: $[n]-> product unique factors.prime n
(1..50) | map => radical
| split.every:10
| map => [join map & 'r -> pad to :string r 4]
| print.lines
print ""
loop [99999, 499999, 999999] 'r ->
print ["Radical of" r "->" radical r]
print ""
(1..1000000) | map => [size unique factors.prime &]
| tally
| loop [k,v] -> print [k ":" v]
- Output:
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 of 99999 -> 33333 Radical of 499999 -> 3937 Radical of 999999 -> 111111 0 : 1 1 : 78734 2 : 288726 3 : 379720 4 : 208034 5 : 42492 6 : 2285 7 : 8
C
#include <stdio.h>
int radical(int num, int *dps_out) {
int rad = 1, dps = 0, div;
if ((num & 1) == 0) {
while ((num & 1) == 0) num >>= 1;
rad *= 2;
dps++;
}
for (div = 3; div <= num; div += 2) {
if (num % div != 0) continue;
rad *= div;
dps++;
while (num % div == 0) num /= div;
}
if (dps_out != NULL) *dps_out = dps;
return rad;
}
void show_first_50() {
int n;
printf("Radicals of 1..50:\n");
for (n = 1; n <= 50; n++) {
printf("%5d", radical(n, NULL));
if (n % 5 == 0) printf("\n");
}
}
void show_rad(int n) {
printf("The radical of %d is %d.\n", n, radical(n, NULL));
}
void find_distribution() {
int n, dps, dist[8] = {0};
for (n = 1; n <= 1000000; n++) {
printf("%d\r", n);
radical(n, &dps);
dist[dps]++;
}
printf("Distribution of radicals:\n");
for (dps = 0; dps < 8; dps++) {
printf("%d: %d\n", dps, dist[dps]);
}
}
int main() {
show_first_50();
printf("\n");
show_rad(99999);
show_rad(499999);
show_rad(999999);
printf("\n");
find_distribution();
return 0;
}
- Output:
Radicals of 1..50: 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 The radical of 99999 is 33333. The radical of 499999 is 3937. The radical of 999999 is 111111. Distribution of radicals: 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
C++
#include <iomanip>
#include <iostream>
#include <map>
#include <vector>
int radical(int n, int& count) {
if (n == 1) {
count = 1;
return 1;
}
int product = 1;
count = 0;
if ((n & 1) == 0) {
product *= 2;
++count;
while ((n & 1) == 0)
n >>= 1;
}
for (int p = 3, sq = 9; sq <= n; p += 2) {
if (n % p == 0) {
product *= p;
++count;
while (n % p == 0)
n /= p;
}
sq += (p + 1) << 2;
}
if (n > 1) {
product *= n;
++count;
}
return product;
}
std::vector<bool> prime_sieve(int limit) {
std::vector<bool> sieve(limit, true);
if (limit > 0)
sieve[0] = false;
if (limit > 1)
sieve[1] = false;
for (int i = 4; i < limit; i += 2)
sieve[i] = false;
for (int p = 3, sq = 9; sq < limit; p += 2) {
if (sieve[p]) {
for (int q = sq; q < limit; q += p << 1)
sieve[q] = false;
}
sq += (p + 1) << 2;
}
return sieve;
}
int main() {
std::cout.imbue(std::locale(""));
std::cout << "Radicals of the first 50 positive integers:\n";
int count = 0;
for (int i = 1; i <= 50; ++i) {
std::cout << std::setw(2) << radical(i, count)
<< (i % 10 == 0 ? '\n' : ' ');
}
std::cout << '\n';
for (int i : {99999, 499999, 999999}) {
std::cout << "Radical of " << std::setw(7) << i << " is "
<< std::setw(7) << radical(i, count) << ".\n";
}
std::map<int, int> prime_factors;
const int limit = 1000000;
for (int i = 1; i <= limit; ++i) {
radical(i, count);
++prime_factors[count];
}
std::cout << "\nRadical factor count breakdown up to " << limit << ":\n";
for (const auto& [count, total] : prime_factors) {
std::cout << count << ": " << std::setw(7) << total << '\n';
}
auto sieve = prime_sieve(limit + 1);
int primes = 0, prime_powers = 0;
for (int i = 1; i <= limit; ++i) {
if (sieve[i]) {
++primes;
int n = limit;
while (i <= n / i) {
++prime_powers;
n /= i;
}
}
}
std::cout << "\nUp to " << limit << ":\n";
std::cout << "Primes: " << std::setw(6) << primes
<< "\nPowers: " << std::setw(6) << prime_powers
<< "\nPlus 1: " << std::setw(6) << 1
<< "\nTotal: " << std::setw(6) << primes + prime_powers + 1
<< '\n';
}
- Output:
Radicals of the 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 of 99,999 is 33,333. Radical of 499,999 is 3,937. Radical of 999,999 is 111,111. Radical factor count breakdown up to 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
CLU
distinct_factors = iter (n: int) yields (int)
if n // 2 = 0 then
yield(2)
while n // 2 = 0 do n := n / 2 end
end
div: int := 3
while div <= n do
if n // div = 0 then
yield(div)
while n // div = 0 do n := n / div end
end
div := div + 2
end
end distinct_factors
radical = proc (n: int) returns (int)
rad: int := 1
for fac: int in distinct_factors(n) do
rad := rad * fac
end
return(rad)
end radical
dpf_count = proc (n: int) returns (int)
count: int := 0
for fac: int in distinct_factors(n) do
count := count + 1
end
return(count)
end dpf_count
show_first_50 = proc (out: stream)
stream$putl(out, "Radicals of 1..50:")
for n: int in int$from_to(1, 50) do
stream$putright(out, int$unparse(radical(n)), 5)
if n // 10 = 0 then stream$putl(out, "") end
end
end show_first_50
show_radical = proc (out: stream, n: int)
stream$puts(out, "The radical of ")
stream$puts(out, int$unparse(n))
stream$puts(out, " is ")
stream$puts(out, int$unparse(radical(n)))
stream$putl(out, ".")
end show_radical
find_distribution = proc (out: stream)
dist: array[int] := array[int]$fill(0, 8, 0)
for n: int in int$from_to(1, 1000000) do
stream$putright(out, int$unparse(n) || "\r", 8)
dpf: int := dpf_count(n)
dist[dpf] := dist[dpf] + 1
end
stream$putl(out, "Distribution of radicals:")
for dpf: int in array[int]$indexes(dist) do
stream$putright(out, int$unparse(dpf), 1)
stream$puts(out, ": ")
stream$putright(out, int$unparse(dist[dpf]), 8)
stream$putl(out, "")
end
end find_distribution
start_up = proc ()
po: stream := stream$primary_output()
show_first_50(po)
stream$putl(po, "")
for n: int in array[int]$elements(array[int]$[99999, 499999, 999999]) do
show_radical(po, n)
end
stream$putl(po, "")
find_distribution(po)
end start_up
- Output:
Radicals of 1..50: 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 The radical of 99999 is 33333. The radical of 499999 is 3937. The radical of 999999 is 111111. Distribution of radicals: 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
Cowgol
include "cowgol.coh";
sub radical(n: uint32): (rad: uint32, dpf: uint8) is
dpf := 0;
rad := 1;
# Handle factors of 2 first
if n & 1 == 0 then
dpf := 1;
rad := 2;
while n > 0 and n & 1 == 0 loop
n := n >> 1;
end loop;
end if;
# Divide out odd prime factors
var div: uint32 := 3;
while div <= n loop
if n % div == 0 then
dpf := dpf + 1;
rad := rad * div;
while n % div == 0 loop
n := n / div;
end loop;
end if;
div := div + 2;
end loop;
end sub;
sub first50() is
print("Radicals of 1..50:\n");
var n: uint32 := 1;
var dpf: uint8;
var rad: uint32;
while n <= 50 loop
(rad, dpf) := radical(n);
print_i32(rad);
if n % 5 != 0 then
print_char('\t');
else
print_nl();
end if;
n := n + 1;
end loop;
end sub;
sub print_rad(n: uint32) is
var dpf: uint8;
var rad: uint32;
print("The radical of ");
print_i32(n);
print(" is ");
(rad, dpf) := radical(n);
print_i32(rad);
print(".\n");
end sub;
sub find_distribution() is
var dist: uint32[8];
MemZero(&dist[0] as [uint8], @bytesof dist);
var n: uint32 := 1;
var dpf: uint8;
var rad: uint32;
while n <= 1000000 loop
print_i32(n);
print_char('\r');
(rad, dpf) := radical(n);
dist[dpf] := dist[dpf] + 1;
n := n + 1;
end loop;
print("Distribution of radicals:\n");
dpf := 0;
while dpf < 8 loop
print_i8(dpf);
print(": ");
print_i32(dist[dpf]);
print_nl();
dpf := dpf + 1;
end loop;
end sub;
first50();
print_nl();
print_rad(99999);
print_rad(499999);
print_rad(999999);
print_nl();
find_distribution();
- Output:
Radicals of 1..50: 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 The radical of 99999 is 33333. The radical of 499999 is 3937. The radical of 999999 is 111111. Distribution of radicals: 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
EasyLang
fastfunc radnf num .
d = 2
while d <= sqrt num
if num mod d = 0
nf += 1
while num mod d = 0
num /= d
.
.
d += 1
.
if d <= num
nf += 1
.
return nf
.
func rad num .
r = 1
d = 2
while d <= sqrt num
if num mod d = 0
r *= d
while num mod d = 0
num /= d
.
.
d += 1
.
if d <= num
r *= num
.
return r
.
proc show50 . .
write "First 50 radicals:"
for n = 1 to 50
write " " & rad n
.
print ""
.
proc show n . .
print "radical(" & n & ") = " & rad n
.
proc dist . .
len dist[] 7
for n = 2 to 1000000
dist[radnf n] += 1
.
print ""
print "Distribution of radicals:"
print "0: 1"
for i = 1 to 7
print i & ": " & dist[i]
.
.
show50
print ""
show 99999
show 499999
show 999999
print ""
dist
FreeBASIC
'#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
~.&.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
Java
import java.util.ArrayList;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.IntStream;
public final class RadicalOfAnInteger {
public static void main(String[] aArgs) {
sievePrimes();
distinctPrimeFactors();
System.out.println("Radical of first 50 positive integers:");
for ( int n = 1; n <= 50; n++ ) {
System.out.print(String.format("%3d%s", radical(n), ( n % 10 == 0 ? "\n" : "" )));
}
System.out.println();
for ( int n : List.of( 99_999, 499_999, 999_999 ) ) {
System.out.println(String.format("%s%7d%s%7d", "Radical for ", n, ":", radical(n)));
}
System.out.println();
System.out.print("Distribution of the first one million positive ");
System.out.println("integers by numbers of distinct prime factors:");
List<Integer> counts = IntStream.rangeClosed(0, 7).boxed().map( i -> 0 ).collect(Collectors.toList());
for ( int n = 1; n <= LIMIT; n++ ) {
counts.set(distinctPrimeFactorCount(n), counts.get(distinctPrimeFactorCount(n)) + 1);
}
for ( int n = 0; n < counts.size(); n++ ) {
System.out.println(String.format("%d%s%7d", n, ":", counts.get(n)));
}
System.out.println();
System.out.println("Number of primes and powers of primes less than or equal to one million:");
int count = 0;
final double logOfLimit = Math.log(LIMIT);
for ( int p : primes ) {
count += logOfLimit / Math.log(p);
}
System.out.println(count);
}
private static int radical(int aNumber) {
return distinctPrimeFactors.get(aNumber).stream().reduce(1, Math::multiplyExact);
}
private static int distinctPrimeFactorCount(int aNumber) {
return distinctPrimeFactors.get(aNumber).size();
}
private static void distinctPrimeFactors() {
distinctPrimeFactors = IntStream.rangeClosed(0, LIMIT).boxed()
.map( i -> new ArrayList<Integer>() ).collect(Collectors.toList());
for ( int p : primes ) {
for ( int n = p; n <= LIMIT; n += p ) {
distinctPrimeFactors.get(n).add(p);
}
}
}
private static void sievePrimes() {
List<Boolean> markedPrime = IntStream.rangeClosed(0, LIMIT).boxed()
.map( i -> true ).collect(Collectors.toList());
for ( int p = 2; p * p <= LIMIT; p++ ) {
if ( markedPrime.get(p) ) {
for ( int i = p * p; i <= LIMIT; i += p ) {
markedPrime.set(i, false);
}
}
}
primes = new ArrayList<Integer>();
for ( int p = 2; p <= LIMIT; p++ ) {
if ( markedPrime.get(p) ) {
primes.add(p);
}
}
}
private static List<Integer> primes;
private static List<List<Integer>> distinctPrimeFactors;
private static final int LIMIT = 1_000_000;
}
- 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 99999: 33333 Radical for 499999: 3937 Radical for 999999: 111111 Distribution of the first one million positive integers by numbers of distinct prime factors: 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8 Number of primes and powers of primes less than or equal to one million: 78734
jq
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, task2
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 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
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
Kotlin
import java.util.BitSet
import java.util.TreeMap
import kotlin.math.sqrt
class Radicals(val limit: Int) {
val primes: List<Int>
val radicals: IntArray // valid indices: 1..limit
val distinctPrimeFactorCounts: IntArray // valid indices: 2..limit
val distinctPrimeFactorCountDistribution: Map<Int, Int>
init {
// Sieve
val composites = BitSet()
for (i in 2..sqrt(limit.toDouble()).toInt()) {
if (!composites[i]) {
for (j in i * i..limit step i) composites.set(j)
}
}
primes = (2..limit).filterNot(composites::get)
// Calculate radicals and distinct prime factor counts
distinctPrimeFactorCounts = IntArray(limit + 1) { 0 }
radicals = IntArray(limit + 1) { 1 }
for (p in primes) {
for (i in p..limit step p) {
distinctPrimeFactorCounts[i]++
radicals[i] *= p
}
}
distinctPrimeFactorCountDistribution =
distinctPrimeFactorCounts.asList()
.subList(2, distinctPrimeFactorCounts.size)
.groupingBy { it }
.eachCountTo(TreeMap())
}
}
fun main() {
with(Radicals(limit = 1_000_000)) {
println("Radicals of first 50 positive integers:")
for (i in 1..41 step 10) {
println((i..i + 9).joinToString(separator = " ") { "%2d".format(radicals[it]) })
}
println()
for (n in listOf(99_999, 499_999, 999_999)) {
println("Radical for %6d: %6d".format(n, radicals[n]))
}
println()
println("Distribution of the first $limit positive integers by numbers of distinct prime factors:")
for ((i, count) in distinctPrimeFactorCountDistribution) {
println("%d: %6d".format(i, count))
}
println()
println("Number of primes and powers of primes less than or equal to $limit:")
val count = primes.sumOf { prime ->
generateSequence(limit) { it / prime }.takeWhile { it >= prime }.count()
}
println(count)
}
}
- Output:
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 Radical for 99999: 33333 Radical for 499999: 3937 Radical for 999999: 111111 Distribution of the first 1000000 positive integers by numbers of distinct prime factors: 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8 Number of primes and powers of primes less than or equal to 1000000: 78734
Lua
do -- find the radicals of some integers - the radical of n is the product
-- of the distinct prime factors of n, the radical of 1 is 1
local maxNumber = 1000000; -- maximum number we will consider
local upfc, radical = {}, {} -- unique prime factor counts and radicals
for i = 1, maxNumber do
upfc[ i ] = 0
radical[ i ] = 1
end
for i = 2, maxNumber do
if upfc[ i ] == 0 then
radical[ i ] = i
upfc[ i ] = 1
for j = i + i, maxNumber, i do
upfc[ j ] = upfc[ j ] + 1
radical[ j ] = radical[ j ] * i
end
end
end
-- show the radicals of the first 50 positive integers
io.write( "Radicals of 1 to 50:\n" )
for i = 1, 50 do
io.write( string.format( "%5d", radical[ i ] ), ( i % 10 == 0 and "\n" or "" ) )
end
-- radicals of some specific numbers
io.write( "\n" )
io.write( "Radical of 99999: ", radical[ 99999 ], "\n" )
io.write( "Radical of 499999: ", radical[ 499999 ], "\n" )
io.write( "Radical of 999999: ", radical[ 999999 ], "\n" )
io.write( "\n" )
-- show the distribution of the unique prime factor counts
local dpfc = {}
for i = 1, maxNumber do
local count = dpfc[ upfc[ i ] ]
dpfc[ upfc[ i ] ] = ( count == nil and 1 or count + 1 )
end
io.write( "Distribution of radicals:\n" )
for i = 0, #dpfc do
io.write( string.format( "%2d", i ), ": ", dpfc[ i ], "\n" )
end
end
- Output:
Radicals of 1 to 50: 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 of 99999: 33333 Radical of 499999: 3937 Radical of 999999: 111111 Distribution of radicals: 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
MAD
NORMAL MODE IS INTEGER
INTERNAL FUNCTION(A,B)
ENTRY TO REM.
FUNCTION RETURN A-A/B*B
END OF FUNCTION
INTERNAL FUNCTION(NN)
ENTRY TO RADIC.
KN = NN
DPF = 0
RAD = 1
WHENEVER KN.E.0, FUNCTION RETURN 0
WHENEVER REM.(KN,2).E.0
DPF = 1
RAD = 2
THROUGH DIV2, FOR KN=KN, 0, REM.(KN,2).NE.0
DIV2 KN = KN / 2
END OF CONDITIONAL
THROUGH ODD, FOR FAC=3, 2, FAC.G.KN
WHENEVER REM.(KN,FAC).E.0
DPF = DPF + 1
RAD = RAD * FAC
THROUGH DIVP, FOR KN=KN, 0, REM.(KN,FAC).NE.0
DIVP KN = KN / FAC
END OF CONDITIONAL
ODD CONTINUE
FUNCTION RETURN RAD
END OF FUNCTION
R PRINT RADICALS FOR 1..50
VECTOR VALUES RADLIN = $10(I5)*$
PRINT COMMENT $ RADICALS OF 1 TO 50$
THROUGH RADL, FOR L=0, 10, L.GE.50
RADL PRINT FORMAT RADLIN,
0 RADIC.(L+1),RADIC.(L+2),RADIC.(L+3),RADIC.(L+4),
1 RADIC.(L+5),RADIC.(L+6),RADIC.(L+7),RADIC.(L+8),
2 RADIC.(L+9),RADIC.(L+10)
R PRINT RADICALS OF TASK NUMBERS
PRINT COMMENT $ $
VECTOR VALUES RADNUM = $14HTHE RADICAL OF,I7,S1,2HIS,I7*$
THROUGH RADN, FOR VALUES OF N = 99999, 499999, 999999
RADN PRINT FORMAT RADNUM,N,RADIC.(N)
R FIND DISTRIBUTION
PRINT COMMENT $ $
DIMENSION DIST(8)
THROUGH ZERO, FOR D = 0, 1, D.G.7
ZERO DIST(D) = 0
THROUGH CNTDST, FOR N = 1, 1, N.G.1000000
RADIC.(N)
CNTDST DIST(DPF) = DIST(DPF) + 1
PRINT COMMENT $ DISTRIBUTION OF RADICALS$
VECTOR VALUES DISTLN = $I1,2H: ,I8*$
THROUGH SHWDST, FOR D = 0, 1, D.G.7
SHWDST PRINT FORMAT DISTLN,D,DIST(D)
END OF PROGRAM
- Output:
RADICALS OF 1 TO 50 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 THE RADICAL OF 99999 IS 33333 THE RADICAL OF 499999 IS 3937 THE RADICAL OF 999999 IS 111111 DISTRIBUTION OF RADICALS 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
Nim
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
Mathematica /Wolfram Language
ClearAll[rad];
primeFactors[n_Integer] := First@Transpose@FactorInteger[n];
rad[n_Integer] := Times @@ primeFactors[n];
Print[TableForm[Partition[rad /@ Range[50], 10] ,
TableAlignments -> Right]];
Print[];
Print[TableForm[
List @@@ Normal@CountsBy[Length[primeFactors[#]] &][Range[10^6]] ,
TableAlignments -> Right]];
- Output:
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 1 78735 2 288726 3 379720 4 208034 5 42492 6 2285 7 8
Pascal
Free Pascal
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
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 %
PascalABC.NET
function Factors(N: integer): List<integer>;
begin
var lst := new List<integer>;
var i := 2;
while i * i <= N do
begin
while N.Divs(i) do
begin
if i not in lst then
lst.Add(i);
N := N div i;
end;
i += 1;
end;
if N >= 2 then
lst.Add(N);
Result := lst;
end;
function Radical(x: integer) := Factors(x).Product;
begin
for var i:=1 to 50 do
begin
Write(Radical(i):3);
if i mod 10 = 0 then
Writeln;
end;
Writeln;
Writeln('Radical(99999) = ',Radical(99999));
Writeln('Radical(499999) = ',Radical(499999));
Writeln('Radical(999999) = ',Radical(999999));
Writeln;
var a := |0| *8;
for var i:=1 to 1000000 do
a[Factors(i).Count] += 1;
for var i:=0 to 7 do
Writeln($'{i}: {a[i]}');
end.
- Output:
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(99999) = 33333 Radical(499999) = 3937 Radical(999999) = 111111 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
Perl
use strict;
use warnings;
use feature <say signatures>;
no warnings 'experimental::signatures';
use List::Util 'uniq';
use ntheory <primes vecprod factor>;
sub comma { reverse ((reverse shift) =~ s/.{3}\K/,/gr) =~ s/^,//r }
sub table (@V) { ( sprintf( ('%3d')x@V, @V) ) =~ s/.{1,30}\K/\n/gr }
sub radical ($n) { vecprod uniq factor $n }
my $limit = 1e6;
my $primes = primes(1,$limit);
my %rad;
$rad{1} = 1;
$rad{ uniq factor $_ }++ for 2..$limit;
my $powers;
my $upto = int sqrt $limit;
for my $p ( grep { $_< $upto} @$primes ) {
for (2..$upto) { ($p ** $_) < $limit ? ++$powers : last }
}
say 'First fifty radicals:';
say table map { radical $_ } 1..50;
printf "Radical for %7s => %7s\n", comma($_), comma radical($_) for 99999, 499999, 999999;
printf "\nRadical factor count breakdown, 1 through %s:\n", comma $limit;
say "$_ => " . comma $rad{$_} for sort keys %rad;
say <<~"END";
Up to @{[comma $limit]}:
Primes: @{[comma 0+@$primes]}
Powers: $powers
Plus 1: 1
Total: @{[comma 1 + $powers + @$primes]}
END
- 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
Phix
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, here, or here.
Python
# int_radicals.py by Xing216
def radical(n):
product = 1
if (n % 2 == 0):
product *= 2
while (n%2 == 0):
n = n/2
for i in range (3, int((n)**0.5), 2):
if (n % i == 0):
product = product * i
while (n%i == 0):
n = n/i
if (n > 2):
product = product * n
return int(product)
def distinctPrimeFactors(N):
factors = []
if (N < 2):
factors.append(-1)
return factors
if N == 2:
factors.append(2)
return factors
visited = {}
i = 2
while(i * i <= N):
while(N % i == 0):
if(i not in visited):
factors.append(i)
visited[i] = 1
N //= i
i+=1
if(N > 2):
factors.append(N)
return factors
if __name__ == "__main__":
print("Radical of first 50 positive integers:")
for i in range(1,51):
print(f"{radical(i):>2}", end=" ")
if (i % 10 == 0):
print()
print()
for n in [99999, 499999, 999999]:
print(f"Radical of {n:>6}: {radical(n)}")
distDict = {1:0,2:0,3:0,4:0,5:0,6:0,7:0}
for i in range(1,1000000):
distinctPrimeFactorCount = len(distinctPrimeFactors(i))
distDict[distinctPrimeFactorCount] += 1
print("\nDistribution of the first one million positive integers by numbers of distinct prime factors:")
for key, value in distDict.items():
print(f"{key}: {value:>6}")
print("\nNumber of primes and powers of primes less than or equal to one million:")
print(distDict[1])
- Output:
Radical of first 50 positive integers: 1 2 3 2 5 6 7 2 9 10 11 6 13 14 15 2 17 18 19 10 21 22 23 6 25 26 3 14 29 30 31 2 33 34 35 18 37 38 39 10 41 42 43 22 15 46 47 6 49 50 Radical of 99999: 33333 Radical of 499999: 3937 Radical of 999999: 111111 Distribution of the first one million positive integers by numbers of distinct prime factors: 1: 78735 2: 288725 3: 379720 4: 208034 5: 42492 6: 2285 7: 8 Number of primes and powers of primes less than or equal to one million: 78735
Raku
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
REXX
Version 1
/* REXX */
arg n
if n = '' then
n = 1000000
Numeric Digits 100
parse version version; say version digits() 'Digits'
say 'Radical of an integer = product of distinct prime factors'
say 'Version 1: tailor-made routines'
say
call setp
/* 1. Radicals from 1 to 50 */
Say 'Radicals of 1..50:'
ol=''
Do i=1 To 50
ol=ol||right(rad(i),5)
If i//10=0 then Do
Say ol
ol=''
End
End
Say ''
Call radt 99999
Call radt 499999
Call radt 999999
Say ''
Say 'Distribution of radicals:'
cnt.=0
m = n/10
Say ''
Say 'Getting distribution list...'
Call time 'R'
Do v=1 To n
p=rad(v,'cnt')
cnt.p+=1
If v//m=0 Then Do
ti=(v%m)*10
say format(ti,3)'%' format(time('E'),4,3) 'seconds'
End
/**************************************
If p=7 Then
Say v':' rad(v,'lstx') '=' rad(v)
**************************************/
End
Say ''
Say 'Results'
Do d=0 To 20
If cnt.d>0 Then
Say d':' right(cnt.d,8)
End
exit
radt:
Parse Arg u
Say 'The radical of' u 'is' rad(u)'.'
Return
rad:
Parse Arg z,what
zz=z
plist=''
exp.=0
Do Until p*p>z
Do pi=1 By 1 Until p*p>z
p=word(primzahlen,pi)
exp.p=0
Do while z//p=0
exp.p+=1
If pos(p,plist)=0 Then
plist=plist p
z=z%p
End
End
End
exp.z+=1
If pos(z,plist)=0 &z<>1 Then
plist=plist z
Select
When what='lst' Then
Return plist
When what='lstx' Then Do
s=''
Do While plist<>''
Parse Var plist p plist
if exp.p=1 Then
s=s'*'p
Else
s=s'*'p'**'exp.p
End
Return strip(s,,'*')
End
When what='cnt' Then
Return words(plist)
Otherwise Do
rad=1
Do while plist>''
Parse Var plist p plist
rad=rad*p
End
Return rad
End
End
setp:
primzahlen=,
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101,
103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197,
199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311,
313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431,
433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557,
563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661,
673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809,
811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937,
941 947 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049
Return
- Output:
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 100 Digits Radical of an integer = product of distinct prime factors Version 1: tailor-made routines Radicals of 1..50: 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 The radical of 99999 is 33333. The radical of 499999 is 3937. The radical of 999999 is 111111. Distribution of radicals: Getting distribution list... 10% 2.747 seconds 20% 6.474 seconds 30% 10.712 seconds 40% 15.524 seconds 50% 20.787 seconds 60% 26.516 seconds 70% 32.870 seconds 80% 39.696 seconds 90% 46.792 seconds 100% 54.265 seconds Results 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8
This program is optimized for the task (all primes up to 1000 are specified...) and works to n = 1000000 max.
Version 2
Libraries: How to use
Libraries: Source code
Below program uses generic optimized procedures for checking and finding factors and primes. It will work for any n (but for say n > 2000000 it becomes very slow).
include Settings
say version; say 'Radical of an integer'; say
arg n
if n = '' then
n = 1000000
numeric digits 100
say 'Radicals for 1..50:'
ol = ''
do i = 1 to 50
ol = ol||Right(Radical(i),5)
if i//10 = 0 then do
say ol; ol = ''
end
end
say
say 'Radicals for:'
say '99999 =' Radical(99999)
say '499999 =' Radical(499999)
say '999999 =' Radical(999999)
say
m = n/10; r = Isqrt(n); radi. = 0
say 'Getting distribution list...'
call Time('r')
do i = 1 to n
call Radical(i)
u = ufac.0
radi.u = radi.u+1
if i//m=0 then do
ti=(i%m)*10
say Format(ti,3)'%' Format(Time('e'),4,3) 'seconds'
end
end
say
say 'Distribution for first' n 'radicals over Number of Factors:'
do i = 0 to 10
if radi.i > 0 then
say Right(i,2)':' Right(radi.i,6)
end
say
say 'Getting Primes up to' n'...'
call Time('r')
pr = Primes(n)
say 'Took' Format(Time('e'),,3) 'seconds'
say
say 'Getting powers of Primes up to' r'...'
call Time('r')
pw = 0
do i = 1
p1 = prim.Prime.i
if p1 > r then
leave
p2 = p1
do forever
p2 = p2*p1
if p2 > n then
leave
pw = pw+1
end
end
say 'Took' Format(Time('e'),,3) 'seconds'
say
say 'Primes' Format(pr,6)
say 'Powers' Format(pw,6)
say ' -----'
say 'Total ' Format(pr+pw,6)
exit
Radical:
/* Radical = product of unique prime factors */
procedure expose ufac.
arg x
/* Get unique factors */
n = Ufactors(x)
/* Calculate product */
y = 1
do i = 1 to n
y = y*ufac.factor.i
end
return y
include Functions
include Numbers
include Sequences
include Abend
- Output:
REXX-ooRexx_5.0.0(MT)_64-bit 6.05 23 Dec 2022 Radical of an integer Radicals for 1..50: 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 Radicals for: 99999 = 33333 499999 = 3937 999999 = 111111 Getting distribution list... 10% 2.581 seconds 20% 5.628 seconds 30% 9.254 seconds 40% 12.769 seconds 50% 16.457 seconds 60% 20.223 seconds 70% 24.439 seconds 80% 28.736 seconds 90% 33.197 seconds 100% 37.532 seconds Distribution for first 1000000 radicals over number of factors: 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8 Getting primes up to 1000000... Took 0.909 seconds Getting powers of primes up to 1000... Took 0.000 seconds Primes 78498 Powers 236 ----- Total 78734
Version 2 is slightly faster than Version 1.
RPL
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 }
Ruby
Adding 3 methods to Integers, but not globally:
module Radical
refine Integer do
require 'prime'
def radical = self == 1 ? 1 : prime_division.map(&:first).inject(&:*)
def num_uniq_prime_factors = prime_division.size
def prime_pow? = prime_division.size == 1
end
end
using Radical
n = 50 # task 1
puts "The radicals for the first #{n} positive integers are:"
(1..n).map(&:radical).each_slice(10){|s| puts "%4d"*s.size % s}
puts # task 2
[99999, 499999 , 999999].each{|n| puts "Radical for %6d: %6d" % [n, n.radical]}
n = 1_000_000 # task 3 & bonus
puts "\nNumbers of distinct prime factors for integers from 1 to #{n}"
(1..n).map(&:num_uniq_prime_factors).tally.each{|kv| puts "%d: %8d" % kv }
puts "\nNumber of primes and powers of primes less than or equal to #{n}: #{(1..n).count(&:prime_pow?)}"
- 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 99999: 33333 Radical for 499999: 3937 Radical for 999999: 111111 Numbers of distinct prime factors for integers from 1 to 1000000 0: 1 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8 Number of primes and powers of primes less than or equal to 1000000: 78734
Sidef
Takes less than one second to run:
with (50) {|n|
say "Radicals of the first #{n} positive integers:"
1..n -> map { .rad }.each_slice(10, {|*s|
say s.map { '%2s' % _ }.join(' ')
})
}
say ''; [99999, 499999, 999999].map {|n|
say "rad(#{n}) = #{rad(n)}"
}
for limit in (1e6, 1e7, 1e8, 1e9) {
say "\nCounting of k-omega primes <= #{limit.commify}:"
for k in (1..Inf) {
break if (pn_primorial(k) > limit)
var c = k.omega_prime_count(limit)
say "#{k}: #{c}"
assert_eq(c, limit.prime_power_count) if (k == 1)
}
}
- Output:
Radicals of the 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 rad(99999) = 33333 rad(499999) = 3937 rad(999999) = 111111 Counting of k-omega primes <= 1,000,000: 1: 78734 2: 288726 3: 379720 4: 208034 5: 42492 6: 2285 7: 8 Counting of k-omega primes <= 10,000,000: 1: 665134 2: 2536838 3: 3642766 4: 2389433 5: 691209 6: 72902 7: 1716 8: 1 Counting of k-omega primes <= 100,000,000: 1: 5762859 2: 22724609 3: 34800362 4: 25789580 5: 9351293 6: 1490458 7: 80119 8: 719 Counting of k-omega primes <= 1,000,000,000: 1: 50851223 2: 206415108 3: 332590117 4: 269536378 5: 114407511 6: 24020091 7: 2124141 8: 55292 9: 138
Wren
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
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