Sequence: nth number with exactly n divisors

From Rosetta Code
Revision as of 01:16, 12 April 2019 by rosettacode>Gerard Schildberger (→‎{{header|REXX}}: elided a dead statement.)
Sequence: nth number with exactly n divisors 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.

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

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.

The remaining terms (up to the 24th) are not particularly large and so are calculated by brute force. <lang go>package main

import (

   "fmt"
   "math/big"

)

const max = 24

var primes = [max]int{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,

   41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89}

func isSmallPrime(n int) bool {

   for _, prime := range primes {
       if prime == n {
           return true
       }
   }
   return false

}

func countDivisors(n int) int {

   count := 0
   for i := 1; i*i <= n; i++ {
       if n%i == 0 {
           if i == n/i {
               count++
           } else {
               count += 2
           }
       }
   }
   return count

}

func main() {

   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 isSmallPrime(i) {
           z.SetInt64(int64(primes[i-1]))
           p.SetInt64(int64(i - 1))
           z.Exp(z, p, nil)
           fmt.Printf("%2d : %d\n", i, z)
       } else {
           count := 0
           for j := 1; ; j++ {
               if countDivisors(j) == i {
                   count++
                   if count == i {
                       fmt.Printf("%2d : %d\n", i, j)
                       break
                   }
               }
           }
       }
   }

}</lang>

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

Java

Translation of: Go

<lang java>import java.math.BigInteger;

public class OEIS_A073916 {

   static int[] primes = new int[]{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89};
   static boolean is_small_prime(int n) {
       for (int prime : primes) {
           if (prime == n) return true;
       }
       return false;
   }
   static int count_divisors(int n) {
       int count = 0;
       for (int i = 1; i * i <= n; ++i) {
           if (n % i == 0) {
               if (i == n / i)
                   count++;
               else
                   count += 2;
           }
       }
       return count;
   }
   public static void main(String[] args) {
       final int max = 24;
       System.out.printf("The first %d terms of the sequence are:\n", max);
       for (int i = 1; i <= max; ++i) {
           if (is_small_prime(i)) {
               BigInteger z = BigInteger.valueOf(primes[i - 1]);
               z = z.pow(i - 1);
               System.out.printf("%2d : %d\n", i, z);
           } else {
               for (int j = 1, count = 0; ; ++j) {
                   if (count_divisors(j) == i) {
                       if (++count == i) {
                           System.out.printf("%2d : %d\n", i, j);
                           break;
                       }
                   }
               }
           }
       }
   }

}</lang>

Output:
The first 24 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
16 : 462
17 : 21559177407076402401757871041
18 : 1044
19 : 740195513856780056217081017732809
20 : 1520
21 : 141376
22 : 84992
23 : 1658509762573818415340429240403156732495289
24 : 1170

Kotlin

Translation of: Go

<lang scala>// Version 1.3.21

import java.math.BigInteger

const val MAX = 24

var primes = intArrayOf(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,

   41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89)

fun isSmallPrime(n: Int) = primes.contains(n)

fun countDivisors(n: Int): Int {

   var count = 0
   var i = 1
   while (i * i <= n) {
       if (n % i == 0) {
           count += if (i == n / i) 1 else 2
       }
       i++
   }
   return count

}

fun main() {

   println("The first $MAX terms in the sequence are:")
   for (i in 1..MAX) {
       if (isSmallPrime(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 (countDivisors(j) == i) {
                   if (++count == i) {
                       System.out.printf("%2d : %d\n", i, j)
                       break
                   }
               }
               j++
           }
       }
   }

}</lang>

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

Perl 6

Works with: Rakudo version 2019.03

Try it online!

<lang perl6>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;
           $_
       }
   }

};</lang>

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 program automatically right justifies the output.
If the output is wider than   100   decimal digits (with commas),   it is left justified. <lang rexx>/*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.*/ say '─divisors─' center("the Nth number with exactly N divisors", 100, '─') /*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;   say center(i, 10)  pPow();       iterate
                                         end
       #= 0                                     /*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); if d\==i  then iterate /*get # divisors;  Is not equal?  Skip.*/
           !.jj=.                               /*mark as having found #divs for this J*/
           #= # + 1                             /*bump number of occurrences for #div. */
           if #\==i  then iterate               /*Not correct occurrence? Keep looking.*/
           say center(i, 10)   right( commas(jj), 100)  /*display Nth number with #divs*/
           if i//4==0  then say                         /*display a separator for eyes.*/
           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 500; _= commas(@.i ** (i-1) ); return right(_, max(100, length(_))) /*──────────────────────────────────────────────────────────────────────────────────────*/

  1. 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  for x%2-3  by 1+odd  while k<y  /*for 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. */</lang>
output   when using the input:     33
─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

Sidef

<lang ruby>func f(n {.is_prime}) {

   n.prime**(n-1)

}

func f(n) {

   n.th { .sigma0 == n }

}

say 20.of { f(_+1) }</lang>

Output:
[1, 3, 25, 14, 14641, 44, 24137569, 70, 1089, 405, 819628286980801, 160, 22563490300366186081, 2752, 9801, 462, 21559177407076402401757871041, 1044, 740195513856780056217081017732809, 1520]