Motzkin numbers
The nth Motzkin number (denoted by M[n]) is the number of different ways of drawing non-intersecting chords between n points on a circle (not necessarily touching every point by a chord).
You are encouraged to solve this task according to the task description, using any language you may know.
- Definition
By convention M[0] = 1.
- Task
Compute and show on this page the first 42 Motzkin numbers or, if your language does not support 64 bit integers, as many such numbers as you can. Indicate which of these numbers are prime.
- See also
- oeis:A001006 Motzkin numbers
ALGOL 68
Uses ALGOL 68G's LONG LONG INT which allows for programmer specifiable precision integers. The default precision is sufficient for this task.
BEGIN # find some Motzkin numbers #
PR read "primes.incl.a68" PR
# returns a table of the Motzkin numbers 0 .. n #
OP MOTZKIN = ( INT n )[]LONG LONG INT:
BEGIN
[ 0 : n ]LONG LONG INT m;
IF n >= 0 THEN
m[ 0 ] := 1;
IF n >= 1 THEN
m[ 1 ] := 1;
FOR i FROM 2 TO UPB m DO
m[ i ] := ( ( m[ i - 1 ] * ( ( 2 * i ) + 1 ) )
+ ( m[ i - 2 ] * ( ( 3 * i ) - 3 ) )
)
OVER ( i + 2 )
OD
FI
FI;
m
END # MOTZKIN # ;
# returns a string representation of n with commas #
PROC commatise = ( LONG LONG INT n )STRING:
BEGIN
STRING result := "";
STRING unformatted = whole( n, 0 );
INT ch count := 0;
FOR c FROM UPB unformatted BY -1 TO LWB unformatted DO
IF ch count <= 2 THEN ch count +:= 1
ELSE ch count := 1; "," +=: result
FI;
unformatted[ c ] +=: result
OD;
result
END; # commatise #
# left-pads a string to at least n characters #
PROC pad left = ( STRING s, INT n )STRING:
IF INT len = ( UPB s - LWB s ) + 1; len >= n THEN s ELSE ( ( n - len ) * " " ) + s FI;
BEGIN # show the Motzkin numbers #
print( ( " n M[n] Prime?", newline ) );
print( ( "-----------------------------------", newline ) );
[]LONG LONG INT m = MOTZKIN 41;
FOR i FROM LWB m TO UPB m DO
print( ( whole( i, -2 )
, pad left( commatise( m[ i ] ), 26 )
, IF is probably prime( m[ i ] ) THEN " prime" ELSE "" FI
, newline
)
)
OD
END
END- Output:
n M[n] Prime? ----------------------------------- 0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2,188 11 5,798 12 15,511 prime 13 41,835 14 113,634 15 310,572 16 853,467 17 2,356,779 18 6,536,382 19 18,199,284 20 50,852,019 21 142,547,559 22 400,763,223 23 1,129,760,415 24 3,192,727,797 25 9,043,402,501 26 25,669,818,476 27 73,007,772,802 28 208,023,278,209 29 593,742,784,829 30 1,697,385,471,211 31 4,859,761,676,391 32 13,933,569,346,707 33 40,002,464,776,083 34 114,988,706,524,270 35 330,931,069,469,828 36 953,467,954,114,363 prime 37 2,750,016,719,520,991 38 7,939,655,757,745,265 39 22,944,749,046,030,949 40 66,368,199,913,921,497 41 192,137,918,101,841,817
Arturo
motzkin: function [n][
result: array.of:n+1 0
result\[0]: 1
result\[1]: 1
loop 2..n 'i ->
result\[i]: ((result\[i-1] * (inc 2*i)) + result\[i-2] * (sub 3*i 3)) / i+2
return result
]
printLine: function [items][
pads: [4 20 7]
loop.with:'i items 'item [
prints pad to :string item pads\[i]
]
print ""
]
motzkins: motzkin 41
printLine ["n" "M[n]" "prime"]
print repeat "=" 31
loop.with:'i motzkins 'm ->
printLine @[i m prime? m]
- Output:
n M[n] prime =============================== 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2188 false 11 5798 false 12 15511 true 13 41835 false 14 113634 false 15 310572 false 16 853467 false 17 2356779 false 18 6536382 false 19 18199284 false 20 50852019 false 21 142547559 false 22 400763223 false 23 1129760415 false 24 3192727797 false 25 9043402501 false 26 25669818476 false 27 73007772802 false 28 208023278209 false 29 593742784829 false 30 1697385471211 false 31 4859761676391 false 32 13933569346707 false 33 40002464776083 false 34 114988706524270 false 35 330931069469828 false 36 953467954114363 true 37 2750016719520991 false 38 7939655757745265 false 39 22944749046030949 false 40 66368199913921497 false 41 192137918101841817 false
AWK
# syntax: GAWK --bignum -f MOTZKIN_NUMBERS.AWK
BEGIN {
print(" n Motzkin[n] prime")
limit = 41
m[0] = m[1] = 1
for (i=2; i<=limit; i++) {
m[i] = (m[i-1]*(2*i+1) + m[i-2]*(3*i-3)) / (i + 2)
}
for (i=0; i<=limit; i++) {
printf("%2d %18d %3d\n",i,m[i],is_prime(m[i]))
}
exit(0)
}
function is_prime(x, i) {
if (x <= 1) {
return(0)
}
for (i=2; i<=int(sqrt(x)); i++) {
if (x % i == 0) {
return(0)
}
}
return(1)
}
- Output:
n Motzkin[n] prime 0 1 0 1 1 0 2 2 1 3 4 0 4 9 0 5 21 0 6 51 0 7 127 1 8 323 0 9 835 0 10 2188 0 11 5798 0 12 15511 1 13 41835 0 14 113634 0 15 310572 0 16 853467 0 17 2356779 0 18 6536382 0 19 18199284 0 20 50852019 0 21 142547559 0 22 400763223 0 23 1129760415 0 24 3192727797 0 25 9043402501 0 26 25669818476 0 27 73007772802 0 28 208023278209 0 29 593742784829 0 30 1697385471211 0 31 4859761676391 0 32 13933569346707 0 33 40002464776083 0 34 114988706524270 0 35 330931069469828 0 36 953467954114363 1 37 2750016719520991 0 38 7939655757745265 0 39 22944749046030949 0 40 66368199913921497 0 41 192137918101841817 0
BASIC256
dim M(42)
M[0] = 1 : M[1] = 1
print rjust("1",18) : print rjust("1",18)
for n = 2 to 41
M[n] = M[n-1]
for i = 0 to n-2
M[n] += M[i]*M[n-2-i]
next i
print rjust(string(M[n]),18); chr(9);
if isPrime(M[n]) then print "is a prime" else print
next n
end
function isPrime(v)
if v <= 1 then return False
for i = 2 to int(sqr(v))
if v % i = 0 then return False
next i
return True
end function
Burlesque
41rz{{2./}c!rz2?*{nr}c!\/2?/{JJ2.*J#Rnr#r\/+.nr.-}m[?*++it+.}m[Jp^#R{~]}5E!{fCL[1==}f[#r- Output:
1
1
2
4
9
21
51
127
323
835
2188
5798
15511
41835
113634
310572
853467
2356779
6536382
18199284
50852019
142547559
400763223
1129760415
3192727797
9043402501
25669818476
73007772802
208023278209
593742784829
1697385471211
4859761676391
13933569346707
40002464776083
114988706524270
330931069469828
953467954114363
2750016719520991
7939655757745265
22944749046030949
66368199913921497
192137918101841817
{2 127 15511 953467954114363}
C
#include <locale.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
bool is_prime(uint64_t n) {
if (n < 2)
return false;
if (n % 2 == 0)
return n == 2;
if (n % 3 == 0)
return n == 3;
for (uint64_t p = 5; p * p <= n; p += 4) {
if (n % p == 0)
return false;
p += 2;
if (n % p == 0)
return false;
}
return true;
}
int main() {
setlocale(LC_ALL, "");
printf(" n M(n) Prime?\n");
printf("-----------------------------------\n");
uint64_t m0 = 1, m1 = 1;
for (uint64_t i = 0; i < 42; ++i) {
uint64_t m =
i > 1 ? (m1 * (2 * i + 1) + m0 * (3 * i - 3)) / (i + 2) : 1;
printf("%2llu%'25llu %s\n", i, m, is_prime(m) ? "true" : "false");
m0 = m1;
m1 = m;
}
}
- Output:
n M(n) Prime? ----------------------------------- 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2,188 false 11 5,798 false 12 15,511 true 13 41,835 false 14 113,634 false 15 310,572 false 16 853,467 false 17 2,356,779 false 18 6,536,382 false 19 18,199,284 false 20 50,852,019 false 21 142,547,559 false 22 400,763,223 false 23 1,129,760,415 false 24 3,192,727,797 false 25 9,043,402,501 false 26 25,669,818,476 false 27 73,007,772,802 false 28 208,023,278,209 false 29 593,742,784,829 false 30 1,697,385,471,211 false 31 4,859,761,676,391 false 32 13,933,569,346,707 false 33 40,002,464,776,083 false 34 114,988,706,524,270 false 35 330,931,069,469,828 false 36 953,467,954,114,363 true 37 2,750,016,719,520,991 false 38 7,939,655,757,745,265 false 39 22,944,749,046,030,949 false 40 66,368,199,913,921,497 false 41 192,137,918,101,841,817 false
C#
Arbitrary precision. Now showing results up to the 80th. Note comment about square root limit, one might think the limit should be a long instead of an int, since the square root of a 35 digit number could be up to 18 digits (too big for an int).
using System;
using BI = System.Numerics.BigInteger;
class Program {
// has multiple factors (other than 1 and x)
static bool hmf(BI x) {
if (x < 4) return x == 1;
if ((x & 1) == 0 || x % 3 == 0) return true;
int l = (int)Math.Sqrt((double)x); // this limit works because the 40th to 80th Motzkin numbers have factors of 2 or 3
for (int j = 5, d = 4; j <= l; j += d = 6 - d)
if (x % j == 0) return x > j;
return false;
}
static void Main(string[] args) {
BI a = 0, b = 1, t;
int n = 1, s = 0, d = 1, c = 0, f = 1;
while (n <= 80)
Console.WriteLine("{0,46:n0} {1}",
t = b / n++,
hmf(t) ? "" : "is prime.",
t = b,
b = ((c += d * 3 + 3) * a +
(f += d * 2 + 3) * b) /
(s += d += 2),
a = t);
}
}
- Output:
1
1
2 is prime.
4
9
21
51
127 is prime.
323
835
2,188
5,798
15,511 is prime.
41,835
113,634
310,572
853,467
2,356,779
6,536,382
18,199,284
50,852,019
142,547,559
400,763,223
1,129,760,415
3,192,727,797
9,043,402,501
25,669,818,476
73,007,772,802
208,023,278,209
593,742,784,829
1,697,385,471,211
4,859,761,676,391
13,933,569,346,707
40,002,464,776,083
114,988,706,524,270
330,931,069,469,828
953,467,954,114,363 is prime.
2,750,016,719,520,991
7,939,655,757,745,265
22,944,749,046,030,949
66,368,199,913,921,497
192,137,918,101,841,817
556,704,809,728,838,604
1,614,282,136,160,911,722
4,684,478,925,507,420,069
13,603,677,110,519,480,289
39,532,221,379,621,112,004
114,956,499,435,014,161,638
334,496,473,194,459,009,429
973,899,740,488,107,474,693
2,837,208,756,709,314,025,578
8,270,140,811,590,103,129,028
24,119,587,499,879,368,045,581
70,380,687,801,729,972,163,737
205,473,381,836,953,330,090,977
600,161,698,382,141,668,958,313
1,753,816,895,177,229,449,263,803
5,127,391,665,653,918,424,581,931
14,996,791,899,280,244,858,336,604
43,881,711,243,248,048,262,611,670
128,453,535,912,993,825,479,057,919
376,166,554,620,363,320,971,336,899
1,101,997,131,244,113,831,001,323,618
3,229,547,920,421,385,142,120,565,580
9,468,017,265,749,942,384,739,441,267
27,766,917,351,255,946,264,000,229,811
81,459,755,507,915,876,737,297,376,646
239,056,762,740,830,735,069,669,439,852
701,774,105,036,927,170,410,592,656,651
2,060,763,101,398,061,220,299,787,957,807
6,053,261,625,552,368,838,017,538,638,577
17,785,981,695,172,350,686,294,020,499,397
52,274,487,460,035,748,810,950,928,411,209
153,681,622,703,766,437,645,990,598,724,233
451,929,928,113,276,686,826,984,901,736,388
1,329,334,277,731,700,374,912,787,442,584,082
3,911,184,337,415,864,255,099,077,969,308,357
11,510,402,374,965,653,734,436,362,305,721,089
33,882,709,435,158,403,490,429,948,661,355,518
99,762,777,233,730,236,158,474,945,885,114,348
C++
#include <cstdint>
#include <iomanip>
#include <iostream>
bool is_prime(uint64_t n) {
if (n < 2)
return false;
if (n % 2 == 0)
return n == 2;
if (n % 3 == 0)
return n == 3;
for (uint64_t p = 5; p * p <= n; p += 4) {
if (n % p == 0)
return false;
p += 2;
if (n % p == 0)
return false;
}
return true;
}
class motzkin_generator {
public:
uint64_t next();
private:
uint64_t n = 0;
uint64_t m0 = 1;
uint64_t m1 = 1;
};
uint64_t motzkin_generator::next() {
uint64_t m = n > 1 ? (m1 * (2 * n + 1) + m0 * (3 * n - 3)) / (n + 2) : 1;
++n;
m0 = m1;
m1 = m;
return m;
}
int main() {
std::cout.imbue(std::locale(""));
std::cout << " n M(n) Prime?\n";
std::cout << "-----------------------------------\n";
std::cout << std::boolalpha;
motzkin_generator mgen;
for (int i = 0; i < 42; ++i) {
uint64_t n = mgen.next();
std::cout << std::setw(2) << i << std::setw(25) << n << " "
<< is_prime(n) << '\n';
}
}
- Output:
n M(n) Prime? ----------------------------------- 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2,188 false 11 5,798 false 12 15,511 true 13 41,835 false 14 113,634 false 15 310,572 false 16 853,467 false 17 2,356,779 false 18 6,536,382 false 19 18,199,284 false 20 50,852,019 false 21 142,547,559 false 22 400,763,223 false 23 1,129,760,415 false 24 3,192,727,797 false 25 9,043,402,501 false 26 25,669,818,476 false 27 73,007,772,802 false 28 208,023,278,209 false 29 593,742,784,829 false 30 1,697,385,471,211 false 31 4,859,761,676,391 false 32 13,933,569,346,707 false 33 40,002,464,776,083 false 34 114,988,706,524,270 false 35 330,931,069,469,828 false 36 953,467,954,114,363 true 37 2,750,016,719,520,991 false 38 7,939,655,757,745,265 false 39 22,944,749,046,030,949 false 40 66,368,199,913,921,497 false 41 192,137,918,101,841,817 false
Dart
import 'dart:math';
void main() {
var M = List<int>.filled(42, 1);
M[0] = 1;
M[1] = 1;
print('1');
print('1');
for (int n = 2; n < 42; ++n) {
M[n] = M[n - 1];
for (int i = 0; i <= n - 2; ++i) {
M[n] += M[i] * M[n - 2 - i];
}
if (isPrime(M[n])) {
print('${M[n]} is a prime');
} else {
print('${M[n]}');
}
}
}
bool isPrime(int n) {
if (n <= 1) return false;
if (n == 2) return true;
for (int i = 2; i <= sqrt(n); ++i) {
if (n % i == 0) return false;
}
return true;
}
- Output:
1 1 2 is a prime 4 9 21 51 127 is a prime 323 835 2188 5798 15511 is a prime 41835 113634 310572 853467 2356779 6536382 18199284 50852019 142547559 400763223 1129760415 3192727797 9043402501 25669818476 73007772802 208023278209 593742784829 1697385471211 4859761676391 13933569346707 40002464776083 114988706524270 330931069469828 953467954114363 is a prime 2750016719520991 7939655757745265 22944749046030956 66368199913921500 192137918101841820
EasyLang
fastfunc isprim num .
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
m[] = [ 1 1 ]
print 1 ; print 1
len m[] 39
arrbase m[] 0
for n = 2 to 38
m[n] = m[n - 1]
for i = 0 to n - 2
m[n] += m[i] * m[n - 2 - i]
.
if isprim m[n] = 1
print m[n] & " is a prime"
else
print m[n]
.
.EMal
fun isPrime = logic by int n
if n <= 1 do return false end
for int i = 2; i <= int!sqrt(n); ++i
if n % i == 0 do return false end
end
return true
end
fun format = text by var n, int max do return " " * (max - length(text!n)) + n end
fun motzkin = List by int n
List m = int[].with(n)
m[0] = 1
m[1] = 1
for int i = 2; i < m.length; ++i
m[i] = (m[i - 1] * (2 * i + 1) + m[i - 2] * (3 * i - 3)) / (i + 2)
end
return m
end
int n = 0
writeLine(format("n", 2) + format("motzkin(n)", 20) + format("prime", 6))
for each int m in motzkin(42)
writeLine(format(n, 2) + format(m, 20) + format(isPrime(m), 4))
++n
end- Output:
n motzkin(n) prime 0 1 ⊥ 1 1 ⊥ 2 2 ⊤ 3 4 ⊥ 4 9 ⊥ 5 21 ⊥ 6 51 ⊥ 7 127 ⊤ 8 323 ⊥ 9 835 ⊥ 10 2188 ⊥ 11 5798 ⊥ 12 15511 ⊤ 13 41835 ⊥ 14 113634 ⊥ 15 310572 ⊥ 16 853467 ⊥ 17 2356779 ⊥ 18 6536382 ⊥ 19 18199284 ⊥ 20 50852019 ⊥ 21 142547559 ⊥ 22 400763223 ⊥ 23 1129760415 ⊥ 24 3192727797 ⊥ 25 9043402501 ⊥ 26 25669818476 ⊥ 27 73007772802 ⊥ 28 208023278209 ⊥ 29 593742784829 ⊥ 30 1697385471211 ⊥ 31 4859761676391 ⊥ 32 13933569346707 ⊥ 33 40002464776083 ⊥ 34 114988706524270 ⊥ 35 330931069469828 ⊥ 36 953467954114363 ⊤ 37 2750016719520991 ⊥ 38 7939655757745265 ⊥ 39 22944749046030949 ⊥ 40 66368199913921497 ⊥ 41 192137918101841817 ⊥
F#
// Motzkin numbers. Nigel Galloway: September 10th., 2021
let M=let rec fN g=seq{yield List.item 1 g; yield! fN(0L::(g|>List.windowed 3|>List.map(List.sum))@[0L;0L])} in fN [0L;1L;0L;0L]
M|>Seq.take 42|>Seq.iter(printfn "%d")
- Output:
1 1 2 4 9 21 51 127 323 835 2188 5798 15511 41835 113634 310572 853467 2356779 6536382 18199284 50852019 142547559 400763223 1129760415 3192727797 9043402501 25669818476 73007772802 208023278209 593742784829 1697385471211 4859761676391 13933569346707 40002464776083 114988706524270 330931069469828 953467954114363 2750016719520991 7939655757745265 22944749046030949 66368199913921497 192137918101841817
Factor
USING: combinators formatting io kernel math math.primes
tools.memory.private ;
MEMO: motzkin ( m -- n )
dup 2 < [
drop 1
] [
{
[ 2 * 1 + ]
[ 1 - motzkin * ]
[ 3 * 3 - ]
[ 2 - motzkin * + ]
[ 2 + /i ]
} cleave
] if ;
" n motzkin(n)\n" print
42 [
dup motzkin [ commas ] keep prime? "prime" "" ?
"%2d %24s %s\n" printf
] each-integer
- Output:
n motzkin(n) 0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2,188 11 5,798 12 15,511 prime 13 41,835 14 113,634 15 310,572 16 853,467 17 2,356,779 18 6,536,382 19 18,199,284 20 50,852,019 21 142,547,559 22 400,763,223 23 1,129,760,415 24 3,192,727,797 25 9,043,402,501 26 25,669,818,476 27 73,007,772,802 28 208,023,278,209 29 593,742,784,829 30 1,697,385,471,211 31 4,859,761,676,391 32 13,933,569,346,707 33 40,002,464,776,083 34 114,988,706,524,270 35 330,931,069,469,828 36 953,467,954,114,363 prime 37 2,750,016,719,520,991 38 7,939,655,757,745,265 39 22,944,749,046,030,949 40 66,368,199,913,921,497 41 192,137,918,101,841,817
Fermat
Array m[42];
m[1]:=1;
m[2]:=2;
!!(1,0); {precompute and print m[0] thru m[2]}
!!(1,0);
!!(2,1);
for n=3 to 41 do
m[n]:=(m[n-1]*(2*n+1) + m[n-2]*(3*n-3))/(n+2);
!!(m[n],Isprime(m[n]));
od;
; {built-in Isprime function returns 0 for 1, 1 for primes, and}
; {the smallest prime factor for composites, so this actually gives}
; {slightly more information than requested}- Output:
1 0 1 0 2 1 4 2 9 3 21 3 51 3 127 1 323 17 835 5 2188 2 5798 2 15511 1 41835 3 113634 2 310572 2 853467 3 2356779 3 6536382 2 18199284 2 50852019 3 142547559 3 400763223 3 1129760415 3 3192727797 3 9043402501 7 25669818476 2 73007772802 2 208023278209 107 593742784829 31 1697385471211 1483 4859761676391 3 13933569346707 3 40002464776083 3 114988706524270 2 330931069469828 2 953467954114363 1 2750016719520991 1601 7939655757745265 5 22944749046030949 1063 66368199913921497 3 192137918101841817 3
FreeBASIC
Use any of the primality testing examples as an include.
#include "isprime.bas"
dim as ulongint M(0 to 41)
M(0) = 1 : M(1) = 1
print "1" : print "1"
for n as integer = 2 to 41
M(n) = M(n-1)
for i as uinteger = 0 to n-2
M(n) += M(i)*M(n-2-i)
next i
print M(n),
if isprime(M(n)) then print "is a prime" else print
next n
- Output:
1 1 2 is a prime 4 9 21 51 127 is a prime 323 835 2188 5798 15511 is a prime 41835 113634 310572 853467 2356779 6536382 18199284 50852019 142547559 400763223 1129760415 3192727797 9043402501 25669818476 73007772802 208023278209 593742784829 1697385471211 4859761676391 13933569346707 40002464776083 114988706524270 330931069469828 953467954114363 is a prime 2750016719520991 7939655757745265 22944749046030949 66368199913921497 192137918101841817
FutureBasic
local fn IsPrime( n as NSUInteger ) as BOOL
BOOL isPrime = YES
NSUInteger i
if n < 2 then exit fn = NO
if n = 2 then exit fn = YES
if n mod 2 == 0 then exit fn = NO
for i = 3 to int(n^.5) step 2
if n mod i == 0 then exit fn = NO
next
end fn = isPrime
local fn IsMotzkin
NSUInteger M(42) : M(0) = 1 : M(1) = 1
NSInteger i, n
printf @" 0.%20ld\n 1.%20ld", 1, 1
for n = 2 to 41
M(n) = M(n-1)
for i = 0 to n-2
M(n) += M(i) * M(n-2-i)
next
printf @"%2ld.%20ld\b", n, M(n)
if fn IsPrime( M(n) ) then print " <-- is a prime" else print
next
end fn
fn IsMotzkin
HandleEvents- Output:
0. 1 1. 1 2. 2 <-- is a prime 3. 4 4. 9 5. 21 6. 51 7. 127 <-- is a prime 8. 323 9. 835 10. 2188 11. 5798 12. 15511 <-- is a prime 13. 41835 14. 113634 15. 310572 16. 853467 17. 2356779 18. 6536382 19. 18199284 20. 50852019 21. 142547559 22. 400763223 23. 1129760415 24. 3192727797 25. 9043402501 26. 25669818476 27. 73007772802 28. 208023278209 29. 593742784829 30. 1697385471211 31. 4859761676391 32. 13933569346707 33. 40002464776083 34. 114988706524270 35. 330931069469828 36. 953467954114363 <-- is a prime 37. 2750016719520991 38. 7939655757745265 39. 22944749046030949 40. 66368199913921497 41. 192137918101841817
Go
Assumes a 64 bit Go compiler.
package main
import (
"fmt"
"rcu"
)
func motzkin(n int) []int {
m := make([]int, n+1)
m[0] = 1
m[1] = 1
for i := 2; i <= n; i++ {
m[i] = (m[i-1]*(2*i+1) + m[i-2]*(3*i-3)) / (i + 2)
}
return m
}
func main() {
fmt.Println(" n M[n] Prime?")
fmt.Println("-----------------------------------")
m := motzkin(41)
for i, e := range m {
fmt.Printf("%2d %23s %t\n", i, rcu.Commatize(e), rcu.IsPrime(e))
}
}
- Output:
n M[n] Prime? ----------------------------------- 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2,188 false 11 5,798 false 12 15,511 true 13 41,835 false 14 113,634 false 15 310,572 false 16 853,467 false 17 2,356,779 false 18 6,536,382 false 19 18,199,284 false 20 50,852,019 false 21 142,547,559 false 22 400,763,223 false 23 1,129,760,415 false 24 3,192,727,797 false 25 9,043,402,501 false 26 25,669,818,476 false 27 73,007,772,802 false 28 208,023,278,209 false 29 593,742,784,829 false 30 1,697,385,471,211 false 31 4,859,761,676,391 false 32 13,933,569,346,707 false 33 40,002,464,776,083 false 34 114,988,706,524,270 false 35 330,931,069,469,828 false 36 953,467,954,114,363 true 37 2,750,016,719,520,991 false 38 7,939,655,757,745,265 false 39 22,944,749,046,030,949 false 40 66,368,199,913,921,497 false 41 192,137,918,101,841,817 false
Haskell
import Control.Monad.Memo (Memo, memo, startEvalMemo)
import Math.NumberTheory.Primes.Testing (isPrime)
import System.Environment (getArgs)
import Text.Printf (printf)
type I = Integer
-- The n'th Motzkin number, where n is assumed to be ≥ 0. We memoize the
-- computations using MonadMemo.
motzkin :: I -> Memo I I I
motzkin 0 = return 1
motzkin 1 = return 1
motzkin n = do
m1 <- memo motzkin (n-1)
m2 <- memo motzkin (n-2)
return $ ((2*n+1)*m1 + (3*n-3)*m2) `div` (n+2)
-- The first n+1 Motzkin numbers, starting at 0.
motzkins :: I -> [I]
motzkins = startEvalMemo . mapM motzkin . enumFromTo 0
-- For i = 0 to n print i, the i'th Motzkin number, and if it's a prime number.
printMotzkins :: I -> IO ()
printMotzkins n = mapM_ prnt $ zip [0 :: I ..] (motzkins n)
where prnt (i, m) = printf "%2d %20d %s\n" i m $ prime m
prime m = if isPrime m then "prime" else ""
main :: IO ()
main = do
[n] <- map read <$> getArgs
printMotzkins n
- Output:
0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2188 11 5798 12 15511 prime 13 41835 14 113634 15 310572 16 853467 17 2356779 18 6536382 19 18199284 20 50852019 21 142547559 22 400763223 23 1129760415 24 3192727797 25 9043402501 26 25669818476 27 73007772802 28 208023278209 29 593742784829 30 1697385471211 31 4859761676391 32 13933569346707 33 40002464776083 34 114988706524270 35 330931069469828 36 953467954114363 prime 37 2750016719520991 38 7939655757745265 39 22944749046030949 40 66368199913921497 41 192137918101841817
J
Implementation:
nextMotzkin=: , (({:*1+2*#) + _2&{*3*#-1:)%2+#
Task example:
(":,.' ',.('prime'#~ 1&p:@{.)"1) ,.nextMotzkin^:40(1 1x)
1
1
2 prime
4
9
21
51
127 prime
323
835
2188
5798
15511 prime
41835
113634
310572
853467
2356779
6536382
18199284
50852019
142547559
400763223
1129760415
3192727797
9043402501
25669818476
73007772802
208023278209
593742784829
1697385471211
4859761676391
13933569346707
40002464776083
114988706524270
330931069469828
953467954114363 prime
2750016719520991
7939655757745265
22944749046030949
66368199913921497
192137918101841817
This might not be a particularly interesting to describe algorithm (which might be why the task description currently fails to provide any details of the algorithm). Still... it's probably worth an attempt:
This is an inductive sequence on integers with X0=1 and X1=1 and for the general case of Xn where n>1, (2+n)Xn = ((1+2n)Xn-1+3(n-1)Xn-2).
So basically we calculate (Xn-1(1+2n)+3Xn-2(n-1)) and divide that by 2+n and append this new value to the list. For n we use the list length. For Xn-1 we use the last element of the list. And, for Xn-2 we use the second to last element of the list. For the task we repeat this list operation 40 times, starting with the list 1 1 and check to see which elements of the resulting list are prime. Because these values get large, we need to use arbitrary precision integers for our list values.
Java
import java.math.BigInteger;
import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
public final class MotzkinNumbers {
public static void main(String[] aArgs) {
final int count = 42;
List<BigInteger> motzkins = motzkinNumbers(count);
NumberFormat ukNumberFormat = NumberFormat.getInstance(Locale.UK);
System.out.println(" n Motzkin[n]");
System.out.println("-----------------------------");
for ( int n = 0; n < count; n++ ) {
BigInteger motzkin = motzkins.get(n);
boolean prime = motzkin.isProbablePrime(PROBABILITY);
System.out.print(String.format("%2d %23s", n, ukNumberFormat.format(motzkin)));
System.out.println( prime ? " prime" : "" );
}
}
private static List<BigInteger> motzkinNumbers(int aSize) {
List<BigInteger> result = new ArrayList<BigInteger>(aSize);
result.add(BigInteger.ONE);
result.add(BigInteger.ONE);
for ( int i = 2; i < aSize; i++ ) {
BigInteger nextMotzkin = result.get(i - 1).multiply(BigInteger.valueOf(2 * i + 1))
.add(result.get(i - 2).multiply(BigInteger.valueOf(3 * i - 3)))
.divide(BigInteger.valueOf(i + 2));
result.add(nextMotzkin);
}
return result;
}
private static final int PROBABILITY = 20;
}
- Output:
n Motzkin[n] ----------------------------- 0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2,188 11 5,798 12 15,511 prime 13 41,835 14 113,634 15 310,572 16 853,467 17 2,356,779 18 6,536,382 19 18,199,284 20 50,852,019 21 142,547,559 22 400,763,223 23 1,129,760,415 24 3,192,727,797 25 9,043,402,501 26 25,669,818,476 27 73,007,772,802 28 208,023,278,209 29 593,742,784,829 30 1,697,385,471,211 31 4,859,761,676,391 32 13,933,569,346,707 33 40,002,464,776,083 34 114,988,706,524,270 35 330,931,069,469,828 36 953,467,954,114,363 prime 37 2,750,016,719,520,991 38 7,939,655,757,745,265 39 22,944,749,046,030,949 40 66,368,199,913,921,497 41 192,137,918,101,841,817
jq
Works with gojq, the Go implementation of jq (*)
(*) The C implementation of jq only produces accurate results up to and including n == 37.
For a suitable implementation of `is_prime`, see e.g. # Erdős-primes#jq.
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
def motzkin:
. as $n
| reduce range(2; $n+1) as $i (
{m: [1,1]};
.m[$i] = (.m[$i-1] * (2*$i + 1) + .m[$i-2] * (3*$i -3))/($i + 2))
| .m ;
" n M[n] Prime?",
"------------------------------",
(41 | motzkin
| range(0;length) as $i
|"\($i|lpad(2)) \(.[$i]|lpad(20)) \(.[$i]|if is_prime then "prime" else "" end)")- Output:
n M[n] Prime? ------------------------------ 0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2188 11 5798 12 15511 prime 13 41835 14 113634 15 310572 16 853467 17 2356779 18 6536382 19 18199284 20 50852019 21 142547559 22 400763223 23 1129760415 24 3192727797 25 9043402501 26 25669818476 27 73007772802 28 208023278209 29 593742784829 30 1697385471211 31 4859761676391 32 13933569346707 33 40002464776083 34 114988706524270 35 330931069469828 36 953467954114363 prime 37 2750016719520991 38 7939655757745265 39 22944749046030949 40 66368199913921497 41 192137918101841817
Julia
using Primes
function motzkin(N)
m = zeros(Int, N)
m[1] = m[2] = 1
for i in 3:N
m[i] = (m[i - 1] * (2i - 1) + m[i - 2] * (3i - 6)) ÷ (i + 1)
end
return m
end
println(" n M[n] Prime?\n-----------------------------------")
for (i, m) in enumerate(motzkin(42))
println(lpad(i - 1, 2), lpad(m, 20), lpad(isprime(m), 8))
end
- Output:
n M[n] Prime? ----------------------------------- 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2188 false 11 5798 false 12 15511 true 13 41835 false 14 113634 false 15 310572 false 16 853467 false 17 2356779 false 18 6536382 false 19 18199284 false 20 50852019 false 21 142547559 false 22 400763223 false 23 1129760415 false 24 3192727797 false 25 9043402501 false 26 25669818476 false 27 73007772802 false 28 208023278209 false 29 593742784829 false 30 1697385471211 false 31 4859761676391 false 32 13933569346707 false 33 40002464776083 false 34 114988706524270 false 35 330931069469828 false 36 953467954114363 true 37 2750016719520991 false 38 7939655757745265 false 39 22944749046030949 false 40 66368199913921497 false 41 192137918101841817 false
Mathematica /Wolfram Language
Table[With[{o = Hypergeometric2F1[(1 - n)/2, -n/2, 2, 4]}, {n, o, PrimeQ[o]}], {n, 0, 41}] // Grid
- Output:
Alignment lost in copying.
0 1 False 1 1 False 2 2 True 3 4 False 4 9 False 5 21 False 6 51 False 7 127 True 8 323 False 9 835 False 10 2188 False 11 5798 False 12 15511 True 13 41835 False 14 113634 False 15 310572 False 16 853467 False 17 2356779 False 18 6536382 False 19 18199284 False 20 50852019 False 21 142547559 False 22 400763223 False 23 1129760415 False 24 3192727797 False 25 9043402501 False 26 25669818476 False 27 73007772802 False 28 208023278209 False 29 593742784829 False 30 1697385471211 False 31 4859761676391 False 32 13933569346707 False 33 40002464776083 False 34 114988706524270 False 35 330931069469828 False 36 953467954114363 True 37 2750016719520991 False 38 7939655757745265 False 39 22944749046030949 False 40 66368199913921497 False 41 192137918101841817 False
Maxima
There are many formulas for these numbers
motzkin[0]:1$
motzkin[1]:1$
motzkin[n]:=((2*n+1)*motzkin[n-1]+(3*n-3)*motzkin[n-2])/(n+2)$
/* Test cases */
makelist(motzkin[i],i,0,41);
sublist(%,primep);
- Output:
[1,1,2,4,9,21,51,127,323,835,2188,5798,15511,41835,113634,310572,853467,2356779,6536382,18199284,50852019,142547559,400763223,1129760415,3192727797,9043402501,25669818476,73007772802,208023278209,593742784829,1697385471211,4859761676391,13933569346707,40002464776083,114988706524270,330931069469828,953467954114363,2750016719520991,7939655757745265,22944749046030949,66368199913921497,192137918101841817] [2,127,15511,953467954114363]
Nim
import strformat, strutils
func isPrime(n: Positive): bool =
if n == 1: return false
if (n and 1) == 0: return n == 2
if n mod 3 == 0: return n == 3
var d = 5
while d * d <= n:
if n mod d == 0: return false
inc d, 2
if n mod d == 0: return false
inc d, 4
result = true
func motzkin(n: Positive): seq[int64] =
result.setLen(n + 1)
result[0] = 1
result[1] = 1
for i in 2..n:
result[i] = (result[i-1] * (2 * i + 1) + result[i-2] * (3 * i - 3)) div (i + 2)
echo " n M[n] Prime?"
echo "-----------------------------------"
var m = motzkin(41)
for i, val in m:
echo &"{i:2} {insertSep($val):>23} {val.isPrime}"
- Output:
n M[n] Prime? ----------------------------------- 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2_188 false 11 5_798 false 12 15_511 true 13 41_835 false 14 113_634 false 15 310_572 false 16 853_467 false 17 2_356_779 false 18 6_536_382 false 19 18_199_284 false 20 50_852_019 false 21 142_547_559 false 22 400_763_223 false 23 1_129_760_415 false 24 3_192_727_797 false 25 9_043_402_501 false 26 25_669_818_476 false 27 73_007_772_802 false 28 208_023_278_209 false 29 593_742_784_829 false 30 1_697_385_471_211 false 31 4_859_761_676_391 false 32 13_933_569_346_707 false 33 40_002_464_776_083 false 34 114_988_706_524_270 false 35 330_931_069_469_828 false 36 953_467_954_114_363 true 37 2_750_016_719_520_991 false 38 7_939_655_757_745_265 false 39 22_944_749_046_030_949 false 40 66_368_199_913_921_497 false 41 192_137_918_101_841_817 false
PARI/GP
M=vector(41)
M[1]=1
M[2]=2
for(n=3, 41, M[n]=(M[n-1]*(2*n+1) + M[n-2]*(3*n-3))/(n+2))
M=concat([1], M)
for(n=1, 42, print(M[n]," ",isprime(M[n])))Perl
#!/usr/bin/perl
use strict; # https://rosettacode.org/wiki/Motzkin_numbers
use warnings;
use ntheory qw( is_prime );
sub motzkin
{
my $N = shift;
my @m = ( 0, 1, 1 );
for my $i ( 3 .. $N )
{
$m[$i] = ($m[$i - 1] * (2 * $i - 1) + $m[$i - 2] * (3 * $i - 6)) / ($i + 1);
}
return splice @m, 1;
}
print " n M[n]\n";
my $count = 0;
for ( motzkin(42) )
{
printf "%3d%25s %s\n", $count++, s/\B(?=(\d\d\d)+$)/,/gr,
is_prime($_) ? 'prime' : '';
}
- Output:
n M[n] 0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2,188 11 5,798 12 15,511 prime 13 41,835 14 113,634 15 310,572 16 853,467 17 2,356,779 18 6,536,382 19 18,199,284 20 50,852,019 21 142,547,559 22 400,763,223 23 1,129,760,415 24 3,192,727,797 25 9,043,402,501 26 25,669,818,476 27 73,007,772,802 28 208,023,278,209 29 593,742,784,829 30 1,697,385,471,211 31 4,859,761,676,391 32 13,933,569,346,707 33 40,002,464,776,083 34 114,988,706,524,270 35 330,931,069,469,828 36 953,467,954,114,363 prime 37 2,750,016,719,520,991 38 7,939,655,757,745,265 39 22,944,749,046,030,949 40 66,368,199,913,921,497 41 192,137,918,101,841,817
Phix
with javascript_semantics include mpfr.e function motzkin(integer n) sequence m = mpz_inits(2,1) -- {1,1} mpz m2 = mpz_init() -- (scratch) for i=3 to n do mpz m1 = mpz_init() -- (a new mpz rqd for every m[i]) mpz_mul_si(m1,m[i-1],2*i-1) -- (nb i is 1-based) mpz_mul_si(m2,m[i-2],3*i-6) -- "" mpz_add(m1,m1,m2) assert(mpz_fdiv_q_ui(m1,m1,i+1)==0) -- (verify rmdr==0) m &= m1 end for return m end function printf(1," n M[n]\n") printf(1,"---------------------------\n") sequence m = motzkin(42) for i=1 to 42 do string mi = mpz_get_str(m[i],comma_fill:=true), pi = iff(mpz_prime(m[i])?"prime":"") printf(1,"%2d %23s %s\n", {i-1, mi, pi}) end for
- Output:
n M[n] --------------------------- 0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2,188 11 5,798 12 15,511 prime 13 41,835 14 113,634 15 310,572 16 853,467 17 2,356,779 18 6,536,382 19 18,199,284 20 50,852,019 21 142,547,559 22 400,763,223 23 1,129,760,415 24 3,192,727,797 25 9,043,402,501 26 25,669,818,476 27 73,007,772,802 28 208,023,278,209 29 593,742,784,829 30 1,697,385,471,211 31 4,859,761,676,391 32 13,933,569,346,707 33 40,002,464,776,083 34 114,988,706,524,270 35 330,931,069,469,828 36 953,467,954,114,363 prime 37 2,750,016,719,520,991 38 7,939,655,757,745,265 39 22,944,749,046,030,949 40 66,368,199,913,921,497 41 192,137,918,101,841,817
native
Normally I'd start with this simpler and perfectly valid 64-bit code, but at the moment am somewhat biased towards things that run online, without limiting things as this does.
with javascript_semantics function motzkin(integer n) sequence m = {1,1} for i=3 to n do m &= (m[i-1]*(2*i-1)+m[i-2]*(3*i-6))/(i+1) end for return m end function printf(1," n M[n]\n") printf(1,"---------------------------\n") integer lim = iff(machine_bits()=32?38:42) sequence m = motzkin(lim) for i=1 to lim do string pi = iff(is_prime(m[i])?"prime":"") printf(1,"%2d %,23d %s\n", {i-1, m[i], pi}) end for
Aside: Note that, under 32-bit and p2js, M[36] and M[37] happen only by chance to be correct: an intermediate value exceeds precision (in this case for M[36] it is a multiple of 2 and hence ok, for M[37] it is rounded to the nearest multiple of 4, and the final divide by (i+1) results in an integer simply because there isn't enough precision to hold any fractional part). Output as above on 64-bit, less four entries under 32-bit and pwa/p2js, since unlike M[36..37], M[38..41] don't quite get away with the precision loss, plus M[39] and above exceed precision and hence is_prime() limits on 32-bit anyway.
Python
""" rosettacode.org/wiki/Motzkin_numbers """
from sympy import isprime
def motzkin(num_wanted):
""" Return list of the first N Motzkin numbers """
mot = [1] * (num_wanted + 1)
for i in range(2, num_wanted + 1):
mot[i] = (mot[i-1]*(2*i+1) + mot[i-2]*(3*i-3)) // (i + 2)
return mot
def print_motzkin_table(N=41):
""" Print table of first N Motzkin numbers, and note if prime """
print(
" n M[n] Prime?\n-----------------------------------")
for i, e in enumerate(motzkin(N)):
print(f'{i : 3}{e : 24,}', isprime(e))
print_motzkin_table()
- Output:
Same as Go example.
Quackery
isprime is defined at Primality by trial division#Quackery.
' [ 1 1 ]
41 times
[ dup -2 split nip do
i^ 2 + 2 * 1 - *
swap i^ 2 + 3 * 6 - * +
i^ 3 + / join ]
behead drop
witheach
[ i^ echo sp dup echo isprime if [ say " prime" ] cr ]- Output:
0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2188 11 5798 12 15511 prime 13 41835 14 113634 15 310572 16 853467 17 2356779 18 6536382 19 18199284 20 50852019 21 142547559 22 400763223 23 1129760415 24 3192727797 25 9043402501 26 25669818476 27 73007772802 28 208023278209 29 593742784829 30 1697385471211 31 4859761676391 32 13933569346707 33 40002464776083 34 114988706524270 35 330931069469828 36 953467954114363 prime 37 2750016719520991 38 7939655757745265 39 22944749046030949 40 66368199913921497 41 192137918101841817
Raku
Using binomial coefficients and Catalan numbers
use Lingua::EN::Numbers;
constant \C = 1, |[\×] (2, 6 … ∞) Z/ 2 .. *;
sub binomial { [×] ($^n … 0) Z/ 1 .. $^p }
my \𝐌 = 1, |(1..∞).map: -> \𝐧 { sum ^𝐧 .map( -> \𝐤 { C[𝐤] × binomial 𝐧, 2×𝐤 } ) };
say " 𝐧 𝐌[𝐧] Prime?";
𝐌[^42].kv.map: { printf "%2d %24s %s\n", $^k, $^v.&comma, $v.is-prime };
- Output:
𝐧 𝐌[𝐧] Prime? 0 1 False 1 1 False 2 2 True 3 4 False 4 9 False 5 21 False 6 51 False 7 127 True 8 323 False 9 835 False 10 2,188 False 11 5,798 False 12 15,511 True 13 41,835 False 14 113,634 False 15 310,572 False 16 853,467 False 17 2,356,779 False 18 6,536,382 False 19 18,199,284 False 20 50,852,019 False 21 142,547,559 False 22 400,763,223 False 23 1,129,760,415 False 24 3,192,727,797 False 25 9,043,402,501 False 26 25,669,818,476 False 27 73,007,772,802 False 28 208,023,278,209 False 29 593,742,784,829 False 30 1,697,385,471,211 False 31 4,859,761,676,391 False 32 13,933,569,346,707 False 33 40,002,464,776,083 False 34 114,988,706,524,270 False 35 330,931,069,469,828 False 36 953,467,954,114,363 True 37 2,750,016,719,520,991 False 38 7,939,655,757,745,265 False 39 22,944,749,046,030,949 False 40 66,368,199,913,921,497 False 41 192,137,918,101,841,817 False
Using recurrence relationship
use Lingua::EN::Numbers;
my \𝐌 = 1, 1, { state $i = 2; ++$i; ($^b × (2 × $i - 1) + $^a × (3 × $i - 6)) ÷ ($i + 1) } … *;
say " 𝐧 𝐌[𝐧] Prime?";
𝐌[^42].kv.map: { printf "%2d %24s %s\n", $^k, $^v.&comma, $v.is-prime };
Same output
REXX
/*REXX program to display the first N Motzkin numbers, and if that number is prime. */
numeric digits 92 /*max number of decimal digits for task*/
parse arg n . /*obtain optional argument from the CL.*/
if n=='' | n=="," then n= 42 /*Not specified? Then use the default.*/
w= length(n) + 1; wm= digits()%4 /*define maximum widths for two columns*/
say center('n', w ) center("Motzkin[n]", wm) center(' primality', 11)
say center('' , w, "─") center('' , wm, "─") center('', 11, "─")
@.= 1 /*define default vale for the @. array.*/
do m=0 for n /*step through indices for Motzkin #'s.*/
if m>1 then @.m= (@(m-1)*(m+m+1) + @(m-2)*(m*3-3))%(m+2) /*calculate a Motzkin #*/
call show /*display a Motzkin number ──► terminal*/
end /*m*/
say center('' , w, "─") center('' , wm, "─") center('', 11, "─")
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
@: parse arg i; return @.i /*return function expression based on I*/
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ?
prime: if isPrime(@.m) then return " prime"; return ''
show: y= commas(@.m); say right(m, w) right(y, max(wm, length(y))) prime(); return
/*──────────────────────────────────────────────────────────────────────────────────────*/
isPrime: procedure expose p?. p#. p_.; parse arg x /*persistent P·· REXX variables.*/
if symbol('P?.#')\=='VAR' then /*1st time here? Then define primes. */
do /*L is a list of some low primes < 100.*/
L= 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101
p?.=0 /* [↓] define P_index, P, P squared.*/
do i=1 for words(L); _= word(L,i); p?._= 1; p#.i= _; p_.i= _*_
end /*i*/; p?.0= x2d(3632c8eb5af3b) /*bypass big ÷*/
p?.n= _ + 4 /*define next prime after last prime. */
p?.#= i - 1 /*define the number of primes found. */
end /* p?. p#. p_ must be unique. */
if x<p?.n then return p?.x /*special case, handle some low primes.*/
if x==p?.0 then return 1 /*save a number of primality divisions.*/
parse var x '' -1 _ /*obtain right─most decimal digit of X.*/
if _==5 then return 0; if x//2 ==0 then return 0 /*X ÷ by 5? X ÷ by 2?*/
if x//3==0 then return 0; if x//7 ==0 then return 0 /*" " " 3? " " " 7?*/
/*weed numbers that're ≥ 11 multiples. */
do j=5 to p?.# while p_.j<=x; if x//p#.j ==0 then return 0
end /*j*/
/*weed numbers that're>high multiple Ps*/
do k=p?.n by 6 while k*k<=x; if x//k ==0 then return 0
if x//(k+2)==0 then return 0
end /*k*/; return 1 /*Passed all divisions? It's a prime.*/
- output when using the default input:
n Motzkin[n] primality ─── ─────────────────────── ─────────── 0 1 1 1 2 2 prime 3 4 4 9 5 21 6 51 7 127 prime 8 323 9 835 10 2,188 11 5,798 12 15,511 prime 13 41,835 14 113,634 15 310,572 16 853,467 17 2,356,779 18 6,536,382 19 18,199,284 20 50,852,019 21 142,547,559 22 400,763,223 23 1,129,760,415 24 3,192,727,797 25 9,043,402,501 26 25,669,818,476 27 73,007,772,802 28 208,023,278,209 29 593,742,784,829 30 1,697,385,471,211 31 4,859,761,676,391 32 13,933,569,346,707 33 40,002,464,776,083 34 114,988,706,524,270 35 330,931,069,469,828 36 953,467,954,114,363 prime 37 2,750,016,719,520,991 38 7,939,655,757,745,265 39 22,944,749,046,030,949 40 66,368,199,913,921,497 41 192,137,918,101,841,817 ─── ─────────────────────── ───────────
RPL
It is unfortunately impossible to check M[36] primality without awaking the emulator's watchdog timer.
| RPL code | Comment |
|---|---|
≪
{ #1 #1 } 2 ROT FOR j
DUP DUP SIZE GET j 2 * 1 + *
OVER DUP SIZE 1 - GET j 3 * 3 - * +
j 2 + /
+ NEXT
≫ 'MOTZK' STO
|
MOTZK ( n -- { M(0)..M(n) } )
Loop from 2 to n
( M[i-1]*(2*i+1)
+ M[i-2]*(3*i-3))
/ (i + 2)
add to list
|
- Input:
41 MOTZK
- Output:
1: { # 1d # 1d # 2d # 4d # 9d # 21d # 51d # 127d # 323d # 835d # 2188d # 5798d # 15511d # 41835d # 113634d # 310572d # 853467d # 2356779d # 6536382d # 18199284d # 50852019d # 142547559d # 400763223d # 1129760415d # 3192727797d # 9043402501d # 25669818476d # 73007772802d # 208023278209d # 593742784829d # 1697385471211d # 4859761676391d # 13933569346707d # 40002464776083d # 114988706524270d # 330931069469828d # 953467954114363d # 2750016719520991d # 7939655757745265d # 22944749046030949d # 66368199913921497d # 192137918101841817d }
Since M(n) is fast to compute, we can sacrify execution time to improve code readability, such as:
| RPL code | Comment |
|---|---|
≪
{ #1 #1 } 2 ROT FOR j
DUP DUP SIZE GET LAST SWAP 1 - GET → m1 m2
≪ '(m1*(2*j+1)+m2*(3*j-3))/(j+2)' EVAL ≫
+ NEXT
≫ 'MOTZK' STO
|
MOTZK ( n -- { M(0)..M(n) } )
Loop from 2 to n
get last 2 values of M(n) from the list
get M(n)
add to list
|
Inputs and output are the same.
Ruby
require "prime"
motzkin = Enumerator.new do |y|
m = [1,1]
m.each{|m| y << m }
2.step do |i|
m << (m.last*(2*i+1) + m[-2]*(3*i-3)) / (i+2)
m.unshift # an arr with last 2 elements is sufficient
y << m.last
end
end
motzkin.take(42).each_with_index do |m, i|
puts "#{'%2d' % i}: #{m}#{m.prime? ? ' prime' : ''}"
end
- Output:
0: 1 1: 1 2: 2 prime 3: 4 4: 9 5: 21 6: 51 7: 127 prime 8: 323 9: 835 10: 2188 11: 5798 12: 15511 prime 13: 41835 14: 113634 15: 310572 16: 853467 17: 2356779 18: 6536382 19: 18199284 20: 50852019 21: 142547559 22: 400763223 23: 1129760415 24: 3192727797 25: 9043402501 26: 25669818476 27: 73007772802 28: 208023278209 29: 593742784829 30: 1697385471211 31: 4859761676391 32: 13933569346707 33: 40002464776083 34: 114988706524270 35: 330931069469828 36: 953467954114363 prime 37: 2750016719520991 38: 7939655757745265 39: 22944749046030949 40: 66368199913921497 41: 192137918101841817
Rust
// [dependencies]
// primal = "0.3"
// num-format = "0.4"
fn motzkin(n: usize) -> Vec<usize> {
let mut m = vec![0; n];
m[0] = 1;
m[1] = 1;
for i in 2..n {
m[i] = (m[i - 1] * (2 * i + 1) + m[i - 2] * (3 * i - 3)) / (i + 2);
}
m
}
fn main() {
use num_format::{Locale, ToFormattedString};
let count = 42;
let m = motzkin(count);
println!(" n M(n) Prime?");
println!("-----------------------------------");
for i in 0..count {
println!(
"{:2} {:>23} {}",
i,
m[i].to_formatted_string(&Locale::en),
primal::is_prime(m[i] as u64)
);
}
}
- Output:
n M(n) Prime? ----------------------------------- 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2,188 false 11 5,798 false 12 15,511 true 13 41,835 false 14 113,634 false 15 310,572 false 16 853,467 false 17 2,356,779 false 18 6,536,382 false 19 18,199,284 false 20 50,852,019 false 21 142,547,559 false 22 400,763,223 false 23 1,129,760,415 false 24 3,192,727,797 false 25 9,043,402,501 false 26 25,669,818,476 false 27 73,007,772,802 false 28 208,023,278,209 false 29 593,742,784,829 false 30 1,697,385,471,211 false 31 4,859,761,676,391 false 32 13,933,569,346,707 false 33 40,002,464,776,083 false 34 114,988,706,524,270 false 35 330,931,069,469,828 false 36 953,467,954,114,363 true 37 2,750,016,719,520,991 false 38 7,939,655,757,745,265 false 39 22,944,749,046,030,949 false 40 66,368,199,913,921,497 false 41 192,137,918,101,841,817 false
Sidef
Built-in:
say 50.of { .motzkin }
Motzkin's triangle (the Motzkin numbers are on the hypotenuse of the triangle):
func motzkin_triangle(num, callback) {
var row = []
{ |n|
row = [0, 1, {|i| row[i+2] + row[i] + row[i+1] }.map(0 .. n-3)..., 0]
callback(row.grep{ .> 0 })
} << 2..(num+1)
}
motzkin_triangle(10, {|row|
row.map { "%4s" % _ }.join(' ').say
})
- Output:
1 1 1 1 2 2 1 3 5 4 1 4 9 12 9 1 5 14 25 30 21 1 6 20 44 69 76 51 1 7 27 70 133 189 196 127 1 8 35 104 230 392 518 512 323 1 9 44 147 369 726 1140 1422 1353 835
Task:
func motzkin_numbers(N) {
var (a, b, n) = (0, 1, 1)
N.of {
var M = b//n #/
n += 1
(a, b) = (b, (3*(n-1)*n*a + (2*n - 1)*n*b) // ((n+1)*(n-1))) #/
M
}
}
motzkin_numbers(42).each_kv {|k,v|
say "#{'%2d' % k}: #{v}#{v.is_prime ? ' prime' : ''}"
}
- Output:
0: 1 1: 1 2: 2 prime 3: 4 4: 9 5: 21 6: 51 7: 127 prime 8: 323 9: 835 10: 2188 11: 5798 12: 15511 prime 13: 41835 14: 113634 15: 310572 16: 853467 17: 2356779 18: 6536382 19: 18199284 20: 50852019 21: 142547559 22: 400763223 23: 1129760415 24: 3192727797 25: 9043402501 26: 25669818476 27: 73007772802 28: 208023278209 29: 593742784829 30: 1697385471211 31: 4859761676391 32: 13933569346707 33: 40002464776083 34: 114988706524270 35: 330931069469828 36: 953467954114363 prime 37: 2750016719520991 38: 7939655757745265 39: 22944749046030949 40: 66368199913921497 41: 192137918101841817
Swift
import Foundation
extension BinaryInteger {
@inlinable
public var isPrime: Bool {
if self == 0 || self == 1 {
return false
} else if self == 2 {
return true
}
let max = Self(ceil((Double(self).squareRoot())))
for i in stride(from: 2, through: max, by: 1) where self % i == 0 {
return false
}
return true
}
}
func motzkin(_ n: Int) -> [Int] {
var m = Array(repeating: 0, count: n)
m[0] = 1
m[1] = 1
for i in 2..<n {
m[i] = (m[i - 1] * (2 * i + 1) + m[i - 2] * (3 * i - 3)) / (i + 2)
}
return m
}
let m = motzkin(42)
for (i, n) in m.enumerated() {
print("\(i): \(n) \(n.isPrime ? "prime" : "")")
}
- Output:
0: 1 1: 1 2: 2 prime 3: 4 4: 9 5: 21 6: 51 7: 127 prime 8: 323 9: 835 10: 2188 11: 5798 12: 15511 prime 13: 41835 14: 113634 15: 310572 16: 853467 17: 2356779 18: 6536382 19: 18199284 20: 50852019 21: 142547559 22: 400763223 23: 1129760415 24: 3192727797 25: 9043402501 26: 25669818476 27: 73007772802 28: 208023278209 29: 593742784829 30: 1697385471211 31: 4859761676391 32: 13933569346707 33: 40002464776083 34: 114988706524270 35: 330931069469828 36: 953467954114363 prime 37: 2750016719520991 38: 7939655757745265 39: 22944749046030949 40: 66368199913921497 41: 192137918101841817
Uiua
Due to limited precision in Uiua only the first 36 terms are generated accurately.
# Uses following formula to generate successive terms:
# m[i] = (m[i-2]*(3*i-3) + m[i-1]*(2*i+1)) / (i + 2)
# Collect the multiplicands from the above formula.
Multiplicands ← [⊃(-3×3|+1×2|÷:1+2)]⧻
# Apply formula to generate next Motzkin number.
Mn ← ×⊃(/+↙2|⊢↙¯1)×⊃(⊂:1↙¯2|Multiplicands)
# Generate list of terms to reasonable limit.
⍢(⊂⟜Mn|<36⧻) [1 1]
# Simple (=slow) prime test.
IsPrime ← ⊙◌⟨⟨=1⧻⍢(↘1|≠0◿⊢)⊂:1+2⇡⌊√.|1⟩=2.|0⟩≤1.
⍉⊟⟜∵IsPrime- Output:
╭─
╷ 1 0
1 0
2 1
4 0
9 0
21 0
51 0
127 1
323 0
835 0
2188 0
5798 0
15511 1
41835 0
113634 0
310572 0
853467 0
2356779 0
6536382 0
18199284 0
50852019 0
142547559 0
400763223 0
1129760415 0
3192727797 0
9043402501 0
25669818476 0
73007772802 0
208023278209 0
593742784829 0
1697385471211 0
4859761676391 0
13933569346707 0
40002464776083 0
114988706524270 0
330931069469828 0
╯
Wren
import "./long" for ULong
import "./fmt" for Fmt
var motzkin = Fn.new { |n|
var m = List.filled(n+1, 0)
m[0] = ULong.one
m[1] = ULong.one
for (i in 2..n) {
m[i] = (m[i-1] * (2*i + 1) + m[i-2] * (3*i -3))/(i + 2)
}
return m
}
System.print(" n M[n] Prime?")
System.print("-----------------------------------")
var m = motzkin.call(41)
for (i in 0..41) {
Fmt.print("$2d $,23i $s", i, m[i], m[i].isPrime)
}
- Output:
n M[n] Prime? ----------------------------------- 0 1 false 1 1 false 2 2 true 3 4 false 4 9 false 5 21 false 6 51 false 7 127 true 8 323 false 9 835 false 10 2,188 false 11 5,798 false 12 15,511 true 13 41,835 false 14 113,634 false 15 310,572 false 16 853,467 false 17 2,356,779 false 18 6,536,382 false 19 18,199,284 false 20 50,852,019 false 21 142,547,559 false 22 400,763,223 false 23 1,129,760,415 false 24 3,192,727,797 false 25 9,043,402,501 false 26 25,669,818,476 false 27 73,007,772,802 false 28 208,023,278,209 false 29 593,742,784,829 false 30 1,697,385,471,211 false 31 4,859,761,676,391 false 32 13,933,569,346,707 false 33 40,002,464,776,083 false 34 114,988,706,524,270 false 35 330,931,069,469,828 false 36 953,467,954,114,363 true 37 2,750,016,719,520,991 false 38 7,939,655,757,745,265 false 39 22,944,749,046,030,949 false 40 66,368,199,913,921,497 false 41 192,137,918,101,841,817 false
Yabasic
dim M(41)
M(0) = 1 : M(1) = 1
print str$(M(0), "%19g")
print str$(M(1), "%19g")
for n = 2 to 41
M(n) = M(n-1)
for i = 0 to n-2
M(n) = M(n) + M(i)*M(n-2-i)
next i
print str$(M(n), "%19.0f"),
if isPrime(M(n)) then print "is a prime" else print : fi
next n
end
sub isPrime(v)
if v < 2 then return False : fi
if mod(v, 2) = 0 then return v = 2 : fi
if mod(v, 3) = 0 then return v = 3 : fi
d = 5
while d * d <= v
if mod(v, d) = 0 then return False else d = d + 2 : fi
wend
return True
end sub