Motzkin numbers
- Definition
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).
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
Factor
<lang 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</lang>
- 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
Julia
<lang 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
</lang>
- 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
Perl
<lang 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' : ;
}</lang>
- 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
Alternative. Note that M[37] happens by chance to be correct under 32-bit and p2js, even though an intermediate value exceeds precision (in this case being rounded to the nearest whole multiple of 16), and the final divide by (i+1) results in a whole integer simply because there isn't enough precision to hold any decimal places. Output as above on 64-bit, less four entries under 32-bit and pwa/p2js.
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
Raku
Using binomial coefficients and Catalan numbers
<lang perl6>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 };</lang>
- 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
<lang perl6>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 };</lang>
Same output
REXX
<lang 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 the Motzkin number ──► term. */
end /*m*/
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 /*is X ÷ by 5? Then return not prime.*/
if x//2 ==0 then return 0 /* " " " " 2? " " " " */
if x//3 ==0 then return 0 /* " " " " 3? " " " " */
if x//7 ==0 then return 0 /* " " " " 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.*/</lang>
- 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 ─── ─────────────────────── ───────────
Rust
<lang 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)
);
}
}</lang>
- 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
Wren
<lang ecmascript>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)
}</lang>
- 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