Sequence: nth number with exactly n divisors: Difference between revisions
m (→{{header|Wren}}: Minor tidy) |
|||
(5 intermediate revisions by 4 users not shown) | |||
Line 14: | Line 14: | ||
:*[[Sequence: smallest number greater than previous term with exactly n divisors]] |
:*[[Sequence: smallest number greater than previous term with exactly n divisors]] |
||
:*[[Sequence: smallest number with exactly n divisors]] |
:*[[Sequence: smallest number with exactly n divisors]] |
||
<br/><br/> |
|||
=={{header|ALGOL 68}}== |
|||
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}} |
|||
Uses Algol 68G's LONG LONG INT (with default precision) to handle the large numbers.<br/> |
|||
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.<br/> |
|||
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.<br/> |
|||
For non-prime n, the divisors are also shown. |
|||
<syntaxhighlight lang="algol68"> |
|||
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 |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre style="font-size: 12px"> |
|||
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</pre> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
{{trans|C++}} |
{{trans|C++}} |
||
< |
<syntaxhighlight lang="c">#include <math.h> |
||
#include <stdbool.h> |
#include <stdbool.h> |
||
#include <stdint.h> |
#include <stdint.h> |
||
Line 135: | Line 210: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>A073916(1) = 1 |
<pre>A073916(1) = 1 |
||
Line 155: | Line 230: | ||
=={{header|C++}}== |
=={{header|C++}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang="cpp">#include <iostream> |
||
#include <vector> |
#include <vector> |
||
Line 246: | Line 321: | ||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>A073916(1) = 1 |
<pre>A073916(1) = 1 |
||
Line 266: | Line 341: | ||
=={{header|D}}== |
=={{header|D}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang="d">import std.bigint; |
||
import std.math; |
import std.math; |
||
import std.stdio; |
import std.stdio; |
||
Line 355: | Line 430: | ||
writeln("A073916(", n, ") = ", OEISA073916(n)); |
writeln("A073916(", n, ") = ", OEISA073916(n)); |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>A073916(1) = 1 |
<pre>A073916(1) = 1 |
||
Line 405: | Line 480: | ||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
This makes use of most of the optimizations discussed in the Go example. |
This makes use of most of the optimizations discussed in the Go example. |
||
< |
<syntaxhighlight lang="factor">USING: combinators formatting fry kernel lists lists.lazy |
||
lists.lazy.examples literals math math.functions math.primes |
lists.lazy.examples literals math math.functions math.primes |
||
math.primes.factors math.ranges sequences ; |
math.primes.factors math.ranges sequences ; |
||
Line 436: | Line 511: | ||
: main ( -- ) 45 [1,b] [ dup fn "%2d : %d\n" printf ] each ; |
: main ( -- ) 45 [1,b] [ dup fn "%2d : %d\n" printf ] each ; |
||
MAIN: main</ |
MAIN: main</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 485: | Line 560: | ||
45 : 327184 |
45 : 327184 |
||
</pre> |
</pre> |
||
=={{header|FreeBASIC}}== |
|||
But it takes an "eternity" to resolve. |
|||
<syntaxhighlight lang="freebasic">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</syntaxhighlight> |
|||
=={{header|Go}}== |
=={{header|Go}}== |
||
Line 490: | Line 590: | ||
The remaining terms (up to the 33rd) are not particularly large and so are calculated by brute force. |
The remaining terms (up to the 33rd) are not particularly large and so are calculated by brute force. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 572: | Line 672: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 620: | Line 720: | ||
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. |
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. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 740: | Line 840: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 792: | Line 892: | ||
</pre> |
</pre> |
||
=={{header|Haskell}}== |
=={{header|Haskell}}== |
||
< |
<syntaxhighlight lang="haskell">import Control.Monad (guard) |
||
import Math.NumberTheory.ArithmeticFunctions (divisorCount) |
import Math.NumberTheory.ArithmeticFunctions (divisorCount) |
||
import Math.NumberTheory.Primes (Prime, unPrime) |
import Math.NumberTheory.Primes (Prime, unPrime) |
||
Line 819: | Line 919: | ||
main :: IO () |
main :: IO () |
||
main = mapM_ print nths</ |
main = mapM_ print nths</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>(1,1) |
<pre>(1,1) |
||
Line 856: | Line 956: | ||
(34,9764864) |
(34,9764864) |
||
(35,446265625)</pre> |
(35,446265625)</pre> |
||
=={{header|J}}== |
|||
Implementation: |
|||
<syntaxhighlight lang="j"> |
|||
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 |
|||
</syntaxhighlight> |
|||
Task example: |
|||
<syntaxhighlight lang="j"> |
|||
(,. 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 |
|||
</syntaxhighlight> |
|||
=={{header|Java}}== |
=={{header|Java}}== |
||
Line 861: | Line 1,008: | ||
Replace translation with Java native implementation. |
Replace translation with Java native implementation. |
||
< |
<syntaxhighlight lang="java"> |
||
import java.math.BigInteger; |
import java.math.BigInteger; |
||
import java.util.ArrayList; |
import java.util.ArrayList; |
||
Line 950: | Line 1,097: | ||
} |
} |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
Line 1,014: | Line 1,161: | ||
'''Preliminaries''' |
'''Preliminaries''' |
||
< |
<syntaxhighlight lang="jq">def count(stream): reduce stream as $i (0; .+1); |
||
# To maintain precision: |
# To maintain precision: |
||
Line 1,037: | Line 1,184: | ||
end |
end |
||
end; |
end; |
||
</syntaxhighlight> |
|||
</lang> |
|||
'''The Task''' |
'''The Task''' |
||
< |
<syntaxhighlight lang="jq"># Emit [n, nth_with_n_divisors] for n in range(1; .+1) |
||
def nth_with_n_divisors: |
def nth_with_n_divisors: |
||
| [limit( .; primes)] as $primes |
| [limit( .; primes)] as $primes |
||
Line 1,065: | Line 1,212: | ||
"The first 33 terms in the sequence are:", |
"The first 33 terms in the sequence are:", |
||
(33 | nth_with_n_divisors)</ |
(33 | nth_with_n_divisors)</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,105: | Line 1,252: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
< |
<syntaxhighlight lang="julia">using Primes |
||
function countdivisors(n) |
function countdivisors(n) |
||
Line 1,134: | Line 1,281: | ||
nthwithndivisors(35) |
nthwithndivisors(35) |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
1 : 1 |
1 : 1 |
||
Line 1,175: | Line 1,322: | ||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
< |
<syntaxhighlight lang="scala">// Version 1.3.21 |
||
import java.math.BigInteger |
import java.math.BigInteger |
||
Line 1,253: | Line 1,400: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{output}} |
{{output}} |
||
Line 1,294: | Line 1,441: | ||
=={{header|Mathematica}} / {{header|Wolfram Language}}== |
=={{header|Mathematica}} / {{header|Wolfram Language}}== |
||
< |
<syntaxhighlight lang="mathematica">d = Table[ |
||
Length[Divisors[n]], {n, |
Length[Divisors[n]], {n, |
||
200000}]; t = {}; n = 0; ok = True; While[ok, n++; |
200000}]; t = {}; n = 0; ok = True; While[ok, n++; |
||
Line 1,300: | Line 1,447: | ||
c = Flatten[Position[d, n, 1, n]]; |
c = Flatten[Position[d, n, 1, n]]; |
||
If[Length[c] >= n, AppendTo[t, c[[n]]], ok = False]]]; |
If[Length[c] >= n, AppendTo[t, c[[n]]], ok = False]]]; |
||
t</ |
t</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>{1,3,25,14,14641,44,24137569,70,1089,405,819628286980801,160,22563490300366186081,2752,9801,462,21559177407076402401757871041,1044,740195513856780056217081017732809,1520,141376,84992,1658509762573818415340429240403156732495289,1170}</pre> |
<pre>{1,3,25,14,14641,44,24137569,70,1089,405,819628286980801,160,22563490300366186081,2752,9801,462,21559177407076402401757871041,1044,740195513856780056217081017732809,1520,141376,84992,1658509762573818415340429240403156732495289,1170}</pre> |
||
Line 1,308: | Line 1,455: | ||
{{libheader|bignum}} |
{{libheader|bignum}} |
||
This is a translation of the fast Go version. It runs in about 23s on our laptop. |
This is a translation of the fast Go version. It runs in about 23s on our laptop. |
||
< |
<syntaxhighlight lang="nim">import math, strformat |
||
import bignum |
import bignum |
||
Line 1,386: | Line 1,533: | ||
k > records[cd].num and (d == 1 or d == 2 and not cd.isOdd): |
k > records[cd].num and (d == 1 or d == 2 and not cd.isOdd): |
||
records[cd].num = k |
records[cd].num = k |
||
inc records[cd].count</ |
inc records[cd].count</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,439: | Line 1,586: | ||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
||
{{trans|Raku}} |
{{trans|Raku}} |
||
< |
<syntaxhighlight lang="perl">use strict; |
||
use warnings; |
use warnings; |
||
use bigint; |
use bigint; |
||
Line 1,459: | Line 1,606: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>First 20 terms of OEIS:A073916 |
<pre>First 20 terms of OEIS:A073916 |
||
Line 1,469: | Line 1,616: | ||
Certainly not the fastest way to do it, hence the relatively small limit of 24, which takes less than 0.4s,<br> |
Certainly not the fastest way to do it, hence the relatively small limit of 24, which takes less than 0.4s,<br> |
||
whereas a limit of 25 would need to invoke factors() 52 million times which would no doubt take a fair while. |
whereas a limit of 25 would need to invoke factors() 52 million times which would no doubt take a fair while. |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
||
<span style="color: #008080;">constant</span> <span style="color: #000000;">LIMIT</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">24</span> |
<span style="color: #008080;">constant</span> <span style="color: #000000;">LIMIT</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">24</span> |
||
Line 1,493: | Line 1,640: | ||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,525: | Line 1,672: | ||
===cheating slightly=== |
===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). |
No real patterns that I could see here, but you can still identify and single out the troublemakers (of which there are about 30). |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
||
Line 1,590: | Line 1,737: | ||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"completed in %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)})</span> |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"completed in %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)})</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 1,699: | Line 1,846: | ||
=={{header|Python}}== |
=={{header|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. |
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. |
||
<syntaxhighlight lang="python"> |
|||
<lang Python> |
|||
def divisors(n): |
def divisors(n): |
||
divs = [1] |
divs = [1] |
||
Line 1,763: | Line 1,910: | ||
for item in sequence(15): |
for item in sequence(15): |
||
print(item) |
print(item) |
||
</syntaxhighlight> |
|||
</lang> |
|||
<b>Output:</b> |
<b>Output:</b> |
||
<syntaxhighlight lang="python"> |
|||
<lang Python> |
|||
1 |
1 |
||
3 |
3 |
||
Line 1,781: | Line 1,928: | ||
2752 |
2752 |
||
9801 |
9801 |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
Line 1,789: | Line 1,936: | ||
[https://tio.run/##dVLbTsJAEH3nKw4IpgW6ATQY2QAao4kvmsijEFPoVpr05nZr2hh@yk/wx@p0CwUf3JfuzJwzc/ZMYyH9cVEk6RqO92ltojRUMJaZia8G6EihUhliBM9FxrzEiqUXCK5rPde3CTwEY1RLPqRirh9F0mSBHU9gzbB09m3Kk4a@SBJk6IDSCHIsc0wppsFwOCiYUmU@p1uzCSPvwzGx0/xdg74NorR9L/AU0UYDXmVutKKEUu9SxNS4VoldH0PGuryxx7xW7BXHGSqE2grYUto5w@Lu6YU6xqlC68GTiTrMUkIGCSIXz/ePi8nt4OriejhucY00qH8FM9k2j4U0JqO1rTbbowftcO@BQbcZLmGHDiVrlSa9WNdFFpcQC8M@asE6XplkiMY40YmhpR0evXvA/6ZIsK0iSRWidzq0PPK0VNo1tbH6Vunr/nwfyTWTueX7JyejyhOKTB2WSI1pWey8/mf4v9BerxRZavmLab/VYb1jXhS/ Try it online!] |
[https://tio.run/##dVLbTsJAEH3nKw4IpgW6ATQY2QAao4kvmsijEFPoVpr05nZr2hh@yk/wx@p0CwUf3JfuzJwzc/ZMYyH9cVEk6RqO92ltojRUMJaZia8G6EihUhliBM9FxrzEiqUXCK5rPde3CTwEY1RLPqRirh9F0mSBHU9gzbB09m3Kk4a@SBJk6IDSCHIsc0wppsFwOCiYUmU@p1uzCSPvwzGx0/xdg74NorR9L/AU0UYDXmVutKKEUu9SxNS4VoldH0PGuryxx7xW7BXHGSqE2grYUto5w@Lu6YU6xqlC68GTiTrMUkIGCSIXz/ePi8nt4OriejhucY00qH8FM9k2j4U0JqO1rTbbowftcO@BQbcZLmGHDiVrlSa9WNdFFpcQC8M@asE6XplkiMY40YmhpR0evXvA/6ZIsK0iSRWidzq0PPK0VNo1tbH6Vunr/nwfyTWTueX7JyejyhOKTB2WSI1pWey8/mf4v9BerxRZavmLab/VYb1jXhS/ Try it online!] |
||
<lang |
<syntaxhighlight lang="raku" line>sub div-count (\x) { |
||
return 2 if x.is-prime; |
return 2 if x.is-prime; |
||
+flat (1 .. x.sqrt.floor).map: -> \d { |
+flat (1 .. x.sqrt.floor).map: -> \d { |
||
Line 1,814: | Line 1,961: | ||
} |
} |
||
} |
} |
||
};</ |
};</syntaxhighlight> |
||
<pre>First 20 terms of OEIS:A073916 |
<pre>First 20 terms of OEIS:A073916 |
||
Line 1,822: | Line 1,969: | ||
Programming note: this REXX version has minor optimization, and all terms of the sequence are determined (found) in order. |
Programming note: this REXX version has minor optimization, and all terms of the sequence are determined (found) in order. |
||
===little optimization=== |
===little optimization=== |
||
< |
<syntaxhighlight lang="rexx">/*REXX program finds and displays the Nth number with exactly N divisors. */ |
||
parse arg N . /*obtain optional argument from the CL.*/ |
parse arg N . /*obtain optional argument from the CL.*/ |
||
if N=='' | N=="," then N= 15 /*Not specified? Then use the default.*/ |
if N=='' | N=="," then N= 15 /*Not specified? Then use the default.*/ |
||
Line 1,881: | Line 2,028: | ||
do j=11 by 6 until j*j>#; if # // j==0 | # // (J+2)==0 then return 0 |
do j=11 by 6 until j*j>#; if # // j==0 | # // (J+2)==0 then return 0 |
||
end /*j*/ /* ___ */ |
end /*j*/ /* ___ */ |
||
return 1 /*Exceeded √ # ? Then # is prime. */</ |
return 1 /*Exceeded √ # ? Then # is prime. */</syntaxhighlight> |
||
{{out|output|text= when using the input: <tt> 45 </tt>}} |
{{out|output|text= when using the input: <tt> 45 </tt>}} |
||
Line 1,952: | Line 2,099: | ||
This REXX version (unlike the 1<sup>st</sup> version), only goes through the numbers once, instead |
This REXX version (unlike the 1<sup>st</sup> version), only goes through the numbers once, instead |
||
of looking for numbers that have specific number of divisors. |
of looking for numbers that have specific number of divisors. |
||
< |
<syntaxhighlight lang="rexx">/*REXX program finds and displays the Nth number with exactly N divisors. */ |
||
parse arg N . /*obtain optional argument from the CL.*/ |
parse arg N . /*obtain optional argument from the CL.*/ |
||
if N=='' | N=="," then N= 15 /*Not specified? Then use the default.*/ |
if N=='' | N=="," then N= 15 /*Not specified? Then use the default.*/ |
||
Line 2,028: | Line 2,175: | ||
if #//(i+2)==0 then return 0 |
if #//(i+2)==0 then return 0 |
||
end /*i*/ /* ___ */ |
end /*i*/ /* ___ */ |
||
return 1 /*Exceeded √ # ? Then # is prime. */</ |
return 1 /*Exceeded √ # ? Then # is prime. */</syntaxhighlight> |
||
{{out|output|text= is identical to the 1<sup>st</sup> REXX version.}} <br><br> |
{{out|output|text= is identical to the 1<sup>st</sup> REXX version.}} <br><br> |
||
=={{header|Ring}}== |
=={{header|Ring}}== |
||
< |
<syntaxhighlight lang="ring"> |
||
load "stdlib.ring" |
load "stdlib.ring" |
||
Line 2,061: | Line 2,208: | ||
see nl + "done..." + nl |
see nl + "done..." + nl |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,086: | Line 2,233: | ||
=={{header|Ruby}}== |
=={{header|Ruby}}== |
||
{{trans|Java}} |
{{trans|Java}} |
||
< |
<syntaxhighlight lang="ruby">def isPrime(n) |
||
return false if n < 2 |
return false if n < 2 |
||
return n == 2 if n % 2 == 0 |
return n == 2 if n % 2 == 0 |
||
Line 2,174: | Line 2,321: | ||
print "A073916(", n, ") = ", OEISA073916(n), "\n" |
print "A073916(", n, ") = ", OEISA073916(n), "\n" |
||
n = n + 1 |
n = n + 1 |
||
end</ |
end</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>A073916(1) = 1 |
<pre>A073916(1) = 1 |
||
Line 2,193: | Line 2,340: | ||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
< |
<syntaxhighlight lang="ruby">func f(n {.is_prime}) { |
||
n.prime**(n-1) |
n.prime**(n-1) |
||
} |
} |
||
Line 2,201: | Line 2,348: | ||
} |
} |
||
say 20.of { f(_+1) }</ |
say 20.of { f(_+1) }</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 2,212: | Line 2,359: | ||
{{libheader|Wren-big}} |
{{libheader|Wren-big}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
< |
<syntaxhighlight lang="wren">import "./math" for Int |
||
import "/big" for BigInt |
import "./big" for BigInt |
||
import "/fmt" for Fmt |
import "./fmt" for Fmt |
||
var MAX = 33 |
var MAX = 33 |
||
Line 2,247: | Line 2,394: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 2,293: | Line 2,440: | ||
[[Extensible prime generator#zkl]] could be used instead. |
[[Extensible prime generator#zkl]] could be used instead. |
||
< |
<syntaxhighlight lang="zkl">var [const] BI=Import("zklBigNum"), pmax=25; // libGMP |
||
p:=BI(1); |
p:=BI(1); |
||
primes:=pmax.pump(List(0), p.nextPrime, "copy"); //-->(0,3,5,7,11,13,17,19,...) |
primes:=pmax.pump(List(0), p.nextPrime, "copy"); //-->(0,3,5,7,11,13,17,19,...) |
||
Line 2,332: | Line 2,479: | ||
} |
} |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
Latest revision as of 11:36, 4 February 2024
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
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
#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++
#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
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
// 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
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
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
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)
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
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
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
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