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

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

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

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

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 %

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.

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

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

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