Numbers whose count of divisors is prime

Numbers whose count of divisors is prime is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Task

Find positive integers   n   which count of divisors is prime,   but not equal to  2,   where   n   <   1,000.

Stretch goal:   (as above),   but where   n   <   100,000.

11l

Translation of: FreeBASIC

F is_prime(a)
   I a == 2
      R 1B
   I a < 2 | a % 2 == 0
      R 0B
   L(i) (3 .. Int(sqrt(a))).step(2)
      I a % i == 0
         R 0B
   R 1B

print(‘Numbers which count of divisors is prime are:’)
V row = 0

L(n) 1..99999
   V num = 0
   L(m) 1 .. n
      I n % m == 0
         num++
   I is_prime(num) & num != 2
      print(‘#6’.format(n), end' ‘ ’)
      row++
      I row % 5 == 0
         print()

print("\n\nFound "row‘ numbers’)

Output:

Numbers which count of divisors is prime are:
     4      9     16     25     49
    64     81    121    169    289
   361    529    625    729    841
   961   1024   1369   1681   1849
  2209   2401   2809   3481   3721
  4096   4489   5041   5329   6241
  6889   7921   9409  10201  10609
 11449  11881  12769  14641  15625
 16129  17161  18769  19321  22201
 22801  24649  26569  27889  28561
 29929  32041  32761  36481  37249
 38809  39601  44521  49729  51529
 52441  54289  57121  58081  59049
 63001  65536  66049  69169  72361
 73441  76729  78961  80089  83521
 85849  94249  96721  97969

Found 79 numbers

Action!

Library: Action! Sieve of Eratosthenes

INCLUDE "H6:SIEVE.ACT"

INT FUNC CountDivisors(INT x)
  INT i,max,count

  count=2 i=2 max=x
  WHILE i<max
  DO
    IF x MOD i=0 THEN
      max=x/i
      IF i<max THEN
        count==+2
      ELSEIF i=max THEN
        count==+1
      FI
    FI
    i==+1
  OD
RETURN (count)

PROC Main()
  DEFINE MAXNUM="999"
  BYTE ARRAY primes(MAXNUM+1)

  INT i,count

  Put(125) PutE() ;clear the screen
  Sieve(primes,MAXNUM+1)
  FOR i=2 TO MAXNUM
  DO
    IF primes(i)=0 THEN
      count=CountDivisors(i)
      IF count>2 AND primes(count)=1 THEN
        PrintF("%I has %I divisors%E",i,count)
      FI
    FI
  OD
RETURN

Output:

Screenshot from Atari 8-bit computer

4 has 3 divisors
9 has 3 divisors
16 has 5 divisors
25 has 3 divisors
49 has 3 divisors
64 has 7 divisors
81 has 5 divisors
121 has 3 divisors
169 has 3 divisors
289 has 3 divisors
361 has 3 divisors
529 has 3 divisors
625 has 5 divisors
729 has 7 divisors
841 has 3 divisors
961 has 3 divisors

ALGOL 68

Counts the divisors without using division.

BEGIN # find numbers with prime divisor counts                         #
    INT max number            := 1 000;
    TO 2 DO
        INT max divisors      := 0;
        # construct a table of the divisor counts                      #
        [ 1 : max number ]INT ndc; FOR i FROM 1 TO UPB ndc DO ndc[ i ] := 1 OD;
        FOR i FROM 2 TO UPB ndc DO
            FOR j FROM i BY i TO UPB ndc DO ndc[ j ] +:= 1 OD
        OD;
        # show the numbers with prime divisor counts                   #
        print( ( "Numbers up to ", whole( max number, 0 ), " with odd prime divisor counts:", newline ) );
        INT p count := 0;
        FOR i TO UPB ndc DO
            INT divisor count = ndc[ i ];
            IF ODD divisor count AND ndc[ divisor count ] = 2 THEN
                print( ( whole( i, -8 ) ) );
                IF ( p count +:= 1 ) MOD 10 = 0 THEN print( ( newline ) ) FI
            FI
        OD;
        print( ( newline ) );
        max number := 100 000
    OD
END

Output:

Numbers up to 1000 with odd prime divisor counts:
       4       9      16      25      49      64      81     121     169     289
     361     529     625     729     841     961
Numbers up to 100000 with odd prime divisor counts:
       4       9      16      25      49      64      81     121     169     289
     361     529     625     729     841     961    1024    1369    1681    1849
    2209    2401    2809    3481    3721    4096    4489    5041    5329    6241
    6889    7921    9409   10201   10609   11449   11881   12769   14641   15625
   16129   17161   18769   19321   22201   22801   24649   26569   27889   28561
   29929   32041   32761   36481   37249   38809   39601   44521   49729   51529
   52441   54289   57121   58081   59049   63001   65536   66049   69169   72361
   73441   76729   78961   80089   83521   85849   94249   96721   97969

AppleScript

Only squares checked, as per the Discussion page.

on countFactors(n) -- Positive ns only.
    if (n < 4) then return 2 - ((n = 1) as integer)
    set factorCount to 2
    set sqrt to n ^ 0.5
    set limit to sqrt div 1
    if (limit = sqrt) then
        set factorCount to 3
        set limit to limit - 1
    end if
    repeat with i from 2 to limit
        if (n mod i = 0) then set factorCount to factorCount + 2
    end repeat
    
    return factorCount
end countFactors

on join(lst, delim)
    set astid to AppleScript's text item delimiters
    set AppleScript's text item delimiters to delim
    set txt to lst as text
    set AppleScript's text item delimiters to astid
    return txt
end join

on task()
    set limit to 100000 - 1
    set nWidth to (count (limit as text))
    set inset to "        "'s text 1 thru (nWidth - 1)
    set output to {"Positive integers < " & (limit + 1) & " whose divisor count is an odd prime:"}
    set row to {}
    set counter to 0
    repeat with sqrt from 2 to (limit ^ 0.5 div 1)
        set n to sqrt * sqrt
        if (countFactors(countFactors(n)) = 2) then
            set counter to counter + 1
            set row's end to (inset & n)'s text -nWidth thru -1
            if ((count row) = 10) then
                set output's end to join(row, "  ")
                set row to {}
            end if
        end if
    end repeat
    if (row ≠ {}) then set output's end to join(row, "  ")
    set output's end to linefeed & counter & " such integers"
    
    return join(output, linefeed)
end task

task()

Output:

"Positive integers < 100000 whose divisor count is an odd prime:
    4      9     16     25     49     64     81    121    169    289
  361    529    625    729    841    961   1024   1369   1681   1849
 2209   2401   2809   3481   3721   4096   4489   5041   5329   6241
 6889   7921   9409  10201  10609  11449  11881  12769  14641  15625
16129  17161  18769  19321  22201  22801  24649  26569  27889  28561
29929  32041  32761  36481  37249  38809  39601  44521  49729  51529
52441  54289  57121  58081  59049  63001  65536  66049  69169  72361
73441  76729  78961  80089  83521  85849  94249  96721  97969

79 such integers"

Arturo

numsWithPrimeNofDivisors: select 1..100000 'x [
    nofDivisors: size factors x
    and? [prime? nofDivisors]
         [nofDivisors <> 2]
]

loop split.every: 5 numsWithPrimeNofDivisors 'x ->
    print map x 's -> pad to :string s 6

Output:

     4      9     16     25     49 
    64     81    121    169    289 
   361    529    625    729    841 
   961   1024   1369   1681   1849 
  2209   2401   2809   3481   3721 
  4096   4489   5041   5329   6241 
  6889   7921   9409  10201  10609 
 11449  11881  12769  14641  15625 
 16129  17161  18769  19321  22201 
 22801  24649  26569  27889  28561 
 29929  32041  32761  36481  37249 
 38809  39601  44521  49729  51529 
 52441  54289  57121  58081  59049 
 63001  65536  66049  69169  72361 
 73441  76729  78961  80089  83521 
 85849  94249  96721  97969

AWK

# syntax: GAWK -f NUMBERS_WHOSE_COUNT_OF_DIVISORS_IS_PRIME.AWK
BEGIN {
    start = 2
    stop = 99999
    stop2 = 999
    for (i=start; i*i<=stop; i++) {
      n = count_divisors(i*i)
      if (n>2 && is_prime(n)) {
        printf("%6d%1s",i*i,++count%10?"":"\n")
        if (i*i <= stop2) {
          count2++
        }
      }
    }
    printf("\nNumbers with odd prime divisor counts %d-%d: %d\n",start,stop2,count2)
    printf("Numbers with odd prime divisor counts %d-%d: %d\n",start,stop,count)
    exit(0)
}
function count_divisors(n,  count,i) {
    for (i=1; i*i<=n; i++) {
      if (n % i == 0) {
        count += (i == n / i) ? 1 : 2
      }
    }
    return(count)
}
function is_prime(x,  i) {
    if (x <= 1) {
      return(0)
    }
    for (i=2; i<=int(sqrt(x)); i++) {
      if (x % i == 0) {
        return(0)
      }
    }
    return(1)
}

Output:

     4      9     16     25     49     64     81    121    169    289
   361    529    625    729    841    961   1024   1369   1681   1849
  2209   2401   2809   3481   3721   4096   4489   5041   5329   6241
  6889   7921   9409  10201  10609  11449  11881  12769  14641  15625
 16129  17161  18769  19321  22201  22801  24649  26569  27889  28561
 29929  32041  32761  36481  37249  38809  39601  44521  49729  51529
 52441  54289  57121  58081  59049  63001  65536  66049  69169  72361
 73441  76729  78961  80089  83521  85849  94249  96721  97969
Numbers with odd prime divisor counts 2-999: 16
Numbers with odd prime divisor counts 2-99999: 79

BASIC

BASIC256

Translation of: FreeBASIC

function isPrime(v)
	if v < 2 then return False
	if (v mod 2) = 0 then return v = 2
	if (v mod 3) = 0 then return v = 3
	d = 5
	while d * d <= v
		if (v mod d) = 0 then return False else d = d + 2
	end while
	return True
end function

row = 0

print "Numbers which count of divisors is prime are:"

for n = 1 to 1000
	num = 0
	for m = 1 to n
		if (n mod m) = 0 then num = num + 1
	next m
	if isPrime(num) and num <> 2 then
		print ""; n; "  ";
		row = row + 1
	end if
next n

print
print "Found "; row; " numbers"
end

FreeBASIC

Function isPrime(Byval ValorEval As Integer) As Boolean
    If ValorEval < 2 Then Return False
    If ValorEval Mod 2 = 0 Then Return ValorEval = 2
    If ValorEval Mod 3 = 0 Then Return ValorEval = 3
    Dim d As Integer = 5
    While d * d <= ValorEval
        If ValorEval Mod d = 0 Then Return False Else d += 2
    Wend 
    Return True
End Function

Dim As Integer row = 0
Dim As Uinteger n, num, m

Print "Numbers which count of divisors is prime are:"

For n = 1 To 100000
    num = 0
    For m = 1 To n
        If n Mod m = 0 Then num += 1 : End If
    Next m
    If isPrime(num) And num <> 2 Then
        Print Using " ##### "; n;
        row += 1
        If row Mod 5 = 0 Then Print : End If
    End If
Next n

Print !"\n\nFound"; row; " numbers"
Sleep

Output:

Numbers which count of divisors is prime are:
     4      9     16     25     49
    64     81    121    169    289
   361    529    625    729    841
   961   1024   1369   1681   1849
  2209   2401   2809   3481   3721
  4096   4489   5041   5329   6241
  6889   7921   9409  10201  10609
 11449  11881  12769  14641  15625
 16129  17161  18769  19321  22201
 22801  24649  26569  27889  28561
 29929  32041  32761  36481  37249
 38809  39601  44521  49729  51529
 52441  54289  57121  58081  59049
 63001  65536  66049  69169  72361
 73441  76729  78961  80089  83521
 85849  94249  96721  97969

Found 79 numbers

QBasic

Works with: QBasic

Works with: QuickBasic

Translation of: FreeBASIC

row = 0

PRINT "Numbers which count of divisors is prime are:"

FOR n = 1 TO 1000
    num = 0
    FOR m = 1 TO n
        IF n MOD m = 0 THEN num = num + 1
    NEXT m
    IF isPrime(num) AND num <> 2 THEN
        PRINT USING " ##### "; n;
        row = row + 1
        IF row MOD 5 = 0 THEN PRINT
    END IF
NEXT n

PRINT : PRINT "Found"; row; " numbers"
Sleep
END

FUNCTION isPrime (ValorEval)
    IF ValorEval < 2 THEN isPrime = False
    IF ValorEval MOD 2 = 0 THEN isPrime = 2
    IF ValorEval MOD 3 = 0 THEN isPrime = 3
    d = 5
    WHILE d * d <= ValorEval
        IF ValorEval MOD d = 0 THEN isPrime = False ELSE d = d + 2
    WEND
    isPrime = True
END FUNCTION

PureBasic

Translation of: FreeBASIC

Procedure isPrime(v.i)
  If     v <= 1    : ProcedureReturn #False
  ElseIf v < 4     : ProcedureReturn #True
  ElseIf v % 2 = 0 : ProcedureReturn #False
  ElseIf v < 9     : ProcedureReturn #True
  ElseIf v % 3 = 0 : ProcedureReturn #False
  Else
    Protected r = Round(Sqr(v), #PB_Round_Down)
    Protected f = 5
    While f <= r
      If v % f = 0 Or v % (f + 2) = 0
        ProcedureReturn #False
      EndIf
      f + 6
    Wend
  EndIf
  ProcedureReturn #True
EndProcedure

OpenConsole()
fila.i = 0

PrintN("Numbers which count of divisors is prime are:")

For n.i = 1 To 100000
  num.i = 0
  For m.i = 1 To n
    If Mod(n, m) = 0 
      num + 1
    EndIf
  Next
  If isPrime(num) And num <> 2
    Print("  " + Str(n) + "  ")
    fila + 1
    If Mod(fila, 5) = 0
      PrintN("")
    EndIf
  EndIf
Next

PrintN(#CRLF$ + "Found " + Str(fila) + " numbers")

Input()
CloseConsole()
End

Yabasic

Translation of: FreeBASIC

row = 0

print "Numbers which count of divisors is prime are:"

for n = 1 to 1000
    num = 0
    for m = 1 to n
        if mod(n, m) = 0 then num = num + 1 : fi
    next m
    if isPrime(num) and num <> 2 then
		print n using "#####",
        row = row + 1
        if mod(row, 5) = 0 then print : fi
    end if
next n

print "\n\nFound ", row, " numbers"
end

sub isPrime(v)
    if v < 2 then return False : fi
    if mod(v, 2) = 0 then return v = 2 : fi
    if mod(v, 3) = 0 then return v = 3 : fi
    d = 5
    while d * d <= v
        if mod(v, d) = 0 then return False else d = d + 2 : fi
    end while
    return True
end sub

Output:

Igual que la entrada de FreeBASIC.

C++

#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>

int divisor_count(int n) {
    int total = 1;
    for (; (n & 1) == 0; n >>= 1)
        ++total;
    for (int p = 3; p * p <= n; p += 2) {
        int count = 1;
        for (; n % p == 0; n /= p)
            ++count;
        total *= count;
    }
    if (n > 1)
        total *= 2;
    return total;
}

bool is_prime(int n) {
    if (n < 2)
        return false;
    if (n % 2 == 0)
        return n == 2;
    if (n % 3 == 0)
        return n == 3;
    for (int p = 5; p * p <= n; p += 4) {
        if (n % p == 0)
            return false;
        p += 2;
        if (n % p == 0)
            return false;
    }
    return true;
}

int main(int argc, char** argv) {
    int limit = 1000;
    switch (argc) {
    case 1:
        break;
    case 2:
        limit = std::strtol(argv[1], nullptr, 10);
        if (limit <= 0) {
            std::cerr << "Invalid limit\n";
            return EXIT_FAILURE;
        }
        break;
    default:
        std::cerr << "usage: " << argv[0] << " [limit]\n";
        return EXIT_FAILURE;
    }
    int width = static_cast<int>(std::ceil(std::log10(limit)));
    int count = 0;
    for (int i = 1;; ++i) {
        int n = i * i;
        if (n >= limit)
            break;
        int divisors = divisor_count(n);
        if (divisors != 2 && is_prime(divisors))
            std::cout << std::setw(width) << n << (++count % 10 == 0 ? '\n' : ' ');
    }
    std::cout << "\nCount: " << count << '\n';
    return EXIT_SUCCESS;
}

Output:

Default input:

  4   9  16  25  49  64  81 121 169 289
361 529 625 729 841 961 
Count: 16

Stretch goal:

    4     9    16    25    49    64    81   121   169   289
  361   529   625   729   841   961  1024  1369  1681  1849
 2209  2401  2809  3481  3721  4096  4489  5041  5329  6241
 6889  7921  9409 10201 10609 11449 11881 12769 14641 15625
16129 17161 18769 19321 22201 22801 24649 26569 27889 28561
29929 32041 32761 36481 37249 38809 39601 44521 49729 51529
52441 54289 57121 58081 59049 63001 65536 66049 69169 72361
73441 76729 78961 80089 83521 85849 94249 96721 97969 
Count: 79

CLU

% Find the amount of divisors for 1..N
div_counts = proc (n: int) returns (sequence[int])
    divs: array[int] := array[int]$fill(1,n,1)
    for d: int in int$from_to(2, n) do
        for m: int in int$from_to_by(d, n, d) do
            divs[m] := divs[m] + 1
        end
    end
    return(sequence[int]$a2s(divs))
end div_counts

% Find maximum element of sequence
max = proc (seq: sequence[int]) returns (int)
    n: int := 0
    for i: int in sequence[int]$elements(seq) do
        if i>n then n:=i end
    end
    return(n)
end max

% Sieve primes up to N
sieve = proc (n: int) returns (sequence[bool])
    prime: array[bool] := array[bool]$fill(1,n,true)
    prime[1] := false
    for p: int in int$from_to(2, n/2) do
        for c: int in int$from_to_by(p*p, n, p) do
            prime[c] := false
        end
    end
    return(sequence[bool]$a2s(prime))
end sieve

start_up = proc ()
    MAX_N = 100000
    po: stream := stream$primary_output()
    div_count: sequence[int] := div_counts(MAX_N)
    prime: sequence[bool] := sieve(max(div_count))
    count: int := 0
    for i: int in int$from_to(1, MAX_N) do
        dc: int := div_count[i]
        if dc ~= 2 cand prime[dc] then
            stream$putright(po, int$unparse(i), 8)
            count := count + 1
            if count//10 = 0 then stream$putl(po, "") end
        end
    end
    stream$putl(po, "\nFound " || int$unparse(count) || " numbers.")
end start_up

Output:

       4       9      16      25      49      64      81     121     169     289
     361     529     625     729     841     961    1024    1369    1681    1849
    2209    2401    2809    3481    3721    4096    4489    5041    5329    6241
    6889    7921    9409   10201   10609   11449   11881   12769   14641   15625
   16129   17161   18769   19321   22201   22801   24649   26569   27889   28561
   29929   32041   32761   36481   37249   38809   39601   44521   49729   51529
   52441   54289   57121   58081   59049   63001   65536   66049   69169   72361
   73441   76729   78961   80089   83521   85849   94249   96721   97969
Found 79 numbers.

Delphi

Works with: Delphi version 6.0

Library: SysUtils,StdCtrls

procedure GetUniqueFactors(N: integer; var IA: TIntegerDynArray);
{Get unique factors of a number}
var I: integer;
begin
SetLength(IA,1);
for I:=2 to N do
 if (N mod I)=0 then
	begin
	SetLength(IA,Length(IA)+1);
	IA[High(IA)]:=I;
	end;
end;


procedure ShowCountPrimes(Memo: TMemo);
{Count the number of Unique factors that are prime}
var I,C,Cnt: integer;
var IA: TIntegerDynArray;
var S: string;
begin
Cnt:=0; S:='';
for I:=1 to 100000-1 do
	begin
	GetUniqueFactors(I,IA);
	C:=Length(IA);
	if (C=2) or (not IsPrime(C)) then continue;
	Inc(Cnt);
	S:=S+Format('%8D',[I]);
	If (Cnt mod 5)=0 then S:=S+CRLF;
	end;
Memo.Lines.Add('Count='+IntToStr(Cnt));
Memo.Lines.Add(S);
end;

Output:

Count=79
       4       9      16      25      49
      64      81     121     169     289
     361     529     625     729     841
     961    1024    1369    1681    1849
    2209    2401    2809    3481    3721
    4096    4489    5041    5329    6241
    6889    7921    9409   10201   10609
   11449   11881   12769   14641   15625
   16129   17161   18769   19321   22201
   22801   24649   26569   27889   28561
   29929   32041   32761   36481   37249
   38809   39601   44521   49729   51529
   52441   54289   57121   58081   59049
   63001   65536   66049   69169   72361
   73441   76729   78961   80089   83521
   85849   94249   96721   97969
Elapsed Time: 17.921 Sec.

EasyLang

fastfunc isprim num .
   if num mod 2 = 0
      return 0
   .
   i = 3
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 2
   .
   return 1
.
for n to 999
   cnt = 0
   for m to n
      cnt += if n mod m = 0
   .
   if cnt > 2 and isprim cnt = 1
      write n & " "
   .
.

Output:

4 9 16 25 49 64 81 121 169 289 361 529 625 729 841 961 

F#

This task uses Extensible Prime Generator (F#)

// Numbers whose divisor count is prime. Nigel Galloway: July 13th., 2021
primes64()|>Seq.takeWhile(fun n->n*n<100000L)|>Seq.collect(fun n->primes32()|>Seq.skip 1|>Seq.map(fun g->pown n (g-1))|>Seq.takeWhile((>)100000L))|>Seq.sort|>Seq.iter(printf "%d "); printfn ""

Output:

4 9 16 25 49 64 81 121 169 289 361 529 625 729 841 961 1024 1369 1681 1849 2209 2401 2809 3481 3721 4096 4489 5041 5329 6241 6889 7921 9409 10201 10609 11449 11881 12769 14641 15625 16129 17161 18769 19321 22201 22801 24649 26569 27889 28561 29929 32041 32761 36481 37249 38809 39601 44521 49729 51529 52441 54289 57121 58081 59049 63001 65536 66049 69169 72361 73441 76729 78961 80089 83521 85849 94249 96721 97969

Factor

Works with: Factor version 0.99 2021-06-02

USING: formatting grouping io kernel math math.primes
math.primes.factors math.ranges sequences sequences.extras ;
FROM: math.extras => integer-sqrt ;

: odd-prime? ( n -- ? ) dup 2 = [ drop f ] [ prime? ] if ;

: pdc-upto ( n -- seq )
    integer-sqrt [1,b]
    [ sq ] [ divisors length odd-prime? ] map-filter ;

100,000 pdc-upto 10 group [ [ "%-8d" printf ] each nl ] each

Output:

4       9       16      25      49      64      81      121     169     289     
361     529     625     729     841     961     1024    1369    1681    1849    
2209    2401    2809    3481    3721    4096    4489    5041    5329    6241    
6889    7921    9409    10201   10609   11449   11881   12769   14641   15625   
16129   17161   18769   19321   22201   22801   24649   26569   27889   28561   
29929   32041   32761   36481   37249   38809   39601   44521   49729   51529   
52441   54289   57121   58081   59049   63001   65536   66049   69169   72361   
73441   76729   78961   80089   83521   85849   94249   96721   97969   

Go

Library: Go-rcu

package main

import (
    "fmt"
    "rcu"
)

func countDivisors(n int) int {
    count := 0
    i := 1
    k := 1
    if n%2 == 1 {
        k = 2
    }
    for ; i*i <= n; i += k {
        if n%i == 0 {
            count++
            j := n / i
            if j != i {
                count++
            }
        }
    }
    return count
}

func main() {
    const limit = 1e5
    var results []int
    for i := 2; i * i < limit; i++ {
        n := countDivisors(i * i)
        if n > 2 && rcu.IsPrime(n) {
            results = append(results, i * i)
        }
    }
    climit := rcu.Commatize(limit)
    fmt.Printf("Positive integers under %7s whose number of divisors is an odd prime:\n", climit)
    under1000 := 0
    for i, n := range results {
        fmt.Printf("%7s", rcu.Commatize(n))
        if (i+1)%10 == 0 {
            fmt.Println()
        }
        if n < 1000 {
            under1000++
        }
    }
    fmt.Printf("\n\nFound %d such integers (%d under 1,000).\n", len(results), under1000)
}

Output:

Positive integers under 100,000 whose number of divisors is an odd prime:
      4      9     16     25     49     64     81    121    169    289
    361    529    625    729    841    961  1,024  1,369  1,681  1,849
  2,209  2,401  2,809  3,481  3,721  4,096  4,489  5,041  5,329  6,241
  6,889  7,921  9,409 10,201 10,609 11,449 11,881 12,769 14,641 15,625
 16,129 17,161 18,769 19,321 22,201 22,801 24,649 26,569 27,889 28,561
 29,929 32,041 32,761 36,481 37,249 38,809 39,601 44,521 49,729 51,529
 52,441 54,289 57,121 58,081 59,049 63,001 65,536 66,049 69,169 72,361
 73,441 76,729 78,961 80,089 83,521 85,849 94,249 96,721 97,969

Found 79 such integers (16 under 1,000).

J

To find the divisors of a number we can use any of the implementations defined at in the Factors of an integer task: foi, factrs, factrst, factorsOfNumber or factors. Here, we will use factors as it's the most efficient:

   (#~ (2&< * 1&p:)@#@factors"0) }. i. 100000
4 9 16 25 49 64 81 121 169 289 361 529 625 729 841 961 1024 1369 1681 1849 2209 2401 2809 3481 3721 4096 4489 5041 5329 6241 6889 7921 9409 10201 10609 11449 11881 12769 14641 15625 16129 17161 18769 19321 22201 22801 24649 26569 27889 28561 29929 32041 32761 36481 37249 38809 39601 44521 49729 51529 52441 54289 57121 58081 59049 63001 65536 66049 69169 72361 73441 76729 78961 80089 83521 85849 94249 96721 97969

(Note that we first conditioned J's environment to show longer lines, using 9!:37(10*9!:36).)

jq

Works with: jq

Works with gojq, the Go implementation of jq

For a suitable definition of `is_prime`, see Erdős-primes#jq.

def add(s): reduce s as $x (null; .+$x);

def count_divisors:
  add(
    if . == 1 then 1
    else . as $n
    | label $out
    | range(1; $n) as $i
    | ($i * $i) as $i2
    | if $i2 > $n then break $out
      else if $i2 == $n
           then 1
           elif ($n % $i) == 0
           then 2
           else empty
  	   end
      end
    end);

1000, 100000
| "\nn with odd prime divisor counts, 1 < n < \(.):",
  (range(1;.) | select(count_divisors | (. > 2 and is_prime)))

Output:

n with odd prime divisor counts, 1 < n < 1000:
4
9
16
25
49
64
81
121
169
289
361
529
625
729
841
961

n with odd prime divisor counts, 1 < n < 100000:
4
9
....
85849
94249
96721
97969

Julia

using Primes

ispdc(n) = (ndivs = prod(collect(values(factor(n))).+ 1); ndivs > 2 && isprime(ndivs))

foreach(p -> print(rpad(p[2], 8), p[1] % 10 == 0 ? "\n" : ""), enumerate(filter(ispdc, 1:100000)))

Output:

4       9       16      25      49      64      81      121     169     289     
361     529     625     729     841     961     1024    1369    1681    1849
2209    2401    2809    3481    3721    4096    4489    5041    5329    6241
6889    7921    9409    10201   10609   11449   11881   12769   14641   15625
16129   17161   18769   19321   22201   22801   24649   26569   27889   28561
29929   32041   32761   36481   37249   38809   39601   44521   49729   51529
52441   54289   57121   58081   59049   63001   65536   66049   69169   72361
73441   76729   78961   80089   83521   85849   94249   96721   97969

Mathematica / Wolfram Language

max = 100000;
maxPrime = NextPrime[Sqrt@max, -1];
maxPower = NextPrime[Log[2, max], -1];
base = NestWhileList[NextPrime, 2, # < maxPrime &];
g = NestWhileList[NextPrime, 3, # < maxPower &] - 1;
ans = Sort@Select[Flatten@Table[base^n, {n, g}], # < max &];
Labeled[Partition[Select[ans, # < 1000 &], UpTo[8]] // 
  TableForm, "Numbers up to 1000 with prime divisor counts:", Top]
Labeled[Partition[ans, UpTo[8]] // 
  TableForm, "Numbers up to 100,000 with prime divisor counts:", Top]

Output:

Numbers up to 1000 with prime divisor counts:
4	9	16	25	49	64	81	121
169	289	361	529	625	729	841	961
Numbers up to 100,000 with prime divisor counts:
4   	9   	16  	25  	49  	64  	81  	121
169 	289 	361 	529 	625 	729 	841 	961
1024	1369	1681	1849	2209	2401	2809	3481
3721	4096	4489	5041	5329	6241	6889	7921
9409	10201	10609	11449	11881	12769	14641	15625
16129	17161	18769	19321	22201	22801	24649	26569
27889	28561	29929	32041	32761	36481	37249	38809
39601	44521	49729	51529	52441	54289	57121	58081
59049	63001	66049	69169	72361	73441	76729	78961
80089	83521	85849	94249	96721	97969				

Nim

Checking only divisors of squares (see discussion).

import math, sequtils, strformat, strutils


func divCount(n: Positive): int =
  var n = n
  for d in 1..n:
    if d * d > n: break
    if n mod d == 0:
      inc result
      if n div d != d:
        inc result


func isOddPrime(n: Positive): bool =
  if n < 3 or n mod 2 == 0: return false
  if n mod 3 == 0: return n == 3
  var d = 5
  while d <= sqrt(n.toFloat).int:
    if n mod d == 0: return false
    inc d, 2
    if n mod d == 0: return false
    inc d, 4
  result = true


iterator numWithOddPrimeDivisorCount(lim: Positive): int =
  for k in 1..sqrt(lim.toFloat).int:
    let n = k * k
    if n.divCount().isOddPrime():
      yield n


var list = toSeq(numWithOddPrimeDivisorCount(1000))

echo &"Found {list.len} numbers between 1 and 999 whose number of divisors is an odd prime:"
echo list.join(" ")
echo()

list = toSeq(numWithOddPrimeDivisorCount(100_000))
echo &"Found {list.len} numbers between 1 and 99_999 whose number of divisors is an odd prime:"
for i, n in list:
  stdout.write &"{n:5}", if (i + 1) mod 10 == 0: '\n' else: ' '
echo()

Output:

Found 16 numbers between 1 and 999 whose number of divisors is an odd prime:
4 9 16 25 49 64 81 121 169 289 361 529 625 729 841 961

Found 79 numbers between 1 and 99_999 whose number of divisors is an odd prime:
    4     9    16    25    49    64    81   121   169   289
  361   529   625   729   841   961  1024  1369  1681  1849
 2209  2401  2809  3481  3721  4096  4489  5041  5329  6241
 6889  7921  9409 10201 10609 11449 11881 12769 14641 15625
16129 17161 18769 19321 22201 22801 24649 26569 27889 28561
29929 32041 32761 36481 37249 38809 39601 44521 49729 51529
52441 54289 57121 58081 59049 63001 65536 66049 69169 72361
73441 76729 78961 80089 83521 85849 94249 96721 97969

PARI/GP

Translation of: Julia

ispdc(n) = {
    local(factors, ndivs);
    factors = factor(n);
    ndivs = numdiv(n);
    ndivs > 2 && isprime(ndivs)
};

{
cnt=0;
for (i=1, 100000,
    if (ispdc(i),
        print1(i, " ");
        cnt=cnt+1;
        if(cnt%10==0,print(""));
    );
);
}

Output:

4 9 16 25 49 64 81 121 169 289 
361 529 625 729 841 961 1024 1369 1681 1849 
2209 2401 2809 3481 3721 4096 4489 5041 5329 6241 
6889 7921 9409 10201 10609 11449 11881 12769 14641 15625 
16129 17161 18769 19321 22201 22801 24649 26569 27889 28561 
29929 32041 32761 36481 37249 38809 39601 44521 49729 51529 
52441 54289 57121 58081 59049 63001 65536 66049 69169 72361 
73441 76729 78961 80089 83521 85849 94249 96721 97969 

Pascal

program FacOfInteger;
{$IFDEF FPC}
//  {$R+,O+} //debuging purpose
  {$MODE DELPHI}
  {$Optimization ON,ALL}
{$ELSE}
   {$APPTYPE CONSOLE}
{$ENDIF}
uses
  sysutils;
//#############################################################################
//Prime decomposition  
type
  tPot = record
           potSoD : Uint64;
           potPrim,
           potMax :Uint32;
         end;

  tprimeFac = record
                 pfPrims : array[0..13] of tPot;
                 pfSumOfDivs : Uint64;
                 pfCnt,
                 pfNum,
                 pfDivCnt: Uint32;
               end;
 
  tSmallPrimes = array[0..6541] of Word;
  tItem     = NativeUint;
  tDivisors = array of tItem;
  tpDivisor = pNativeUint;
var
  SmallPrimes: tSmallPrimes;

 
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 InitSmallPrimes;
var
  pr,testPr,j,maxprimidx,delta: Uint32;
  isPrime : boolean;
Begin
  SmallPrimes[0] := 2;
  SmallPrimes[1] := 3;  
  delta := 2;
  maxprimidx := 1;  
  pr := 5;
  repeat
    isprime := true;
    j := 0;
    repeat
      testPr := SmallPrimes[j];
      IF testPr*testPr > pr then
        break;
      If pr mod testPr = 0 then
      Begin
        isprime := false;
        break;
      end;
      inc(j);
    until false;
 
    if isprime then
    Begin
      inc(maxprimidx);
      SmallPrimes[maxprimidx]:= pr;
    end;
    inc(pr,delta);
    delta := 2+4-delta;
  until pr > 1 shl 16 -1;
end;

function isPrime(n:Uint32):boolean;
var
  pr,idx: NativeInt;
begin
  result := n in [2,3];
  if NOT(result) AND (n>4) AND (n AND 1 <> 0 ) then
  begin
    idx := 1;
    repeat
      pr := SmallPrimes[idx];
      result := (n mod pr) <>0;
      inc(idx);
    until NOT(result) or (sqr(pr)>n) or (idx > High(SmallPrimes)); 
  end;
end; 

procedure PrimeFacOut(const primeDecomp:tprimeFac;proper:Boolean=true);
var
  i,k : Int32;
begin
  with primeDecomp do
  Begin
    write(pfNum,' = ');
    k := pfCnt-1;
    For i := 0 to k-1 do
      with pfPrims[i] do
        If potMax = 1 then
          write(potPrim,'*')
        else
          write(potPrim,'^',potMax,'*');
    with pfPrims[k] do
      If potMax = 1 then
        write(potPrim)
      else
        write(potPrim,'^',potMax);
    if proper then
      writeln(' got ',pfDivCnt-1,' proper divisors with sum : ',pfSumOfDivs-pfNum)
    else
      writeln(' got ',pfDivCnt,' divisors with sum : ',pfSumOfDivs);
  end;
end;

procedure PrimeDecomposition(var res:tprimeFac;n:Uint32);
var
  DivSum,fac:Uint64;
  i,pr,cnt,DivCnt,quot{to minimize divisions} : NativeUint;
Begin
  if SmallPrimes[0] <> 2 then
    InitSmallPrimes;
  res.pfNum := n;
  cnt := 0;
  DivCnt := 1;
  DivSum := 1;
  i := 0;
  if n <= 1 then
  Begin
    with res.pfPrims[0] do
    Begin
      potPrim := n;
      potMax  := 1;
    end;
    cnt := 1;
  end
  else
  repeat
    pr := SmallPrimes[i];
    IF pr*pr>n then
      Break;
 
    quot := n div pr;
    IF pr*quot = n then
      with res do
      Begin
        with pfPrims[Cnt] do
        Begin
          potPrim := pr;
          potMax := 0;
          fac := pr;
          repeat
            n := quot;
            quot := quot div pr;
            inc(potMax);
            fac *= pr;
          until pr*quot <> n;
          DivCnt *= (potMax+1);
          DivSum *= (fac-1)DIV (pr-1);          
        end;
        inc(Cnt);
      end;
     inc(i);
  until false;
  //a big prime left over?
  IF n > 1 then
    with res do
    Begin
      with pfPrims[Cnt] do
      Begin
        potPrim := n;
        potMax := 1;
      end;
      inc(Cnt);
      DivCnt *= 2;
      DivSum *= n+1;      
    end;
  with res do
  Begin  
   pfCnt:= cnt;
   pfDivCnt := DivCnt;
   pfSumOfDivs := DivSum;
  end; 
end;

function isAbundant(const pD:tprimeFac):boolean;inline;
begin
  with pd do
    result := pfSumOfDivs-pfNum > pfNum;
end;

function DivCount(const pD:tprimeFac):NativeUInt;inline;
begin
  result := pD.pfDivCnt;
end;
 
function SumOfDiv(const primeDecomp:tprimeFac):NativeUInt;inline;
begin
  result := primeDecomp.pfSumOfDivs;
end;
 
procedure GetDivs(var pD:tprimeFac;var Divs:tDivisors);
var
  pDivs : tpDivisor;
  i,len,j,l,p,pPot,k: NativeInt;
Begin
  i := DivCount(pD);
  IF i > Length(Divs) then
    setlength(Divs,i);
  pDivs := @Divs[0];
  pDivs[0] := 1;
  len := 1;
  l := len;
  For i := 0 to pD.pfCnt-1 do
    with pD.pfPrims[i] do
    Begin
      //Multiply every divisor before with the new primefactors
      //and append them to the list
      k := potMax-1;
      p := potPrim;
      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;
  //Sort. Insertsort much faster than QuickSort in this special case
  InsertSort(pDivs,0,len-1);
end;
 
Function GetDivisors(var pD:tprimeFac;n:Uint32;var Divs:tDivisors):Int32;
var
  i:Int32;
Begin
  if pD.pfNum <> n then
    PrimeDecomposition(pD,n);
  i := DivCount(pD);
  IF i > Length(Divs) then
    setlength(Divs,i+1);
  GetDivs(pD,Divs);
  result := DivCount(pD);
end;
 
procedure AllFacsOut(var pD:tprimeFac;n: Uint32;Divs:tDivisors;proper:boolean=true);
var
  k,j: Int32;
Begin
  k := GetDivisors(pD,n,Divs)-1;// zero based
  PrimeFacOut(pD,proper);
  IF proper then
    dec(k);
  IF k > 0 then
  Begin
    For j := 0 to k-1 do
      write(Divs[j],',');
    writeln(Divs[k]);
  end;
end;
//Prime decomposition  
//#############################################################################
procedure SpeedTest(var pD: tprimeFac;Limit:Uint32);
var
  Ticks : Int64;
  number,numSqr,Cnt: UInt32;
Begin
  Ticks := GetTickCount64;
  Cnt := 0;
  number := 1;
  numSqr:=1;
  repeat
    number += 1;
    numSqr := sqr(number);
    PrimeDecomposition(pD,numSqr);
    IF DivCount(pD)>2 then
      if isPrime(DivCount(pD)) then
        inc(cnt);//writeln(number:5,numSqr:10,DivCount(pD):5);
  until numSqr>= Limit;
  writeln('SpeedTest ',(GetTickCount64-Ticks)/1000:0:3,' secs for 1..',Limit,' found ',Cnt);
  writeln;
end;
 
var
  pD: tprimeFac;
  Divisors : tDivisors;
  numroot,num,cnt : Uint32;
BEGIN
  InitSmallPrimes;
  setlength(Divisors,1);
  
  write('':4);
  for cnt := 1 to 10 do
    write(cnt:7);
  writeln;  
    
  cnt := 0;
  write(cnt:3,':');
  For numroot := 2 to 1000 do
  begin
    num := sqr(numroot);
    PrimeDecomposition(pD,num);
    IF DivCount(pD)>2 then
      if isPrime(DivCount(pD)) then
      begin
        write(num:7);
        inc(cnt);      
        if cnt MOD 10 =0 then
        Begin
          writeln;write(cnt:3,':');
        end;  
      end;    
  end;      
  if cnt MOD 8 <>0 then
    writeln;
  writeln;
  SpeedTest(pD,4000*1000*1000);
END.

Output:

          1      2      3      4      5      6      7      8      9     10
  0:      4      9     16     25     49     64     81    121    169    289
 10:    361    529    625    729    841    961   1024   1369   1681   1849
 20:   2209   2401   2809   3481   3721   4096   4489   5041   5329   6241
 30:   6889   7921   9409  10201  10609  11449  11881  12769  14641  15625
 40:  16129  17161  18769  19321  22201  22801  24649  26569  27889  28561
 50:  29929  32041  32761  36481  37249  38809  39601  44521  49729  51529
 60:  52441  54289  57121  58081  59049  63001  65536  66049  69169  72361
 70:  73441  76729  78961  80089  83521  85849  94249  96721  97969 100489
 80: 109561 113569 117649 120409 121801 124609 128881 130321 134689 139129
 90: 143641 146689 151321 157609 160801 167281 175561 177241 185761 187489
100: 192721 196249 201601 208849 212521 214369 218089 229441 237169 241081
110: 249001 253009 259081 262144 271441 273529 279841 292681 299209 310249
120: 316969 323761 326041 332929 344569 351649 358801 361201 368449 375769
130: 380689 383161 398161 410881 413449 418609 426409 434281 436921 452929
140: 458329 466489 477481 491401 502681 516961 528529 531441 537289 546121
150: 552049 564001 573049 579121 591361 597529 619369 635209 654481 657721
160: 674041 677329 683929 687241 703921 707281 727609 734449 737881 744769
170: 769129 776161 779689 786769 822649 829921 844561 863041 877969 885481
180: 896809 908209 923521 935089 942841 954529 966289 982081 994009

SpeedTest 0.230 secs for 1..4000000000 found 6417

Perl

Library: ntheory

use strict;
use warnings;
use ntheory <is_prime divisors>;

push @matches, $_**2 for grep { is_prime divisors $_**2 } 1..int sqrt 1e5;
print @matches . " matching:\n" . (sprintf "@{['%6d' x @matches]}", @matches) =~ s/(.{72})/$1\n/gr;

Output:

79 matching:
     4     9    16    25    49    64    81   121   169   289   361   529
   625   729   841   961  1024  1369  1681  1849  2209  2401  2809  3481
  3721  4096  4489  5041  5329  6241  6889  7921  9409 10201 10609 11449
 11881 12769 14641 15625 16129 17161 18769 19321 22201 22801 24649 26569
 27889 28561 29929 32041 32761 36481 37249 38809 39601 44521 49729 51529
 52441 54289 57121 58081 59049 63001 65536 66049 69169 72361 73441 76729
 78961 80089 83521 85849 94249 96721 97969

Phix

atom t0 = time()
function pd(integer n) n = length(factors(n,1)) return n!=2 and is_prime(n) end function
for n in {1e3,1e5} do
    sequence res = filter(tagset(n),pd)
    printf(1,"%d < %,d found: %v\n",{length(res),n,shorten(res,"",5)})
end for
?elapsed(time()-t0)

Output:

16 < 1,000 found: {4,9,16,25,49,"...",529,625,729,841,961}
79 < 100,000 found: {4,9,16,25,49,"...",83521,85849,94249,96721,97969}
"0.4s"

smarter, faster

As per Nigel's analysis on the talk page, valid numbers must be of the form a^(b-1) where a is prime and b is an odd prime, with no further checking required.

atom t1 = time()
for n in {1e3,1e5,4e9} do
    sequence r = {}
    for a in get_primes_le(floor(sqrt(n))) do
        integer b = 2
        while true do
            atom c = power(a,get_prime(b)-1)
            if c>n then exit end if
            r &= c
            b += 1
        end while
    end for
    r = sort(deep_copy(r))
    printf(1,"%d < %,d found: %v\n",{length(r),n,shorten(r,"",5)})
end for
?elapsed(time()-t1)

Output:

16 < 1,000 found: {4,9,16,25,49,"...",529,625,729,841,961}
79 < 100,000 found: {4,9,16,25,49,"...",83521,85849,94249,96721,97969}
6417 < 4,000,000,000 found: {4,9,16,25,49,"...",3991586041.0,3993860809.0,3994113601.0,3995630521.0,3999424081.0}
"0.0s"

Quackery

Noting that only perfect squares can fit the criteria (see discussion).

factors and isqrt are defined at Factors of an integer#Quackery.

  [ dip number$
    over size -
    space swap of
    swap join echo$ ] is recho ( n n -->   )

  [ []
    100000 isqrt drop times
      [ i^ dup * factors size
        dup 2 = iff drop done
        factors size 2 = while
        i^ dup * join ] ]

  [ dup [] != while
    10 split swap
    witheach [ 6 recho ]
    cr again ]
   drop

Output:

     4     9    16    25    49    64    81   121   169   289
   361   529   625   729   841   961  1024  1369  1681  1849
  2209  2401  2809  3481  3721  4096  4489  5041  5329  6241
  6889  7921  9409 10201 10609 11449 11881 12769 14641 15625
 16129 17161 18769 19321 22201 22801 24649 26569 27889 28561
 29929 32041 32761 36481 37249 38809 39601 44521 49729 51529
 52441 54289 57121 58081 59049 63001 65536 66049 69169 72361
 73441 76729 78961 80089 83521 85849 94249 96721 97969

Raku

use Prime::Factor;

my $ceiling = ceiling sqrt 1e5;

say display :10cols, :fmt('%6d'), (^$ceiling)».² .grep: { .&divisors.is-prime };

sub display ($list, :$cols = 10, :$fmt = '%6d', :$title = "{+$list} matching:\n" )   {
    cache $list;
    $title ~ $list.batch($cols)».fmt($fmt).join: "\n"
}

Output:

79 matching:
     4      9     16     25     49     64     81    121    169    289
   361    529    625    729    841    961   1024   1369   1681   1849
  2209   2401   2809   3481   3721   4096   4489   5041   5329   6241
  6889   7921   9409  10201  10609  11449  11881  12769  14641  15625
 16129  17161  18769  19321  22201  22801  24649  26569  27889  28561
 29929  32041  32761  36481  37249  38809  39601  44521  49729  51529
 52441  54289  57121  58081  59049  63001  65536  66049  69169  72361
 73441  76729  78961  80089  83521  85849  94249  96721  97969

REXX

/*REXX pgm finds positive integers N whose # of divisors is prime (& ¬=2), where N<1000.*/
parse arg hi cols .                              /*obtain optional arguments from the CL*/
if   hi=='' |   hi==","  then   hi= 1000         /*Not specified?  Then use the defaults*/
if cols=='' | cols==","  then cols=   10         /* "      "         "   "   "      "   */
call genP                                        /*build array of semaphores for primes.*/
w= 10                                            /*W:  the maximum width of any column. */
title= ' positive integers  N  whose number of divisors is prime (and not equal to 2), ' ,
       "where  N < "            commas(hi)
say ' index │'center(title, 1 + cols*(w+1)     )
say '───────┼'center(""   , 1 + cols*(w+1), '─')
finds= 0;                  idx= 1;            $=
         do j=2; jj= j*j;  if jj>=hi  then leave /*process positive square ints in range*/
         n= nDivs(jj);     if n==2  then iterate /*get number of divisors of composite J*/
                           if \!.n  then iterate /*Number divisors prime?  No, then skip*/
         finds= finds + 1                        /*bump the number of found numbers.    */
         $= $  right( commas(j),  w)             /*add a positive integer  ──►  $ list. */
         if finds//cols\==0 then iterate         /*have we populated a line of output?  */
         say center(idx, 7)'│'  substr($, 2); $= /*display what we have so far  (cols). */
         idx= idx + cols                         /*bump the  index  count for the output*/
         end   /*j*/                             /* [↑]   process a range of integers.  */

if $\==''  then say center(idx, 7)"│"  substr($, 2)  /*possible display residual output.*/
say '───────┴'center(""   , 1 + cols*(w+1), '─')
say
say 'Found '     commas(finds)      title
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg ?;  do jc=length(?)-3  to 1  by -3; ?=insert(',', ?, jc); end;  return ?
/*──────────────────────────────────────────────────────────────────────────────────────*/
nDivs: procedure; parse arg x;  d= 2;  if x==1  then return 1 /*handle special case of 1*/
       odd= x // 2                               /* [↓]  process EVEN or ODD ints.   ___*/
              do j=2+odd  by 1+odd  while j*j<x  /*divide by all the integers up to √ x */
              if x//j==0  then d= d + 2          /*÷?  Add two divisors to the total.   */
              end   /*j*/                        /* [↑]  %  ≡  integer division.        */
       if j*j==x  then  return  d + 1            /*Was X a square?  Then add 1 to total.*/
                        return  d                /*return the total.                    */
/*──────────────────────────────────────────────────────────────────────────────────────*/
genP:        @.1=2; @.2=3; @.3=5; @.4=7;  @.5=11 /*define some low primes.              */
      !.=0;  !.2=1; !.3=1; !.5=1; !.7=1;  !.11=1 /*   "     "   "    "     semaphores.  */
                           #=5;   s.#= @.# **2   /*number of primes so far;     prime². */
        do j=@.#+2  by 2  to hi-1                /*find odd primes from here on.        */
        parse var j '' -1 _; if     _==5  then iterate  /*J divisible by 5?  (right dig)*/
                             if j// 3==0  then iterate  /*"     "      " 3?             */
                             if j// 7==0  then iterate  /*"     "      " 7?             */
               do k=5  while s.k<=j              /* [↓]  divide by the known odd primes.*/
               if j // @.k == 0  then iterate j  /*Is  J ÷ X?  Then not prime.     ___  */
               end   /*k*/                       /* [↑]  only process numbers  ≤  √ J   */
        #= #+1;    @.#= j;    s.#= j*j;   !.j= 1 /*bump # of Ps; assign next P;  P²; P# */
        end          /*j*/;   return

output   when using the default inputs:

 index │        positive integers  N  whose number of divisors is prime (and not equal to 2),  where  N <  1,000
───────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1   │          4          9         16         25         49         64         81        121        169        289
  11   │        361        529        625        729        841        961
───────┴───────────────────────────────────────────────────────────────────────────────────────────────────────────────

Found  16  positive integers  N  whose number of divisors is prime (and not equal to 2),  where  N <  1,000

output   when using the input of:     100000

 index │       positive integers  N  whose number of divisors is prime (and not equal to 2),  where  N <  100,000
───────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1   │          4          9         16         25         49         64         81        121        169        289
  11   │        361        529        625        729        841        961      1,024      1,369      1,681      1,849
  21   │      2,209      2,401      2,809      3,481      3,721      4,096      4,489      5,041      5,329      6,241
  31   │      6,889      7,921      9,409     10,201     10,609     11,449     11,881     12,769     14,641     15,625
  41   │     16,129     17,161     18,769     19,321     22,201     22,801     24,649     26,569     27,889     28,561
  51   │     29,929     32,041     32,761     36,481     37,249     38,809     39,601     44,521     49,729     51,529
  61   │     52,441     54,289     57,121     58,081     59,049     63,001     65,536     66,049     69,169     72,361
  71   │     73,441     76,729     78,961     80,089     83,521     85,849     94,249     96,721     97,969
───────┴───────────────────────────────────────────────────────────────────────────────────────────────────────────────

Found  79  positive integers  N  whose number of divisors is prime (and not equal to 2),  where  N <  100,000

Ring

load "stdlib.ring"
row = 0

see "working..." + nl
see "Numbers which count of divisors is prime are:" + nl

for n = 1 to 1000
    num = 0
    for m = 1 to n
        if n%m = 0
           num++
        ok
    next
    if isprime(num) and num != 2
       see "" + n + " "
       row++
       if row%5 = 0
          see nl
       ok
    ok  
next

see nl + "Found " + row + " numbers" + nl
see "done..." + nl

Output:

working...
Numbers which count of divisors is prime are:
4 9 16 25 49 
64 81 121 169 289 
361 529 625 729 841 
961 
Found 16 numbers
done...

RPL

Works with: HP version 49

≪ { } 
   1 999 FOR j
     j DIVIS SIZE
     IF DUP 2 == THEN DROP ELSE
        IF ISPRIME? THEN j + END
     END
   NEXT
≫ 'TASK' STO 

Output:

1; { 4 9 16 25 49 64 81 121 169 289 361 529 625 729 841 961 }

Runs in 2 minutes 16 seconds on a HP-50g.

Rust

fn count_divisors(number : u64 ) -> u64 {
   let mut divisors : u64 = 0 ;
   for n in 1..=number {
      if number % n == 0 {
         divisors += 1 ;
      }
   }
   divisors 
}

fn main() {
   let mut the_numbers : Vec<u64> = Vec::new( ) ;
   let mut count : i32 = 0 ;
   for n in 1..100000 {
      let divis = count_divisors( n as u64 ) ;
      if divis > 2 && primal::is_prime( divis ) {
         the_numbers.push( n ) ;
      }
   }
   for n in the_numbers {
      count += 1 ;
      print!(" {:6}" , n ) ;
      if count % 8 == 0 {
         println!( ) ;
      }
   }
   println!( ) ;
   println!("count is {}" , count) ;
}

Output:

      4      9     16     25     49     64     81    121
    169    289    361    529    625    729    841    961
   1024   1369   1681   1849   2209   2401   2809   3481
   3721   4096   4489   5041   5329   6241   6889   7921
   9409  10201  10609  11449  11881  12769  14641  15625
  16129  17161  18769  19321  22201  22801  24649  26569
  27889  28561  29929  32041  32761  36481  37249  38809
  39601  44521  49729  51529  52441  54289  57121  58081
  59049  63001  65536  66049  69169  72361  73441  76729
  78961  80089  83521  85849  94249  96721  97969
count is 79

Ruby

Testing squares only, according to observation on the discussion page.

require 'prime'

def tau(n) = n.prime_division.inject(1){|res, (d, exp)| res *= exp+1}

res = (1..Integer.sqrt(100_000)).filter_map{|n| sqr = n*n; sqr if tau(sqr).prime? }
res.each_slice(10){|slice| puts "%10d"*slice.size % slice}

Output:

         4         9        16        25        49        64        81       121       169       289
       361       529       625       729       841       961      1024      1369      1681      1849
      2209      2401      2809      3481      3721      4096      4489      5041      5329      6241
      6889      7921      9409     10201     10609     11449     11881     12769     14641     15625
     16129     17161     18769     19321     22201     22801     24649     26569     27889     28561
     29929     32041     32761     36481     37249     38809     39601     44521     49729     51529
     52441     54289     57121     58081     59049     63001     65536     66049     69169     72361
     73441     76729     78961     80089     83521     85849     94249     96721     97969

Sidef

var limit = 100_000
say "Positive integers under #{limit.commify} whose number of divisors is an odd prime:"

1..limit -> grep { !.is_prime && .sigma0.is_prime }.each_slice(10, {|*a|
    say a.map{'%6s' % _}.join(' ')
})

Output:

Positive integers under 100,000 whose number of divisors is an odd prime:
     4      9     16     25     49     64     81    121    169    289
   361    529    625    729    841    961   1024   1369   1681   1849
  2209   2401   2809   3481   3721   4096   4489   5041   5329   6241
  6889   7921   9409  10201  10609  11449  11881  12769  14641  15625
 16129  17161  18769  19321  22201  22801  24649  26569  27889  28561
 29929  32041  32761  36481  37249  38809  39601  44521  49729  51529
 52441  54289  57121  58081  59049  63001  65536  66049  69169  72361
 73441  76729  78961  80089  83521  85849  94249  96721  97969

Wren

Library: Wren-math

Library: Wren-fmt

import "./math" for Int
import "./fmt" for Fmt

var limit = 1e5
var results = []
var i = 2
while (i * i < limit) {
    var n = Int.divisors(i * i).count
    if (n > 2 && Int.isPrime(n)) results.add(i * i)
    i = i + 1
}
Fmt.print("Positive integers under $,7d whose number of divisors is an odd prime:", limit)
Fmt.tprint("$,7d", results, 10)
var under1000 = results.count { |r| r < 1000 }
System.print("\nFound %(results.count) such integers (%(under1000) under 1,000).")

Output:

Positive integers under 100,000 whose number of divisors is an odd prime:
      4       9      16      25      49      64      81     121     169     289
    361     529     625     729     841     961   1,024   1,369   1,681   1,849
  2,209   2,401   2,809   3,481   3,721   4,096   4,489   5,041   5,329   6,241
  6,889   7,921   9,409  10,201  10,609  11,449  11,881  12,769  14,641  15,625
 16,129  17,161  18,769  19,321  22,201  22,801  24,649  26,569  27,889  28,561
 29,929  32,041  32,761  36,481  37,249  38,809  39,601  44,521  49,729  51,529
 52,441  54,289  57,121  58,081  59,049  63,001  65,536  66,049  69,169  72,361
 73,441  76,729  78,961  80,089  83,521  85,849  94,249  96,721  97,969

Found 79 such integers (16 under 1,000).

XPL0

func IsPrime(N);        \Return 'true' if N is a prime number
int  N, I;
[if N <= 1 then return false;
for I:= 2 to sqrt(N) do
    if rem(N/I) = 0 then return false;
return true;
];

func Divisors(N);       \Return number of unique divisors of N
int  N, SN, Count, D;
[SN:= sqrt(N);          \N must be a perfect square to get an odd (prime>2) count
if SN*SN # N then return 0;
Count:= 3;              \SN, 1 and N are unique divisors of N >= 4
for D:= 2 to SN-1 do
    if rem(N/D) = 0 then Count:= Count+2;
return Count;
];

int N, Count;
[Count:= 0;
for N:= 4 to 100_000-1 do
    if IsPrime(Divisors(N)) then
        [Count:= Count+1;
        IntOut(0, N);
        if rem(Count/10) = 0 then CrLf(0) else ChOut(0, 9\tab\);
        ];
]

Output:

4       9       16      25      49      64      81      121     169     289
361     529     625     729     841     961     1024    1369    1681    1849
2209    2401    2809    3481    3721    4096    4489    5041    5329    6241
6889    7921    9409    10201   10609   11449   11881   12769   14641   15625
16129   17161   18769   19321   22201   22801   24649   26569   27889   28561
29929   32041   32761   36481   37249   38809   39601   44521   49729   51529
52441   54289   57121   58081   59049   63001   65536   66049   69169   72361
73441   76729   78961   80089   83521   85849   94249   96721   97969