Catalan 数之Python演示

Catalan 数之Python演示

关于Catalan 数，英文的下面网页可以参考：

https://brilliant.org/wiki/catalan-numbers/

以下卡特兰数（Catalan number）简介来自下面网页，Python演示代码为原创。

http://blog.miskcoo.com/2015/07/catalan-number

卡特兰数（Catalan number）是组合数学中一个重要的计数数列，在很多看似毫不相关地方都能见到它的身影

带限制条件的路径总数

首先我们来看一个问题：

在一个平面直角坐标系中，只能往右或往上走一个单位长度，问有多少种不同的路径可以从左下角 $(1,1)$ 走到右上角 $(n,n)$，并且要求路径不能经过直线 $y=x$ 上方的点，下图中的路径都是合法的（图片来源 Wikipedia）

如果没有限制条件，那么从左下角走到右上角一共有 $2n$ 步，有 $n$ 步是往右，另外 $n$ 步是往上，那么路径方案数就是 $2n$ 步中选择 $n$步往右，一共有 $\left(\genfrac{}{}{0ex}{}{2n}{n}\right)$ （即 $C_{2n}^n$）种方案。

那么我们考虑一下这里面有多少种方案是不合法的。

首先对于每一种不合法的方案，它的路径一定与 $y=x+1$ 有交。我们找到它与 $y=x+1$ 的第一个交点，然后将这个点后面部分的路径关于$y=x+1$ 做一个对称。由于原来路径到达 $(n,n)$，新的对称之后的路径就会到达 $(n-1,n+1)$。这样我们把每一种不合法方案都对应到了一条从 $(1,1)$ 到 $(n-1,n+1)$的路径。

现在再来看是否每一条这样的路径都能对应到一种不合法方案，如果是，那么这就建立了一个一一映射的关系，也就是它们的方案总数相同。这是肯定的，因为每一条这样的路径必定与$y=x+1$ 有交，那么对称回去，就得到一条不合法方案。

由于从 $(1,1)$ 到 $(n-1,n+1)$ 的路径有 $\left(\genfrac{}{}{0ex}{}{n-1+n+1}{n-1}\right)$
条，那么合法的方案就是

$$C_n=\left(\genfrac{}{}{0ex}{}{2n}{n}\right)-\left(\genfrac{}{}{0ex}{}{2n}{n-1}\right)=\frac{1}{n+1}\left(\genfrac{}{}{0ex}{}{2n}{n}\right)$$

我们把这个方案数记为 $C_n$，这就是著名的 Catalan 数

我们来看看它的前几项（$n$ 从 $0$ 开始）

1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012, 742900,
2674440, 9694845, 35357670, 129644790, 477638700, 1767263190

Python演示代码说明

合法路径要求不越过对角线的本质意义，就是要求合法路径的每一步的坐标$x\ge y$，也即向右的步数必须大于或等于向上的步数。类比到计算机动态内存管理中，就是进栈动作数量必须大于或等于出栈动作数量，这是必然的。在括号配对问题中，则是左括号必须必须大于或等于右括号的数量，lisp语言中使用一系列的括号来区分代码的层次结构，非常能体现Catalan 数的意义。语言结构层次可以通过括号来区分，也可以通过二叉树来展示，括号配对与二叉树是可以互相映射的。

为了简单，代码演示中，选择了最简单的“田”字格，即$n=2$，从(0,0)出发，只能向左或向上，一共有 $\left(\genfrac{}{}{0ex}{}{2\ast 2}{2}\right)=6$种路径方案到达(2,2)，但仅有2种方案未越过对角线$x=y$，属于合法路径。这不难理解。

Catalan 数理解中有难度的地方，就是如何计算非法路径的数量。如上面转引文字所述，就是把每一条非法路径的越线部分，进行翻转，使终点到达$(n-1,n+1)$，就是走“目”字格。目字格，显然只有4条路径，我们可以通过的横线从下到上编号$1、2、3、4$。数学家们巧妙的地方就是发现目字格的4条路径与田字格的4条非法路径是一一对应的。

演示图中，天蓝色的斜线是翻转对称轴$y=x+1$。普通路径画为绿色，需要翻转的部分画为黄色，对称翻转后画为紫色。

第1条、第2条是合法路径，
$$[[0,0],[1,0],[2,0],[2,1],[2,2]]$$
$$[[0,0],[1,0],[1,1],[2,1],[2,2]]$$
每一步都是横坐标大于纵坐标，没有越过对角线。

第3条是非法路径,
$$[[0,0],[1,0],[1,1],[1,2],[2,2]]$$
因为第3步，就是第4个坐标点，横坐标小于纵坐标，$1<2$。从该点以后的路径需要以$y=x+1$为对称轴，进行翻转，以映射到目字格第1条路径。

第4条是非法路径,
$$[[0,0],[0,1],[1,1],[2,1],[2,2]]$$
因为第1步，就是第2个坐标点，横坐标小于纵坐标，$0<1$。从该点以后的路径需要以$y=x+1$为对称轴，进行翻转，以映射到目字格第4条路径。

第5条是非法路径,
$$[[0,0],[0,1],[1,1],[1,2],[2,2]]$$
因为第1步，就是第2个坐标点，横坐标小于纵坐标，$0<1$。从该点以后的路径需要以$y=x+1$为对称轴，进行翻转，以映射到目字格第3条路径。

第6条是非法路径,
$$[[0,0],[0,1],[0,2],[1,2],[2,2]]$$
因为第1步，就是第2个坐标点，横坐标小于纵坐标，$0<1$。从该点以后的路径需要以$y=x+1$为对称轴，进行翻转，以映射到目字格第2条路径。

求一点关于直线的对称点公式见下面网页，

https://tieba.baidu.com/p/5198944956

若点$A_1(x_1,y_1)$关于直线$y=kx+b$成轴对称，则对应点$A_2(x_2,y_2)$的坐标为：

$$x_2=\frac{2\cdot k\cdot \left(y_1-b\right)-x_1\cdot \left(k^2-1\right)}{k^2+1}$$

$$y_2=\frac{2\cdot k\cdot \left(b+x_{1}k\right)+y_1\cdot \left(k^2-1\right)}{k^2+1}$$

由于本问题中对称轴$y=x+1$特殊，实际计算对称点的代码很简单：
endc0=(end[1]-1)
endc1=(1+end[0])

Python代码

# coding=utf-8

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import MultipleLocator 

def plt_arrow(x_begin,y_begin,x_end,y_end,color,plt):#画带箭头的线段

    plt.arrow(x_begin, y_begin, x_end - x_begin, y_end - y_begin,
             length_includes_head=True,     # 增加的长度包含箭头部分
             head_width = 0.15, head_length =0.15, fc=color, ec=color)
    
def route1(path, plt,color):#处理一条路径，就是画一个子图
    print(path)
    lpath =len(path)
    isCantelan = True
    beginc0=0#保存越过限制线后箭头的起点坐标
    beginc1=0
    for i in range(lpath-1):
        begin=path[i]
        end=path[i+1]
        #print(begin,end)
        if(isCantelan):
            if(begin[1]>begin[0]):#该点以后越过限制线：y=x，换颜色
                isCantelan = False
                plt_arrow(begin[0], begin[1], end[0], end[1], "orange",plt)
                
                #求end点关于y=x+1的对称点
                #print(end)
                endc0=(end[1]-1)
                endc1=(1+end[0])
                #print("映射对称点坐标：",endc0,endc1)
                plt_arrow(begin[0], begin[1], endc0, endc1, "magenta",plt)#画映射后的路径
                #保存对称点坐标，为下一个箭头的起点
                beginc0=endc0
                beginc1=endc1
            else:
                plt_arrow(begin[0], begin[1], end[0], end[1], color,plt)#未越过限制线的路径
            
        else:
            plt_arrow(begin[0], begin[1], end[0], end[1], "orange",plt)
            
            #求end点关于y=x+1的对称点
            #print(end)
            endc0=(end[1]-1)
            endc1=(1+end[0])
            #print("映射对称点坐标：",endc0,endc1)
            plt_arrow(beginc0, beginc1, endc0, endc1, "magenta",plt)#画映射后的路径
            #保存对称点坐标，为下一个箭头的起点
            beginc0=endc0
            beginc1=endc1
            
    return isCantelan #返回该路径是否属于卡特兰路径


def buildpath(n):    
    #构造n为参数的全部路径
    allpath=[]#全部路径
    p =[]      
    p.append([0,0])
    needstep=n

    allpath.append(p)
    for i in range(4):
        newallpath=[]#新全部路径
        for j,path in enumerate(allpath):
            #print("\n当前路径：",j+1,path)

            endpoint=path[-1][:]#取已有最后坐标作为下一步的基础 
            pxstep=endpoint[0]
            pystep=endpoint[1]

            if pxstep<needstep:    
                newpointx=endpoint[:]
                newpointx[0]=newpointx[0]+1
                pathx=path[:]           
                pathx.append(newpointx)
                #print("扩展后路径：",pathx)
                newallpath.append(pathx)

            if pystep<needstep: 
                newpointy=endpoint[:]
                newpointy[1]=newpointy[1]+1
                pathy=path[:]           
                pathy.append(newpointy)
                #print("扩展后路径：",pathy)
                newallpath.append(pathy)

        #print("\n第",i+1,"步所有路径",len(newallpath),"条：")
        #[print(x) for x in newallpath]
        allpath=newallpath[:]
    return allpath

def drawallpath(allpath):

    fig,axes = plt.subplots(2,3,figsize=(4.4,4),sharex=True,sharey=True)
    plt.subplots_adjust(wspace=0.2,hspace=0.2)
    xmajorLocator  = MultipleLocator(1) #将x主刻度标签设置为1的倍数
    ymajorLocator  = MultipleLocator(1) #将y轴主刻度标签设置为1的倍数 

    for i in range(2):
        for j in range(3):
            id=i*3+j
            #print(id)
            axes[i,j].set_xlim(-0.1, 2.1)#指定画图范围
            axes[i,j].set_ylim(-0.1, 3.1)
            axes[i,j].xaxis.set_major_locator(xmajorLocator) 
            axes[i,j].yaxis.set_major_locator(ymajorLocator) 
            axes[i,j].xaxis.grid(True, which='major') #x坐标轴的网格使用主刻度 
            axes[i,j].yaxis.grid(True, which='major') #y坐标轴的网格使用次刻度 
            
            axes[i,j].grid(ls=":",c='gray',)#打开坐标网格
            
            axes[i,j].spines['right'].set_visible(False)#去掉边框
            axes[i,j].spines['top'].set_visible(False)
            axes[i,j].spines['left'].set_visible(False)
            axes[i,j].spines['bottom'].set_visible(False)
            
            axes[i,j].scatter(2, 2, s = 30, color = 'r',marker='o')#目标终点
            axes[i,j].scatter(1, 3, s = 30, color = 'purple',marker='o')#映射后的终点
            
            lc=[0,3]#限制线坐标
            axes[i,j].plot(lc,lc,color='blue', linewidth=0.7, linestyle='dashed') # 限制线
            
            scx=[0,2]# 翻转对称线x坐标
            scy=[1,3]# 翻转对称线y坐标
            axes[i,j].plot(scx,scy,color='cyan', linewidth=0.7, linestyle='dashed') # 翻转对称线
            
            isCantelan=route1(allpath[id],axes[i,j],"green")  
                       
            if isCantelan:
                axes[i,j].set_title(str(id+1)+' Cantelan')
            else:
                axes[i,j].set_title(str(id+1)+' Not')
                
        plt.savefig(r"data\Catalan.png",dpi=72) #保存分辨率为72的图片

allpath=buildpath(2)           
drawallpath(allpath)