Sequence: nth number with exactly n divisors

From Rosetta Code
Task
Sequence: nth number with exactly n divisors
You are encouraged to solve this task according to the task description, using any language you may know.

Calculate the sequence where each term an is the nth that has n divisors.

Task

Show here, on this page, at least the first 15 terms of the sequence.

See also
Related tasks




ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Uses Algol 68G's LONG LONG INT (with default precision) to handle the large numbers.
Builds a table of proper divisor counts for testing non-prime n (As with other samples uses the fact that for prime n, the required number is the nth prime to the power n - 1.
Using a table is practical for smallish n, (a table of 500 000 elements is used here, big enough for up to n = 24. For n = 25, a table of over 52 000 000 would be required.
For non-prime n, the divisors are also shown.

BEGIN
    INT max number = 500 000; # maximum number we will count the divisors of #
    # form a table of proper divisor counts                                  #
    [ 1 : max number ]INT pdc; FOR i TO UPB pdc DO pdc[ i ] := 1 OD;
    pdc[ 1 ] := 0;
    FOR i FROM 2 TO UPB pdc DO
        FOR j FROM i + i BY i TO UPB pdc DO pdc[ j ] +:= 1 OD
    OD;
    # find the first few primes                                              #
    [ 1 : 30 ]INT prime;
    INT p count := 0;
    FOR i WHILE p count < UPB prime DO
        IF pdc[ i ] = 1 THEN prime[ p count +:= 1 ] := i FI
    OD;
    # show the nth number with n divisors                                    #
    INT w = -43;  # width to print the numbers (negative means no leading +) #
    print( ( " 1: ", whole( 1, w ), " | 1", newline ) );
    FOR i FROM 2 TO 23 DO
        print( ( whole( i, -2 ), ": " ) );
        IF pdc( i ) = 1 THEN
            print( ( whole( LENG LENG prime[ i ] ^ ( i - 1 ), w ), newline ) )
        ELSE
            INT c := 0;
            FOR j TO UPB pdc WHILE c < i DO
                IF pdc[ j ] = i - 1 THEN
                    c +:= 1;
                    IF c = i THEN
                        print( ( whole( j, w ), " | 1" ) );
                        FOR d FROM 2 TO j OVER 2 DO
                            IF j MOD d = 0 THEN print( ( " ", whole( d, 0 ) ) ) FI
                        OD;
                        print( ( " ", whole( j, 0 ), newline ) )
                    FI
                FI
            OD
        FI
    OD

END
Output:
 1:                                           1 | 1
 2:                                           3
 3:                                          25
 4:                                          14 | 1 2 7 14
 5:                                       14641
 6:                                          44 | 1 2 4 11 22 44
 7:                                    24137569
 8:                                          70 | 1 2 5 7 10 14 35 70
 9:                                        1089 | 1 3 9 11 33 99 121 363 1089
10:                                         405 | 1 3 5 9 15 27 45 81 135 405
11:                             819628286980801
12:                                         160 | 1 2 4 5 8 10 16 20 32 40 80 160
13:                        22563490300366186081
14:                                        2752 | 1 2 4 8 16 32 43 64 86 172 344 688 1376 2752
15:                                        9801 | 1 3 9 11 27 33 81 99 121 297 363 891 1089 3267 9801
16:                                         462 | 1 2 3 6 7 11 14 21 22 33 42 66 77 154 231 462
17:               21559177407076402401757871041
18:                                        1044 | 1 2 3 4 6 9 12 18 29 36 58 87 116 174 261 348 522 1044
19:           740195513856780056217081017732809
20:                                        1520 | 1 2 4 5 8 10 16 19 20 38 40 76 80 95 152 190 304 380 760 1520
21:                                      141376 | 1 2 4 8 16 32 47 64 94 188 376 752 1504 2209 3008 4418 8836 17672 35344 70688 141376
22:                                       84992 | 1 2 4 8 16 32 64 83 128 166 256 332 512 664 1024 1328 2656 5312 10624 21248 42496 84992
23: 1658509762573818415340429240403156732495289

C

Translation of: C++
#include <math.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>

#define LIMIT 15
int smallPrimes[LIMIT];

static void sieve() {
    int i = 2, j;
    int p = 5;

    smallPrimes[0] = 2;
    smallPrimes[1] = 3;

    while (i < LIMIT) {
        for (j = 0; j < i; j++) {
            if (smallPrimes[j] * smallPrimes[j] <= p) {
                if (p % smallPrimes[j] == 0) {
                    p += 2;
                    break;
                }
            } else {
                smallPrimes[i++] = p;
                p += 2;
                break;
            }
        }
    }
}

static bool is_prime(uint64_t n) {
    uint64_t i;

    for (i = 0; i < LIMIT; i++) {
        if (n % smallPrimes[i] == 0) {
            return n == smallPrimes[i];
        }
    }

    i = smallPrimes[LIMIT - 1] + 2;
    for (; i * i <= n; i += 2) {
        if (n % i == 0) {
            return false;
        }
    }

    return true;
}

static uint64_t divisor_count(uint64_t n) {
    uint64_t count = 1;
    uint64_t d;

    while (n % 2 == 0) {
        n /= 2;
        count++;
    }

    for (d = 3; d * d <= n; d += 2) {
        uint64_t q = n / d;
        uint64_t r = n % d;
        uint64_t dc = 0;
        while (r == 0) {
            dc += count;
            n = q;
            q = n / d;
            r = n % d;
        }
        count += dc;
    }

    if (n != 1) {
        return count *= 2;
    }
    return count;
}

static uint64_t OEISA073916(size_t n) {
    uint64_t count = 0;
    uint64_t result = 0;
    size_t i;

    if (is_prime(n)) {
        return (uint64_t)pow(smallPrimes[n - 1], n - 1);
    }

    for (i = 1; count < n; i++) {
        if (n % 2 == 1) {
            //  The solution for an odd (non-prime) term is always a square number
            uint64_t root = (uint64_t)sqrt(i);
            if (root * root != i) {
                continue;
            }
        }
        if (divisor_count(i) == n) {
            count++;
            result = i;
        }
    }

    return result;
}

int main() {
    size_t n;

    sieve();

    for (n = 1; n <= LIMIT; n++) {
        if (n == 13) {
            printf("A073916(%lu) = One more bit needed to represent result.\n", n);
        } else {
            printf("A073916(%lu) = %llu\n", n, OEISA073916(n));
        }
    }

    return 0;
}
Output:
A073916(1) = 1
A073916(2) = 3
A073916(3) = 25
A073916(4) = 14
A073916(5) = 14641
A073916(6) = 44
A073916(7) = 24137569
A073916(8) = 70
A073916(9) = 1089
A073916(10) = 405
A073916(11) = 819628286980801
A073916(12) = 160
A073916(13) = One more bit needed to represent result.
A073916(14) = 2752
A073916(15) = 9801

C++

Translation of: Java
#include <iostream>
#include <vector>

std::vector<int> smallPrimes;

bool is_prime(size_t test) {
    if (test < 2) {
        return false;
    }
    if (test % 2 == 0) {
        return test == 2;
    }
    for (size_t d = 3; d * d <= test; d += 2) {
        if (test % d == 0) {
            return false;
        }
    }
    return true;
}

void init_small_primes(size_t numPrimes) {
    smallPrimes.push_back(2);

    int count = 0;
    for (size_t n = 3; count < numPrimes; n += 2) {
        if (is_prime(n)) {
            smallPrimes.push_back(n);
            count++;
        }
    }
}

size_t divisor_count(size_t n) {
    size_t count = 1;
    while (n % 2 == 0) {
        n /= 2;
        count++;
    }
    for (size_t d = 3; d * d <= n; d += 2) {
        size_t q = n / d;
        size_t r = n % d;
        size_t dc = 0;
        while (r == 0) {
            dc += count;
            n = q;
            q = n / d;
            r = n % d;
        }
        count += dc;
    }
    if (n != 1) {
        count *= 2;
    }
    return count;
}

uint64_t OEISA073916(size_t n) {
    if (is_prime(n)) {
        return (uint64_t) pow(smallPrimes[n - 1], n - 1);
    }

    size_t count = 0;
    uint64_t result = 0;
    for (size_t i = 1; count < n; i++) {
        if (n % 2 == 1) {
            //  The solution for an odd (non-prime) term is always a square number
            size_t root = (size_t) sqrt(i);
            if (root * root != i) {
                continue;
            }
        }
        if (divisor_count(i) == n) {
            count++;
            result = i;
        }
    }
    return result;
}

int main() {
    const int MAX = 15;
    init_small_primes(MAX);
    for (size_t n = 1; n <= MAX; n++) {
        if (n == 13) {
            std::cout << "A073916(" << n << ") = One more bit needed to represent result.\n";
        } else {
            std::cout << "A073916(" << n << ") = " << OEISA073916(n) << '\n';
        }
    }

    return 0;
}
Output:
A073916(1) = 1
A073916(2) = 3
A073916(3) = 25
A073916(4) = 14
A073916(5) = 14641
A073916(6) = 44
A073916(7) = 24137569
A073916(8) = 70
A073916(9) = 1089
A073916(10) = 405
A073916(11) = 819628286980801
A073916(12) = 160
A073916(13) = One more bit needed to represent result.
A073916(14) = 2752
A073916(15) = 9801

D

Translation of: Java
import std.bigint;
import std.math;
import std.stdio;

bool isPrime(long test) {
    if (test == 2) {
        return true;
    }
    if (test % 2 == 0) {
        return false;
    }
    for (long d = 3 ; d * d <= test; d += 2) {
        if (test % d == 0) {
            return false;
        }
    }
    return true;
}

int[] calcSmallPrimes(int numPrimes) {
    int[] smallPrimes;
    smallPrimes ~= 2;

    int count = 0;
    int n = 3;
    while (count < numPrimes) {
        if (isPrime(n)) {
            smallPrimes ~= n;
            count++;
        }
        n += 2;
    }

    return smallPrimes;
}

immutable MAX = 45;
immutable smallPrimes = calcSmallPrimes(MAX);

int getDivisorCount(long n) {
    int count = 1;
    while (n % 2 == 0) {
        n /= 2;
        count += 1;
    }
    for (long d = 3; d * d <= n; d += 2) {
        long q = n / d;
        long r = n % d;
        int dc = 0;
        while (r == 0) {
            dc += count;
            n = q;
            q = n / d;
            r = n % d;
        }
        count += dc;
    }
    if (n != 1) {
        count *= 2;
    }
    return count;
}

BigInt OEISA073916(int n) {
    if (isPrime(n) ) {
        return BigInt(smallPrimes[n-1]) ^^ (n - 1);
    }
    int count = 0;
    int result = 0;
    for (int i = 1; count < n; i++) {
        if (n % 2 == 1) {
            //  The solution for an odd (non-prime) term is always a square number
            int root = cast(int) sqrt(cast(real) i);
            if (root * root != i) {
                continue;
            }
        }
        if (getDivisorCount(i) == n) {
            count++;
            result = i;
        }
    }
    return BigInt(result);
}

void main() {
    foreach (n; 1 .. MAX + 1) {
        writeln("A073916(", n, ") = ", OEISA073916(n));
    }
}
Output:
A073916(1) = 1
A073916(2) = 3
A073916(3) = 25
A073916(4) = 14
A073916(5) = 14641
A073916(6) = 44
A073916(7) = 24137569
A073916(8) = 70
A073916(9) = 1089
A073916(10) = 405
A073916(11) = 819628286980801
A073916(12) = 160
A073916(13) = 22563490300366186081
A073916(14) = 2752
A073916(15) = 9801
A073916(16) = 462
A073916(17) = 21559177407076402401757871041
A073916(18) = 1044
A073916(19) = 740195513856780056217081017732809
A073916(20) = 1520
A073916(21) = 141376
A073916(22) = 84992
A073916(23) = 1658509762573818415340429240403156732495289
A073916(24) = 1170
A073916(25) = 52200625
A073916(26) = 421888
A073916(27) = 52900
A073916(28) = 9152
A073916(29) = 1116713952456127112240969687448211536647543601817400964721
A073916(30) = 6768
A073916(31) = 1300503809464370725741704158412711229899345159119325157292552449
A073916(32) = 3990
A073916(33) = 12166144
A073916(34) = 9764864
A073916(35) = 446265625
A073916(36) = 5472
A073916(37) = 11282036144040442334289838466416927162302790252609308623697164994458730076798801
A073916(38) = 43778048
A073916(39) = 90935296
A073916(40) = 10416
A073916(41) = 1300532588674810624476094551095787816112173600565095470117230812218524514342511947837104801
A073916(42) = 46400
A073916(43) = 635918448514386699807643535977466343285944704172890141356181792680152445568879925105775366910081
A073916(44) = 240640
A073916(45) = 327184

Factor

This makes use of most of the optimizations discussed in the Go example.

USING: combinators formatting fry kernel lists lists.lazy
lists.lazy.examples literals math math.functions math.primes
math.primes.factors math.ranges sequences ;
IN: rosetta-code.nth-n-div

CONSTANT: primes $[ 100 nprimes ]

: prime ( m -- n ) 1 - [ primes nth ] [ ^ ] bi ;

: (non-prime) ( m quot -- n )
    '[
        [ 1 - ] [ drop @ ] [ ] tri '[ divisors length _ = ]
        lfilter swap [ cdr ] times car
    ] call ; inline

: non-prime ( m quot -- n )
    {
        { [ over 2 = ] [ 2drop 3 ] }
        { [ over 10 = ] [ 2drop 405 ] }
        [ (non-prime) ]
    } cond ; inline

: fn ( m -- n )
    {
        { [ dup even? ] [ [ evens ] non-prime ] }
        { [ dup prime? ] [ prime ] }
        [ [ squares ] non-prime ]
    } cond ;

: main ( -- ) 45 [1,b] [ dup fn "%2d : %d\n" printf ] each ;

MAIN: main
Output:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14641
 6 : 44
 7 : 24137569
 8 : 70
 9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170
25 : 52200625
26 : 421888
27 : 52900
28 : 9152
29 : 1116713952456127112240969687448211536647543601817400964721
30 : 6768
31 : 1300503809464370725741704158412711229899345159119325157292552449
32 : 3990
33 : 12166144
34 : 9764864
35 : 446265625
36 : 5472
37 : 11282036144040442334289838466416927162302790252609308623697164994458730076798801
38 : 43778048
39 : 90935296
40 : 10416
41 : 1300532588674810624476094551095787816112173600565095470117230812218524514342511947837104801
42 : 46400
43 : 635918448514386699807643535977466343285944704172890141356181792680152445568879925105775366910081
44 : 240640
45 : 327184


FreeBASIC

But it takes an "eternity" to resolve.

Dim As Integer n, num, pnum
Dim As Ulongint m, p
Dim As Ulongint limit = 18446744073709551615

Print "The first 15 terms of OEIS:A073916 are:"
For n = 1 To 15
    num = 0
    For m = 1 To limit
        pnum = 0
        For p = 1 To limit
            If (m Mod p = 0) Then pnum += 1
        Next p
        If pnum = n Then num += 1
        If num = n Then
            Print Using "## : &"; n; m
            Exit  For
        End If 
    Next m
Next n
Sleep


Go

This makes use of the relationship: a[p] = prime[p]^(p-1) if p is prime, mentioned in the blurb for A073916 (and also on the talk page) to calculate the larger terms, some of which require big.Int in Go. It also makes use of another hint on the talk page that all odd terms are square numbers.

The remaining terms (up to the 33rd) are not particularly large and so are calculated by brute force.

package main

import (
    "fmt"
    "math"
    "math/big"
)

var bi = new(big.Int)

func isPrime(n int) bool {
    bi.SetUint64(uint64(n))
    return bi.ProbablyPrime(0)
}

func generateSmallPrimes(n int) []int {
    primes := make([]int, n)
    primes[0] = 2
    for i, count := 3, 1; count < n; i += 2 {
        if isPrime(i) {
            primes[count] = i
            count++
        }
    }
    return primes
}

func countDivisors(n int) int {
    count := 1
    for n%2 == 0 {
        n >>= 1
        count++
    }
    for d := 3; d*d <= n; d += 2 {
        q, r := n/d, n%d
        if r == 0 {
            dc := 0
            for r == 0 {
                dc += count
                n = q
                q, r = n/d, n%d
            }
            count += dc
        }
    }
    if n != 1 {
        count *= 2
    }
    return count
}

func main() {
    const max = 33
    primes := generateSmallPrimes(max)
    z := new(big.Int)
    p := new(big.Int)
    fmt.Println("The first", max, "terms in the sequence are:")
    for i := 1; i <= max; i++ {
        if isPrime(i) {
            z.SetUint64(uint64(primes[i-1]))
            p.SetUint64(uint64(i - 1))
            z.Exp(z, p, nil)
            fmt.Printf("%2d : %d\n", i, z)
        } else {
            count := 0
            for j := 1; ; j++ {
                if i%2 == 1 {
                    sq := int(math.Sqrt(float64(j)))
                    if sq*sq != j {
                        continue
                    }
                }
                if countDivisors(j) == i {
                    count++
                    if count == i {
                        fmt.Printf("%2d : %d\n", i, j)
                        break
                    }
                }
            }
        }
    }
}
Output:
The first 33 terms in the sequence are:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14641
 6 : 44
 7 : 24137569
 8 : 70
 9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170
25 : 52200625
26 : 421888
27 : 52900
28 : 9152
29 : 1116713952456127112240969687448211536647543601817400964721
30 : 6768
31 : 1300503809464370725741704158412711229899345159119325157292552449
32 : 3990
33 : 12166144

The following much faster version (runs in less than 90 seconds on my 1.6GHz Celeron) uses three further optimizations:

1. Apart from the 2nd and 10th terms, all the even terms are themselves even.

2. A sieve is used to generate all prime divisors needed. This doesn't take up much time or memory but speeds up the counting of all divisors considerably.

3. While searching for the nth number with exactly n divisors, where feasible a record is kept of any numbers found to have exactly k divisors (k > n) so that the search for these numbers can start from a higher base.

package main

import (
    "fmt"
    "math"
    "math/big"
)

type record struct{ num, count int }

var (
    bi     = new(big.Int)
    primes = []int{2}
)

func isPrime(n int) bool {
    bi.SetUint64(uint64(n))
    return bi.ProbablyPrime(0)
}

func sieve(limit int) {
    c := make([]bool, limit+1) // composite = true
    // no need to process even numbers
    p := 3
    for {
        p2 := p * p
        if p2 > limit {
            break
        }
        for i := p2; i <= limit; i += 2 * p {
            c[i] = true
        }
        for {
            p += 2
            if !c[p] {
                break
            }
        }
    }
    for i := 3; i <= limit; i += 2 {
        if !c[i] {
            primes = append(primes, i)
        }
    }
}

func countDivisors(n int) int {
    count := 1
    for i, p := 0, primes[0]; p*p <= n; i, p = i+1, primes[i+1] {
        if n%p != 0 {
            continue
        }
        n /= p
        count2 := 1
        for n%p == 0 {
            n /= p
            count2++
        }
        count *= (count2 + 1)
        if n == 1 {
            return count
        }
    }
    if n != 1 {
        count *= 2
    }
    return count
}

func isOdd(x int) bool {
    return x%2 == 1
}

func main() {
    sieve(22000)
    const max = 45
    records := [max + 1]record{}
    z := new(big.Int)
    p := new(big.Int)
    fmt.Println("The first", max, "terms in the sequence are:")
    for i := 1; i <= max; i++ {
        if isPrime(i) {
            z.SetUint64(uint64(primes[i-1]))
            p.SetUint64(uint64(i - 1))
            z.Exp(z, p, nil)
            fmt.Printf("%2d : %d\n", i, z)
        } else {
            count := records[i].count
            if count == i {
                fmt.Printf("%2d : %d\n", i, records[i].num)
                continue
            }
            odd := isOdd(i)
            k := records[i].num
            l := 1
            if !odd && i != 2 && i != 10 {
                l = 2
            }
            for j := k + l; ; j += l {
                if odd {
                    sq := int(math.Sqrt(float64(j)))
                    if sq*sq != j {
                        continue
                    }
                }
                cd := countDivisors(j)
                if cd == i {
                    count++
                    if count == i {
                        fmt.Printf("%2d : %d\n", i, j)
                        break
                    }
                } else if cd > i && cd <= max && records[cd].count < cd &&
                    j > records[cd].num && (l == 1 || (l == 2 && !isOdd(cd))) {
                    records[cd].num = j
                    records[cd].count++
                }
            }
        }
    }
}
Output:
The first 45 terms in the sequence are:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14641
 6 : 44
 7 : 24137569
 8 : 70
 9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170
25 : 52200625
26 : 421888
27 : 52900
28 : 9152
29 : 1116713952456127112240969687448211536647543601817400964721
30 : 6768
31 : 1300503809464370725741704158412711229899345159119325157292552449
32 : 3990
33 : 12166144
34 : 9764864
35 : 446265625
36 : 5472
37 : 11282036144040442334289838466416927162302790252609308623697164994458730076798801
38 : 43778048
39 : 90935296
40 : 10416
41 : 1300532588674810624476094551095787816112173600565095470117230812218524514342511947837104801
42 : 46400
43 : 635918448514386699807643535977466343285944704172890141356181792680152445568879925105775366910081
44 : 240640
45 : 327184

Haskell

import           Control.Monad                         (guard)
import           Math.NumberTheory.ArithmeticFunctions (divisorCount)
import           Math.NumberTheory.Primes              (Prime, unPrime)
import           Math.NumberTheory.Primes.Testing      (isPrime)

calc :: Integer -> [(Integer, Integer)]
calc n = do
  x <- [1..]
  guard (even n || odd n && f x == x)
  [(x, divisorCount x)]
 where f n = floor (sqrt $ realToFrac n) ^ 2

havingNthDivisors :: Integer -> [(Integer, Integer)]
havingNthDivisors n = filter ((==n) . snd) $ calc n

nths :: [(Integer, Integer)]
nths = do
  n <- [1..35] :: [Integer]
  if isPrime n then
    pure (n, nthPrime (fromIntegral n) ^ pred n)
  else
    pure (n, f n)
 where
  f n = fst (havingNthDivisors n !! pred (fromIntegral n))
  nthPrime n = unPrime (toEnum n :: Prime Integer)

main :: IO ()
main = mapM_ print nths
Output:
(1,1)
(2,3)
(3,25)
(4,14)
(5,14641)
(6,44)
(7,24137569)
(8,70)
(9,1089)
(10,405)
(11,819628286980801)
(12,160)
(13,22563490300366186081)
(14,2752)
(15,9801)
(16,462)
(17,21559177407076402401757871041)
(18,1044)
(19,740195513856780056217081017732809)
(20,1520)
(21,141376)
(22,84992)
(23,1658509762573818415340429240403156732495289)
(24,1170)
(25,52200625)
(26,421888)
(27,52900)
(28,9152)
(29,1116713952456127112240969687448211536647543601817400964721)
(30,6768)
(31,1300503809464370725741704158412711229899345159119325157292552449)
(32,3990)
(33,12166144)
(34,9764864)
(35,446265625)

J

Implementation:

A073916=: {{
  if.1 p: y do. (p:^x:)y-1 return.
  elseif.1|y do.
    f= *:
  else.
    f=. ]
  end. r=.i.0
  off=. 1
  while. y>#r do.
    r=. r,f off+I.y=*/|:1+_ q:f off+i.y
    off=. off+y
  end.
  (y-1){r
}}"0

Task example:

   (,. A073916) 1+i.20
 1                                 1
 2                                 3
 3                                25
 4                                14
 5                             14641
 6                                44
 7                          24137569
 8                                70
 9                              1089
10                               405
11                   819628286980801
12                               160
13              22563490300366186081
14                              2752
15                              9801
16                               462
17     21559177407076402401757871041
18                              1044
19 740195513856780056217081017732809
20                              1520

Java

Replace translation with Java native implementation.

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class SequenceNthNumberWithExactlyNDivisors {

    public static void main(String[] args) {
        int max = 45;
        smallPrimes(max);
        for ( int n = 1; n <= max ; n++ ) {
            System.out.printf("A073916(%d) = %s%n", n, OEISA073916(n));
        }
    }
    
    private static List<Integer> smallPrimes = new ArrayList<>();
    
    private static void smallPrimes(int numPrimes) {
        smallPrimes.add(2);
        for ( int n = 3, count = 0 ; count < numPrimes ; n += 2 ) {
            if ( isPrime(n) ) {
                smallPrimes.add(n);
                count++;
            }
        }
    }
    
    private static final boolean isPrime(long test) {
        if ( test == 2 ) {
            return true;
        }
        if ( test % 2 == 0 ) {
            return false;
        }
        for ( long d = 3 ; d*d <= test ; d += 2 ) {
            if ( test % d == 0 ) {
                return false;
            }
        }
        return true;
    }

    private static int getDivisorCount(long n) {
        int count = 1;
        while ( n % 2 == 0 ) {
            n /= 2;
            count += 1;
        }
        for ( long d = 3 ; d*d <= n ; d += 2 ) {
            long q = n / d;
            long r = n % d;
            int dc = 0;
            while ( r == 0 ) {
                dc += count;
                n = q;
                q = n / d;
                r = n % d;
            }
            count += dc;
        }
        if ( n != 1 ) {
            count *= 2;
        }
        return count;
    }
    
    private static BigInteger OEISA073916(int n) {
        if ( isPrime(n) ) {
            return BigInteger.valueOf(smallPrimes.get(n-1)).pow(n - 1);
        }
        int count = 0;
        int result = 0;
        for ( int i = 1 ; count < n ; i++ ) {
            if ( n % 2 == 1 ) {
                //  The solution for an odd (non-prime) term is always a square number
                int sqrt = (int) Math.sqrt(i);
                if ( sqrt*sqrt != i ) {
                    continue;
                }
            }
            if ( getDivisorCount(i) == n ) {
                count++;
                result = i;
            }
        }
        return BigInteger.valueOf(result);
    }

}
Output:
A073916(1) = 1
A073916(2) = 3
A073916(3) = 25
A073916(4) = 14
A073916(5) = 14641
A073916(6) = 44
A073916(7) = 24137569
A073916(8) = 70
A073916(9) = 1089
A073916(10) = 405
A073916(11) = 819628286980801
A073916(12) = 160
A073916(13) = 22563490300366186081
A073916(14) = 2752
A073916(15) = 9801
A073916(16) = 462
A073916(17) = 21559177407076402401757871041
A073916(18) = 1044
A073916(19) = 740195513856780056217081017732809
A073916(20) = 1520
A073916(21) = 141376
A073916(22) = 84992
A073916(23) = 1658509762573818415340429240403156732495289
A073916(24) = 1170
A073916(25) = 52200625
A073916(26) = 421888
A073916(27) = 52900
A073916(28) = 9152
A073916(29) = 1116713952456127112240969687448211536647543601817400964721
A073916(30) = 6768
A073916(31) = 1300503809464370725741704158412711229899345159119325157292552449
A073916(32) = 3990
A073916(33) = 12166144
A073916(34) = 9764864
A073916(35) = 446265625
A073916(36) = 5472
A073916(37) = 11282036144040442334289838466416927162302790252609308623697164994458730076798801
A073916(38) = 43778048
A073916(39) = 90935296
A073916(40) = 10416
A073916(41) = 1300532588674810624476094551095787816112173600565095470117230812218524514342511947837104801
A073916(42) = 46400
A073916(43) = 635918448514386699807643535977466343285944704172890141356181792680152445568879925105775366910081
A073916(44) = 240640
A073916(45) = 327184

jq

Adapted from Wren
Works with gojq, the Go implementation of jq

See Erdős-primes#jq for a suitable implementation of `is_prime`.

The precision of the integer arithmetic of the C implementation of jq is only precise enough for computing the n-th value up to and including [16,462]. Accordingly gojq was used to produce the output shown below.

Preliminaries

def count(stream): reduce stream as $i (0; .+1);

# To maintain precision:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);

def primes: 2, (range(3; infinite; 2) | select(is_prime));

# divisors as an unsorted stream
def divisors:
  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 $i
         elif ($n % $i) == 0
         then $i, ($n/$i)
         else empty
         end
    end
  end;

The Task

# Emit [n, nth_with_n_divisors] for n in range(1; .+1)
def nth_with_n_divisors:
  | [limit( .; primes)] as $primes
  | range( 1; 1 + .) as $i
  | if $i | is_prime
    then [$i, ($primes[$i-1]|power($i-1))]
    else {count: 0, j: 1}
    | until(.count == $i ;
        .cont = false
        | if ($i % 2) == 1 then (.j|sqrt|floor) as $sq
          | if ($sq * $sq) != .j then .j += 1 | .cont = true else . end
	  else .
	  end
          | if .cont == false
            then if (.j | count(divisors)) == $i
                 then .count += 1
                 else . 
                 end
            | if .count != $i then .j += 1 else . end
            else .
            end )

    | [ $i, .j]
    end;

"The first 33 terms in the sequence are:",
(33 | nth_with_n_divisors)
Output:
The first 33 terms in the sequence are:
[1,1]
[2,3]
[3,25]
[4,14]
[5,14641]
[6,44]
[7,24137569]
[8,70]
[9,1089]
[10,405]
[11,819628286980801]
[12,160]
[13,22563490300366186081]
[14,2752]
[15,9801]
[16,462]
[17,21559177407076402401757871041]
[18,1044]
[19,740195513856780056217081017732809]
[20,1520]
[21,141376]
[22,84992]
[23,1658509762573818415340429240403156732495289]
[24,1170]
[25,52200625]
[26,421888]
[27,52900]
[28,9152]
[29,1116713952456127112240969687448211536647543601817400964721]
[30,6768]
[31,1300503809464370725741704158412711229899345159119325157292552449]
[32,3990]
[33,12166144]

Julia

using Primes

function countdivisors(n)
    f = [one(n)]
    for (p, e) in factor(n)
        f = reduce(vcat, [f * p ^ j for j in 1:e], init = f)
    end
    length(f)
end

function nthwithndivisors(N)
    parray = findall(primesmask(100 * N))
    for i = 1:N
        if isprime(i)
            println("$i : ", BigInt(parray[i])^(i-1))
        else
            k = 0
            for j in 1:100000000000
                if (iseven(i) || Int(floor(sqrt(j)))^2 == j) &&
                    i == countdivisors(j) && (k += 1) == i
                    println("$i : $j")
                    break
                end
            end
        end
    end
end

nthwithndivisors(35)
Output:
1 : 1
2 : 3
3 : 25
4 : 14
5 : 14641
6 : 44
7 : 24137569
8 : 70
9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170
25 : 52200625
26 : 421888
27 : 52900
28 : 9152
29 : 1116713952456127112240969687448211536647543601817400964721
30 : 6768
31 : 1300503809464370725741704158412711229899345159119325157292552449
32 : 3990
33 : 12166144
34 : 9764864
35 : 446265625

Kotlin

Translation of: Go
// Version 1.3.21

import java.math.BigInteger
import kotlin.math.sqrt

const val MAX = 33

fun isPrime(n: Int) = BigInteger.valueOf(n.toLong()).isProbablePrime(10)

fun generateSmallPrimes(n: Int): List<Int> {
    val primes = mutableListOf<Int>()
    primes.add(2)
    var i = 3
    while (primes.size < n) {
        if (isPrime(i)) {
            primes.add(i)
        }
        i += 2
    }
    return primes
}

fun countDivisors(n: Int): Int {
    var nn = n
    var count = 1
    while (nn % 2 == 0) {
        nn = nn shr 1
        count++
    }
    var d = 3
    while (d * d <= nn) {
        var q = nn / d
        var r = nn % d
        if (r == 0) {
            var dc = 0
            while (r == 0) {
                dc += count
                nn = q
                q = nn / d
                r = nn % d
            }
            count += dc
        }
        d += 2
    }
    if (nn != 1) count *= 2
    return count
}

fun main() {
    var primes = generateSmallPrimes(MAX)
    println("The first $MAX terms in the sequence are:")
    for (i in 1..MAX) {
        if (isPrime(i)) {
            var z = BigInteger.valueOf(primes[i - 1].toLong())
            z = z.pow(i - 1)
            System.out.printf("%2d : %d\n", i, z)
        } else {
            var count = 0
            var j = 1
            while (true) {
                if (i % 2 == 1) {
                    val sq = sqrt(j.toDouble()).toInt()
                    if (sq * sq != j) {
                        j++
                        continue
                    }
                }
                if (countDivisors(j) == i) {
                    if (++count == i) {
                        System.out.printf("%2d : %d\n", i, j)
                        break
                    }
                }
                j++
            }
        }
    }
}
Output:
The first 33 terms in the sequence are:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14641
 6 : 44
 7 : 24137569
 8 : 70
 9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170
25 : 52200625
26 : 421888
27 : 52900
28 : 9152
29 : 1116713952456127112240969687448211536647543601817400964721
30 : 6768
31 : 1300503809464370725741704158412711229899345159119325157292552449
32 : 3990
33 : 12166144

Mathematica / Wolfram Language

d = Table[
  Length[Divisors[n]], {n, 
   200000}]; t = {}; n = 0; ok = True; While[ok, n++; 
 If[PrimeQ[n], AppendTo[t, Prime[n]^(n - 1)], 
  c = Flatten[Position[d, n, 1, n]]; 
  If[Length[c] >= n, AppendTo[t, c[[n]]], ok = False]]]; 
t
Output:
{1,3,25,14,14641,44,24137569,70,1089,405,819628286980801,160,22563490300366186081,2752,9801,462,21559177407076402401757871041,1044,740195513856780056217081017732809,1520,141376,84992,1658509762573818415340429240403156732495289,1170}

Nim

Translation of: Go
Library: bignum

This is a translation of the fast Go version. It runs in about 23s on our laptop.

import math, strformat
import bignum

type Record = tuple[num, count: Natural]

template isOdd(n: Natural): bool =
  (n and 1) != 0

func isPrime(n: int): bool =
  let bi = newInt(n)
  result = bi.probablyPrime(25) != 0

proc findPrimes(limit: Natural): seq[int] {.compileTime.} =
  result = @[2]
  var isComposite = newSeq[bool](limit + 1)
  var p = 3
  while true:
    let p2 = p * p
    if p2 > limit: break
    for i in countup(p2, limit, 2 * p):
      isComposite[i] = true
    while true:
      inc p, 2
      if not isComposite[p]: break
  for n in countup(3, limit, 2):
    if not isComposite[n]:
      result.add n

const Primes = findPrimes(22_000)

proc countDivisors(n: Natural): int =
  result = 1
  var n = n
  for i, p in Primes:
    if p * p > n: break
    if n mod p != 0: continue
    n = n div p
    var count = 1
    while n mod p == 0:
      n = n div p
      inc count
    result *= count + 1
    if n == 1: return
  if n != 1: result *= 2

const Max = 45
var records: array[0..Max, Record]
echo &"The first {Max} terms in the sequence are:"

for n in 1..Max:

  if n.isPrime:
    var z = newInt(Primes[n - 1])
    z = pow(z, culong(n - 1))
    echo &"{n:2}: {z}"

  else:
    var count = records[n].count
    if count == n:
      echo &"{n:2}: {records[n].num}"
      continue
    let odd = n.isOdd
    let d = if odd or n == 2 or n == 10: 1 else: 2
    var k = records[n].num
    while true:
      inc k, d
      if odd:
        let sq = sqrt(k.toFloat).int
        if sq * sq != k: continue
      let cd = k.countDivisors()
      if cd == n:
        inc count
        if count == n:
          echo &"{n:2}: {k}"
          break
      elif cd in (n + 1)..Max and records[cd].count < cd and
           k > records[cd].num and (d == 1 or d == 2 and not cd.isOdd):
        records[cd].num = k
        inc records[cd].count
Output:
The first 45 terms in the sequence are:
 1: 1
 2: 3
 3: 25
 4: 14
 5: 14641
 6: 44
 7: 24137569
 8: 70
 9: 1089
10: 405
11: 819628286980801
12: 160
13: 22563490300366186081
14: 2752
15: 9801
16: 462
17: 21559177407076402401757871041
18: 1044
19: 740195513856780056217081017732809
20: 1520
21: 141376
22: 84992
23: 1658509762573818415340429240403156732495289
24: 1170
25: 52200625
26: 421888
27: 52900
28: 9152
29: 1116713952456127112240969687448211536647543601817400964721
30: 6768
31: 1300503809464370725741704158412711229899345159119325157292552449
32: 3990
33: 12166144
34: 9764864
35: 446265625
36: 5472
37: 11282036144040442334289838466416927162302790252609308623697164994458730076798801
38: 43778048
39: 90935296
40: 10416
41: 1300532588674810624476094551095787816112173600565095470117230812218524514342511947837104801
42: 46400
43: 635918448514386699807643535977466343285944704172890141356181792680152445568879925105775366910081
44: 240640
45: 327184

Perl

Library: ntheory
Translation of: Raku
use strict;
use warnings;
use bigint;
use ntheory <nth_prime is_prime divisors>;

my $limit = 20;

print "First $limit terms of OEIS:A073916\n";

for my $n (1..$limit) {
    if ($n > 4 and is_prime($n)) {
        print nth_prime($n)**($n-1) . ' ';
    } else {
        my $i = my $x = 0;
        while (1) {
            my $nn = $n%2 ? ++$x**2 : ++$x;
            next unless $n == divisors($nn) and ++$i == $n;
            print "$nn " and last;
      }
    }
}
Output:
First 20 terms of OEIS:A073916
1 3 25 14 14641 44 24137569 70 1089 405 819628286980801 160 22563490300366186081 2752 9801 462 21559177407076402401757871041 1044 740195513856780056217081017732809 1520

Phix

Library: Phix/mpfr

simple

Certainly not the fastest way to do it, hence the relatively small limit of 24, which takes less than 0.4s,
whereas a limit of 25 would need to invoke factors() 52 million times which would no doubt take a fair while.

with javascript_semantics 
constant LIMIT = 24
include mpfr.e
mpz z = mpz_init()
 
sequence fn = repeat(0,LIMIT) fn[1] = 1
integer k = 1
printf(1,"The first %d terms in the sequence are:\n",LIMIT)
for i=1 to LIMIT do
    if is_prime(i) then
        mpz_ui_pow_ui(z,get_prime(i),i-1)
        printf(1,"%2d : %s\n",{i,mpz_get_str(z)})
    else
        while fn[i]<i do
            k += 1
            integer l = length(factors(k,1))
            if l<=LIMIT and fn[l]<l then
                fn[l] = iff(fn[l]+1<l?fn[l]+1:k)
            end if
        end while
        printf(1,"%2d : %d\n",{i,fn[i]})
    end if
end for
Output:
The first 24 terms in the sequence are:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14641
 6 : 44
 7 : 24137569
 8 : 70
 9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170

cheating slightly

No real patterns that I could see here, but you can still identify and single out the troublemakers (of which there are about 30).

with javascript_semantics 
include mpfr.e
atom t0 = time()
constant LIMIT = 100
include mpfr.e
include primes.e
mpz z = mpz_init(),
    p = mpz_init() 
string mz
sequence fn = repeat(0,LIMIT), dx;  fn[1] = 1
integer k = 1, idx, p1, p2
printf(1,"The first %d terms in the sequence are:\n",LIMIT)
for i=1 to LIMIT do
    if is_prime(i) or i=1 then
        mpz_ui_pow_ui(z,get_prime(i),i-1)
        mz = mpz_get_str(z)
    else
        sequence f = prime_factors(i,1)
        if length(f)=2 and f[1]=2 and f[2]>7 then
            mz = sprintf("%d",power(2,f[2]-1)*get_prime(i+1))
        elsif length(f)=2 and f[1]>2 then
            if f[1]=f[2] then
                mz = sprintf("%d",power(f[1]*get_prime(f[1]+2),f[1]-1))
            else -- deal with some tardy ones...
                dx = {15,21,33,35,39,51,55,57,65,69,77,85,87,91,93,95}; idx = find(i,dx)
                p1 = { 3, 2, 2, 5, 2, 2, 2, 2, 2, 2, 7, 2, 2, 7, 2, 2}[idx]
                p2 = { 5,15,29, 6,35,49,34,56,45,69, 7,65,88, 7,94,77}[idx]
                mpz_ui_pow_ui(z,p1,f[2]-1)
                mpz_ui_pow_ui(p,get_prime(p2),f[1]-1)
                mpz_mul(z,z,p)
                mz = mpz_get_str(z)
            end if
        elsif (length(f)=3 and i>50) or (length(f)=4 and (f[1]=3 or f[4]>7)) then 
            if i=99 then    -- (oops, messed that one up!)
                mz = sprintf("%d",4*power(3,10)*31*31)
            elsif i=63 then -- (and another!)
                mz = sprintf("%d",power(2,8)*power(5,6))
            else
                dx = {52,66,68,70,75,76,78,92,98,81,88}; idx = find(i,dx)
                p1 = { 7, 3, 1, 5, 3, 5, 5,13, 3,35,35}[idx]
                p2 = { 1, 2, 1, 4, 4, 1, 2, 1, 1, 2, 1}[idx]
                mpz_ui_pow_ui(z,2,f[$]-1)
                mpz_ui_pow_ui(p,p1,p2)
                mpz_mul(z,z,p)
                p1 = {13,37, 4, 9,34,22,19,12, 4,11,13}[idx]
                p2 = { 1, 1, 3, 1, 2, 1, 1, 1, 6, 2, 1}[idx]
                mpz_ui_pow_ui(p,get_prime(p1),p2)
                mpz_mul(z,z,p)
                mz = mpz_get_str(z)
            end if
        else
            while fn[i]<i do
                k += 1
                integer l = length(factors(k,1))
                if l<=LIMIT and fn[l]<l then
                    fn[l] = iff(fn[l]+1<l?fn[l]+1:k)
                end if
            end while
            mz = sprintf("%d",fn[i])
        end if
    end if
    printf(1,"%3d : %s\n",{i,mz})
end for
printf(1,"completed in %s\n",{elapsed(time()-t0)})
Output:
The first 100 terms in the sequence are:
  1 : 1
  2 : 3
  3 : 25
  4 : 14
  5 : 14641
  6 : 44
  7 : 24137569
  8 : 70
  9 : 1089
 10 : 405
 11 : 819628286980801
 12 : 160
 13 : 22563490300366186081
 14 : 2752
 15 : 9801
 16 : 462
 17 : 21559177407076402401757871041
 18 : 1044
 19 : 740195513856780056217081017732809
 20 : 1520
 21 : 141376
 22 : 84992
 23 : 1658509762573818415340429240403156732495289
 24 : 1170
 25 : 52200625
 26 : 421888
 27 : 52900
 28 : 9152
 29 : 1116713952456127112240969687448211536647543601817400964721
 30 : 6768
 31 : 1300503809464370725741704158412711229899345159119325157292552449
 32 : 3990
 33 : 12166144
 34 : 9764864
 35 : 446265625
 36 : 5472
 37 : 11282036144040442334289838466416927162302790252609308623697164994458730076798801
 38 : 43778048
 39 : 90935296
 40 : 10416
 41 : 1300532588674810624476094551095787816112173600565095470117230812218524514342511947837104801
 42 : 46400
 43 : 635918448514386699807643535977466343285944704172890141356181792680152445568879925105775366910081
 44 : 240640
 45 : 327184
 46 : 884998144
 47 : 82602452843197830915655434062758747152610200533183747995128511868250464749389571755574391210629602061883161
 48 : 10296
 49 : 17416274304961
 50 : 231984
 51 : 3377004544
 52 : 1175552
 53 : 7326325566540660915295202005885275873916026034616342139474905237555535331121749053330837020397976615915057535109963186790081
 54 : 62208
 55 : 382260265984
 56 : 63168
 57 : 18132238336
 58 : 74356621312
 59 : 4611334279555550707926152839105934955536765902552873727962394200823974159354935875908492026570361080937000929065119751494662472171586496615769
 60 : 37200
 61 : 1279929743416851311019131209907830943453757487243270654630811620734985849511676634764875391422075025095805774223361200187655617244608064273703030801
 62 : 329638739968
 63 : 4000000
 64 : 41160
 65 : 6169143218176
 66 : 1446912
 67 : 20353897784481135224502113429729640062994484338530413467091588021107086251737634020247647652000753728181181145357697865506347474542010115076391004870941216126804332281
 68 : 22478848
 69 : 505031950336
 70 : 920000
 71 : 22091712217028661091647719716134154062183987922906664635563029317259865249987461330814689139636373404600637581380931231750650949001643115899851798743405544731506806491024751606849
 72 : 48300
 73 : 45285235038445046669368642612544904396805516154393281169675637706411327508046898517381759728413013085702957690245765106506995874808813788844198933536768701568785385215106907990288684161
 74 : 26044681682944
 75 : 25040016
 76 : 103546880
 77 : 6818265813529681
 78 : 6860800
 79 : 110984176612396876252402058909207317796166059426692518840795949938301678339569859458072604697803922487329059012193474923358078243829751108364014428972188856355641430510895584045477184112155202949344511201
 80 : 96720
 81 : 4708900
 82 : 473889511571456
 83 : 1064476683917919713953093000677954858036756167846865592483240200233630032347646244510522542053167377047784795269272961130616738371982635464615430562192693194769301221853619917764723198332349478419665523610384617408161
 84 : 225216
 85 : 629009610244096
 86 : 1974722883485696
 87 : 56062476550144
 88 : 1469440
 89 : 2544962774801294304714624882135254894108219227449639770372304502957346499018390075803907657903246999131414158076182409047363202723848127272231619125736007088495905384436604400674375401897829996007586872027878808309385140119563002941281
 90 : 352512
 91 : 334095024862954369
 92 : 2017460224
 93 : 258858752671744
 94 : 35114003344654336
 95 : 6002585119227904
 96 : 112860
 97 : 69969231567692157576407845029145070949540195647704307603423555494283752374775631665902846216473259715737953596002226233187827382886325202177640164868195792546734599315840795700630834939445407388277880586442087150607690134279001258366485550281200590593848327041
 98 : 22588608
 99 : 226984356
100 : 870000
completed in 4.4s

Python

This implementation exploits the fact that terms corresponding to a prime value for n are always the nth prime to the (n-1)th power.

def divisors(n):
    divs = [1]
    for ii in range(2, int(n ** 0.5) + 3):
        if n % ii == 0:
            divs.append(ii)
            divs.append(int(n / ii))
    divs.append(n)
    return list(set(divs))


def is_prime(n):
    return len(divisors(n)) == 2


def primes():
    ii = 1
    while True:
        ii += 1
        if is_prime(ii):
            yield ii


def prime(n):
    generator = primes()
    for ii in range(n - 1):
        generator.__next__()
    return generator.__next__()


def n_divisors(n):
    ii = 0
    while True:
        ii += 1
        if len(divisors(ii)) == n:
            yield ii


def sequence(max_n=None):
    if max_n is not None:
        for ii in range(1, max_n + 1):
            if is_prime(ii):
                yield prime(ii) ** (ii - 1)
            else:
                generator = n_divisors(ii)
                for jj, out in zip(range(ii - 1), generator):
                    pass
                yield generator.__next__()
    else:
        ii = 1
        while True:
            ii += 1
            if is_prime(ii):
                yield prime(ii) ** (ii - 1)
            else:
                generator = n_divisors(ii)
                for jj, out in zip(range(ii - 1), generator):
                    pass
                yield generator.__next__()


if __name__ == '__main__':
    for item in sequence(15):
        print(item)

Output:

1
3
25
14
14641
44
24137569
70
1089
405
819628286980801
160
22563490300366186081
2752
9801

Raku

(formerly Perl 6)

Works with: Rakudo version 2019.03

Try it online!

sub div-count (\x) {
    return 2 if x.is-prime;
    +flat (1 .. x.sqrt.floor).map: -> \d {
        unless x % d { my \y = x div d; y == d ?? y !! (y, d) }
    }
}

my $limit = 20;

my @primes = grep { .is-prime }, 1..*;
@primes[$limit]; # prime the array. SCNR

put "First $limit terms of OEIS:A073916";
put (1..$limit).hyper(:2batch).map: -> $n {
    ($n > 4 and $n.is-prime) ??
    exp($n - 1, @primes[$n - 1]) !!
    do {
        my $i = 0;
        my $iterator = $n %% 2 ?? (1..*) !! (1..*).map: *²;
        $iterator.first: {
            next unless $n == .&div-count;
            next unless ++$i == $n;
            $_
        }
    }
};
First 20 terms of OEIS:A073916
1 3 25 14 14641 44 24137569 70 1089 405 819628286980801 160 22563490300366186081 2752 9801 462 21559177407076402401757871041 1044 740195513856780056217081017732809 1520

REXX

Programming note:   this REXX version has minor optimization, and all terms of the sequence are determined (found) in order.

little optimization

/*REXX program  finds and displays  the    Nth  number   with exactly   N   divisors.   */
parse arg N .                                    /*obtain optional argument from the CL.*/
if N=='' | N==","  then N= 15                    /*Not specified?  Then use the default.*/
if N>=50  then numeric digits 10                 /*use more decimal digits for large N. */
w= 50                                            /*W:  width of the 2nd column of output*/
say '─divisors─'  center("the Nth number with exactly N divisors", w, '─')      /*title.*/
@.1= 2;                                   Ps= 1  /*1st prime;  number of primes (so far)*/
        do p=3  until Ps==N                      /* [↓]  gen N primes, store in @ array.*/
        if \isPrime(p)  then iterate;     Ps= Ps + 1;        @.Ps= p
        end   /*gp*/
!.=                                              /*the  !  array is used for memoization*/
        do i=1  for N;      odd= i//2            /*step through a number of divisors.   */
        if odd  then  if isPrime(i)  then do;  _= pPow();            w= max(w, length(_) )
                                               call tell  commas(_);              iterate
                                          end
        #= 0;            even= \odd              /*the number of occurrences for #div.  */
            do j=1;      jj= j                   /*now, search for a number that ≡ #divs*/
            if odd  then jj= j*j                 /*Odd and non-prime?  Calculate square.*/
            if !.jj==.  then iterate             /*has this number already been found?  */
            d= #divs(jj)                         /*get # divisors;  Is not equal?  Skip.*/
            if even  then if d<i  then do;  !.j=.;  iterate;  end   /*Too low?  Flag it.*/
            if d\==i  then iterate               /*Is not equal?  Then skip this number.*/
            #= # + 1                             /*bump number of occurrences for #div. */
            if #\==i  then iterate               /*Not correct occurrence? Keep looking.*/
            call tell  commas(jj)                /*display Nth number with #divs*/
            leave                                /*found a number, so now get the next I*/
            end   /*j*/
        end       /*i*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _;  do j=length(_)-3  to 1  by -3; _=insert(',', _, j); end;    return _
pPow:   numeric digits 1000;  return @.i**(i-1)  /*temporarily increase decimal digits. */
tell: parse arg _; say center(i,10) right(_,max(w,length(_))); if i//5==0 then say; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
#divs: procedure; parse arg x 1 y                /*X and Y:  both set from 1st argument.*/
       if x<7  then do                           /*handle special cases for numbers < 7.*/
                    if x<3   then return x       /*   "      "      "    "  one and two.*/
                    if x<5   then return x - 1   /*   "      "      "    "  three & four*/
                    if x==5  then return 2       /*   "      "      "    "  five.       */
                    if x==6  then return 4       /*   "      "      "    "  six.        */
                    end
       odd= x // 2                               /*check if   X   is  odd  or not.      */
       if odd  then do;  #= 1;             end   /*Odd?   Assume  Pdivisors  count of 1.*/
               else do;  #= 3;    y= x%2;  end   /*Even?     "        "        "    " 3.*/
                                                 /* [↑]   start with known num of Pdivs.*/
                  do k=3  by 1+odd  while k<y    /*when doing odd numbers,  skip evens. */
                  if x//k==0  then do            /*if no remainder, then found a divisor*/
                                   #=#+2;  y=x%k /*bump  #  Pdivs,  calculate limit  Y. */
                                   if k>=y  then do;  #= #-1;  leave;  end      /*limit?*/
                                   end                                          /*  ___ */
                              else if k*k>x  then leave        /*only divide up to √ x  */
                  end   /*k*/                    /* [↑]  this form of DO loop is faster.*/
       return #+1                                /*bump "proper divisors" to "divisors".*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
isPrime: procedure; parse arg #;         if wordpos(#, '2 3 5 7 11 13')\==0  then return 1
         if #<2  then return 0;    if #//2==0 | #//3==0 | #//5==0 | #//7==0  then return 0
                                         if # // 2==0 | # // 3    ==0  then return 0
           do j=11  by 6  until j*j>#;   if # // j==0 | # // (J+2)==0  then return 0
           end   /*j*/                           /*           ___                       */
         return 1                                /*Exceeded  √ #  ?    Then # is prime. */
output   when using the input:     45

(Shown at   3/4   size.)

─divisors─ ───────────────────────────────────────────the Nth number with exactly N divisors──────────────────────────────────────────────
    1                                                                                                                                    1
    2                                                                                                                                    3
    3                                                                                                                                   25
    4                                                                                                                                   14
    5                                                                                                                               14,641

    6                                                                                                                                   44
    7                                                                                                                           24,137,569
    8                                                                                                                                   70
    9                                                                                                                                1,089
    10                                                                                                                                 405

    11                                                                                                                 819,628,286,980,801
    12                                                                                                                                 160
    13                                                                                                          22,563,490,300,366,186,081
    14                                                                                                                               2,752
    15                                                                                                                               9,801

    16                                                                                                                                 462
    17                                                                                              21,559,177,407,076,402,401,757,871,041
    18                                                                                                                               1,044
    19                                                                                         740,195,513,856,780,056,217,081,017,732,809
    20                                                                                                                               1,520

    21                                                                                                                             141,376
    22                                                                                                                              84,992
    23                                                                           1,658,509,762,573,818,415,340,429,240,403,156,732,495,289
    24                                                                                                                               1,170
    25                                                                                                                          52,200,625

    26                                                                                                                             421,888
    27                                                                                                                              52,900
    28                                                                                                                               9,152
    29                                                       1,116,713,952,456,127,112,240,969,687,448,211,536,647,543,601,817,400,964,721
    30                                                                                                                               6,768

    31                                               1,300,503,809,464,370,725,741,704,158,412,711,229,899,345,159,119,325,157,292,552,449
    32                                                                                                                               3,990
    33                                                                                                                          12,166,144
    34                                                                                                                           9,764,864
    35                                                                                                                         446,265,625

    36                                                                                                                               5,472
    37                          11,282,036,144,040,442,334,289,838,466,416,927,162,302,790,252,609,308,623,697,164,994,458,730,076,798,801
    38                                                                                                                          43,778,048
    39                                                                                                                          90,935,296
    40                                                                                                                              10,416

    41           1,300,532,588,674,810,624,476,094,551,095,787,816,112,173,600,565,095,470,117,230,812,218,524,514,342,511,947,837,104,801
    42                                                                                                                              46,400
    43     635,918,448,514,386,699,807,643,535,977,466,343,285,944,704,172,890,141,356,181,792,680,152,445,568,879,925,105,775,366,910,081
    44                                                                                                                             240,640
    45                                                                                                                             327,184

more optimization

Programming note:   this REXX version has major optimization, and the logic flow is:

  •   build a table of prime numbers (this also helps winnow the numbers being tested).
  •   the generation of the sequence is broken into three parts:
  •   odd prime numbers.
  •   odd non-prime numbers.
  •   even numbers.

This REXX version (unlike the 1st version),   only goes through the numbers once, instead of looking for numbers that have specific number of divisors.

/*REXX program  finds and displays  the    Nth  number   with exactly   N   divisors.   */
parse arg N .                                    /*obtain optional argument from the CL.*/
if N=='' | N==","  then N= 15                    /*Not specified?  Then use the default.*/
if N>=50  then numeric digits 10                 /*use more decimal digits for large N. */
@.1= 2;               Ps= 1;    !.= 0;    !.1= 2 /*1st prime;  number of primes (so far)*/
        do p=3  until Ps==N**3                   /* [↓]  gen N primes, store in @ array.*/
        if \isPrime(p)  then iterate;     Ps= Ps + 1;    if Ps<=N  then  @.Ps= p;   !.p= 1
        end   /*p*/

zfin.= 0;    zcnt. = 0;  znum.1= 1;  znum.2= 3   /*completed;   index;   count of items.*/
w= 50                                            /*──────────handle odd primes──────────*/
     do j=3  by 2  to N;  if \!.j  then iterate  /*Not prime?  Then skip this odd number*/
     zfin.j= 1;   zcnt.j= j;   znum.j= pPow();   /*compute # divisors for this odd prime*/
     w= max(w, length( commas( znum.j) ) )       /*the last prime will be the biggest #.*/
     end   /*j*/                                 /*process a small number of primes ≤ N.*/
dd.=;                     mx= 200000             /*──────────handle odd non─primes──────*/
     do j=3  by 2  to N;  if !.j  then iterate   /*Is a prime?  Then skip this odd prime*/
        do sq=6;  _= sq*sq                       /*step through squares starting at  36.*/
        if dd._\=='' then d= dd._                /*maybe use a pre─computed # divisors. */
                     else d= #divs(_)            /*Not defined?  Then calculate # divs. */
        if _<=mx  then dd._= d                   /*use memoization for the  evens  loop.*/
        if d\==j  then iterate                   /*if not the right D, then skip this sq*/
        zcnt.d= zcnt.d+1;         if zcnt.d==d  then zfin.d= 1;        znum.d= _
        if zfin.d  then iterate j                /*if all were found,  then do next odd#*/
        end   /*sq*/
     end      /*j*/
                                                 /*──────────handle even numbers.───────*/
     do j=4  by 2; if dd.j\=='' then d= dd.j     /*maybe use a pre─computed # divisors. */
                                else d= #divs(j) /*Not defined?  Then calculate # divs. */
     if d>N       then iterate                   /*Divisors greater than N?  Then skip. */
     if zfin.d    then iterate                   /*Already populated?          "    "   */
                  else do; zcnt.d= zcnt.d+1;  if zcnt.d==d  then zfin.d= 1;  znum.d= j
                           if done()  then leave  /*j*/    /*Are the even #'s all done? */
                       end
     end       /*j*/

say '─divisors─'  center("the Nth number with exactly N divisors", w, '─')      /*title.*/
     do s=1  for N;  call tell  s,commas(znum.s) /*display  Nth  number with number divs*/
     end   /*s*/
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg _;  do c=length(_)-3  to 1  by -3; _=insert(',', _, c); end;    return _
done:      do f=N  by -1  for N-3;      if \zfin.f  then return 0;        end;    return 1
pPow:   numeric digits 2000;  return @.j**(j-1)  /*temporarily increase decimal digits. */
tell: parse arg _; say center(i,10) right(_,max(w,length(_))); if i//5==0 then say; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
#divs: procedure; parse arg x 1 y                /*X and Y:  both set from 1st argument.*/
       if x<7  then do                           /*handle special cases for numbers < 7.*/
                    if x<3   then return x       /*   "      "      "    "  one and two.*/
                    if x<5   then return x - 1   /*   "      "      "    "  three & four*/
                    if x==5  then return 2       /*   "      "      "    "  five.       */
                    if x==6  then return 4       /*   "      "      "    "  six.        */
                    end
       odd= x // 2                               /*check if   X   is  odd  or not.      */
       if odd  then do;  #= 1;             end   /*Odd?   Assume  Pdivisors  count of 1.*/
               else do;  #= 3;    y= x%2;  end   /*Even?     "        "        "    " 3.*/
                                                 /* [↑]   start with known num of Pdivs.*/
                  do k=3  by 1+odd  while k<y    /*when doing odd numbers,  skip evens. */
                  if x//k==0  then do            /*if no remainder, then found a divisor*/
                                   #=#+2;  y=x%k /*bump  #  Pdivs,  calculate limit  Y. */
                                   if k>=y  then do;  #= #-1;  leave;  end      /*limit?*/
                                   end                                          /*  ___ */
                              else if k*k>x  then leave        /*only divide up to √ x  */
                  end   /*k*/                    /* [↑]  this form of DO loop is faster.*/
       return #+1                                /*bump "proper divisors" to "divisors".*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
isPrime: procedure; parse arg # . '' -1 _
         if #<31  then do;   if wordpos(#, '2 3 5 7 11 13 17 19 23 29')\==0  then return 1
                             if #<2  then return 0
                       end
         if #// 2==0 then return 0; if #// 3==0  then return 0; if     _==5  then return 0
         if #// 7==0 then return 0; if #//11==0  then return 0; if #//11==0  then return 0
         if #//13==0 then return 0; if #//17==0  then return 0; if #//19==0  then return 0
                               do i=23  by 6  until i*i>#;   if #// i   ==0  then return 0
                                                             if #//(i+2)==0  then return 0
                               end   /*i*/       /*           ___                       */
         return 1                                /*Exceeded  √ #  ?    Then # is prime. */
output   is identical to the 1st REXX version.


Ring

load "stdlib.ring"

num = 0
limit = 22563490300366186081

see "working..." + nl
see "the first 15 terms of the sequence are:" + nl

for n = 1 to 15
    num = 0
    for m = 1 to limit
        pnum = 0
        for p = 1 to limit
            if (m % p = 0)
               pnum = pnum + 1
            ok
        next
        if pnum = n
           num = num + 1
           if num = n
              see "" + n + ": " + m + " " + nl
              exit 
           ok 
        ok
     next
next

see nl + "done..." + nl
Output:
working...
the first 15 terms of the sequence are:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14641
 6 : 44
 7 : 24137569
 8 : 70
 9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
done...

Ruby

Translation of: Java
def isPrime(n)
    return false if n < 2
    return n == 2 if n % 2 == 0
    return n == 3 if n % 3 == 0

    k = 5
    while k * k <= n
        return false if n % k == 0
        k = k + 2
    end

    return true
end

def getSmallPrimes(numPrimes)
    smallPrimes = [2]
    count = 0
    n = 3
    while count < numPrimes
        if isPrime(n) then
            smallPrimes << n
            count = count + 1
        end
        n = n + 2
    end
    return smallPrimes
end

def getDivisorCount(n)
    count = 1
    while n % 2 == 0
        n = (n / 2).floor
        count = count + 1
    end

    d = 3
    while d * d <= n
        q = (n / d).floor
        r = n % d
        dc = 0
        while r == 0
            dc = dc + count
            n = q
            q = (n / d).floor
            r = n % d
        end
        count = count + dc
        d = d + 2
    end
    if n != 1 then
        count = 2 * count
    end
    return count
end

MAX = 15
@smallPrimes = getSmallPrimes(MAX)

def OEISA073916(n)
    if isPrime(n) then
        return @smallPrimes[n - 1] ** (n - 1)
    end

    count = 0
    result = 0
    i = 1
    while count < n
        if n % 2 == 1 then
            # The solution for an odd (non-prime) term is always a square number
            root = Math.sqrt(i)
            if root * root != i then
                i = i + 1
                next
            end
        end
        if getDivisorCount(i) == n then
            count = count + 1
            result = i
        end
        i = i + 1
    end
    return result
end

n = 1
while n <= MAX
    print "A073916(", n, ") = ", OEISA073916(n), "\n"
    n = n + 1
end
Output:
A073916(1) = 1
A073916(2) = 3
A073916(3) = 25
A073916(4) = 14
A073916(5) = 14641
A073916(6) = 44
A073916(7) = 24137569
A073916(8) = 70
A073916(9) = 1089
A073916(10) = 405
A073916(11) = 819628286980801
A073916(12) = 160
A073916(13) = 22563490300366186081
A073916(14) = 2752
A073916(15) = 9801

Sidef

func f(n {.is_prime}) {
    n.prime**(n-1)
}

func f(n) {
    n.th { .sigma0 == n }
}

say 20.of { f(_+1) }
Output:
[1, 3, 25, 14, 14641, 44, 24137569, 70, 1089, 405, 819628286980801, 160, 22563490300366186081, 2752, 9801, 462, 21559177407076402401757871041, 1044, 740195513856780056217081017732809, 1520]

Wren

Translation of: Kotlin
Library: Wren-math
Library: Wren-big
Library: Wren-fmt
import "./math" for Int
import "./big" for BigInt
import "./fmt" for Fmt

var MAX = 33
var primes = Int.primeSieve(MAX * 5)
System.print("The first %(MAX) terms in the sequence are:")
for (i in 1..MAX) {
    if (Int.isPrime(i)) {
        var z = BigInt.new(primes[i-1]).pow(i-1)
        Fmt.print("$2d : $i", i, z)
    } else {
        var count = 0
        var j = 1
        while (true) {
            var cont = false
            if (i % 2 == 1) {
                var sq = j.sqrt.floor
                if (sq * sq != j) {
                    j = j + 1
                    cont = true
                }
            }
            if (!cont) {
                if (Int.divisors(j).count == i) {
                    count = count + 1
                    if (count == i) {
                        Fmt.print("$2d : $d", i, j)
                        break
                    }
                }
                j = j + 1
            }
        }
    }
}
Output:
The first 33 terms in the sequence are:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14641
 6 : 44
 7 : 24137569
 8 : 70
 9 : 1089
10 : 405
11 : 819628286980801
12 : 160
13 : 22563490300366186081
14 : 2752
15 : 9801
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170
25 : 52200625
26 : 421888
27 : 52900
28 : 9152
29 : 1116713952456127112240969687448211536647543601817400964721
30 : 6768
31 : 1300503809464370725741704158412711229899345159119325157292552449
32 : 3990
33 : 12166144

zkl

Translation of: Go

Using GMP (GNU Multiple Precision Arithmetic Library, probabilistic primes), because it is easy and fast to generate primes.

Extensible prime generator#zkl could be used instead.

var [const] BI=Import("zklBigNum"), pmax=25;  // libGMP
p:=BI(1);
primes:=pmax.pump(List(0), p.nextPrime, "copy");  //-->(0,3,5,7,11,13,17,19,...)
 
fcn countDivisors(n){
   count:=1;
   while(n%2==0){ n/=2; count+=1; }
   foreach d in ([3..*,2]){
      q,r := n/d, n%d;
      if(r==0){
	 dc:=0;
	 while(r==0){
	    dc+=count;
	    n,q,r = q, n/d, n%d;
	 }
	 count+=dc;
      }
      if(d*d > n) break;
   }
   if(n!=1) count*=2;
   count
}

println("The first ", pmax, " terms in the sequence are:");
foreach i in ([1..pmax]){
   if(BI(i).probablyPrime()) println("%2d : %,d".fmt(i,primes[i].pow(i-1)));
   else{
      count:=0;
      foreach j in ([1..*]){
         if(i%2==1 and j != j.toFloat().sqrt().toInt().pow(2)) continue;
	 if(countDivisors(j) == i){
	    count+=1;
	    if(count==i){
	       println("%2d : %,d".fmt(i,j));
	       break;
	    }
	 }
      }
   }
}
Output:
The first 25 terms in the sequence are:
 1 : 1
 2 : 3
 3 : 25
 4 : 14
 5 : 14,641
 6 : 44
 7 : 24,137,569
 8 : 70
 9 : 1,089
10 : 405
11 : 819,628,286,980,801
12 : 160
13 : 22,563,490,300,366,186,081
14 : 2,752
15 : 9,801
16 : 462
17 : 21,559,177,407,076,402,401,757,871,041
18 : 1,044
19 : 740,195,513,856,780,056,217,081,017,732,809
20 : 1,520
21 : 141,376
22 : 84,992
23 : 1,658,509,762,573,818,415,340,429,240,403,156,732,495,289
24 : 1,170
25 : 52,200,625