Radical of an integer

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

• Square-free integers
• Almost prime

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

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

Translation of: XPL0

'#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

Works with: jq

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

Lua

Translation of: ALGOL 68

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

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 %

Perl

Library: ntheory

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

/* REXX */
call setp
Numeric Digits 100
/* 1. Radicals from 1 to 50 */
Say 'Radicals of 1..50:'
ol=''
Do n=1 To 50
  ol=ol||right(rad(n),5)
  If n//10=0 then Do
    Say ol
    ol=''
    End
  End

Say ''
Call radt 99999
Call radt 499999
Call radt 999999

Say ''
Say 'Distribution of radicals:'
Call time 'R'
cnt.=0
Do v=1 To 1000000
  If v//100000=1 Then Do
    ti=v%100000
    etime.ti=format(time('E'),4,1)
    End
  p=rad(v,'cnt')
  cnt.p+=1
  /**************************************
  If p=7 Then
    Say v':' rad(v,'lstx') '=' rad(v)
  **************************************/
  End
etime.10=format(time('E'),4,1)
Do d=0 To 20
  If cnt.d>0 Then
    Say d':' right(cnt.d,8)
  End
Say ''
Say 'Timings:'
Do ti=1 To 10
  Say right(ti,2) etime.ti 'seconds'
  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:

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

Timings:
 1    4.2 seconds
 2   10.4 seconds
 3   18.0 seconds
 4   26.1 seconds
 5   35.1 seconds
 6   44.6 seconds
 7   55.9 seconds
 8   68.5 seconds
 9   83.6 seconds
10  100.6 seconds

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

Library: Wren-math

Library: Wren-seq

Library: Wren-fmt

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