字数太多,上图 Mahjong

麻将的游戏是玩s种花色的牌.每张牌也有一个在范围1..n内的数字,对于每个花色/数字组合,正好有四张相同花色和数字的牌。(真正的麻将游戏也包含其他奖赏牌,但这些在此问题中不考虑。)

糊牌的一手是3t+2张牌(其中t是固定整数),它可以安排为t个 三重和一对,其中:

  • 三重是顺子或碰
  • 顺子是三张相同花色和连续数字的牌
  • 碰是三张相同的牌(相同花色和相同数字)
  • 对是两张相同的牌(相同花色和相同数字)

例如,下面是一只胜利的手,它包括两个顺子,两个碰和一对:

请注意,有时相同的切片集合可以以多种方式表示为三重和一对。这只能算作一只糊牌的手。例如,下面被认为是与上面相同的糊牌手,因为它们由相同的牌组成:

设 w(n,s,t)是t个三重和一对形成不同的糊牌手的数量,其中牌的花色数为s,数字最大为n。

例如,对于单个花色和牌数字最大为 4,我们有w(4,1,1)=20 :其中有12只糊牌的手,包括一个碰和一对,另外8个糊牌的手,包括一只顺子和一对。已知

w(9,1,4)=13259 w(9,3,4)=5237550 w(1000,1000,5)≡107662178(mod 1000000007)

求:w(10^8 ,10^8 ,30)给出mod 1000000007的答案.

用sql模拟w(4,1,1)

with p as
(select level||','||level||',' n, power(10,level)*2 flg from dual connect by level<=9)
,c as
(select level||','||level||','||level||',' n , power(10,level)*3 flg from dual connect by level<=9
union
select level||','||(level+1)||','||(level+2)||',' n, power(10,level)*111 flg from dual connect by level<=9-2)
,t(lv,s,flg) as(
select 0,n s,flg from p
union all
select lv+1,s||n,t.flg+c.flg from t,c where lv<4 and regexp_instr(t.flg+c.flg,'[5-9]')=0)
select count(*),count(distinct flg) from t where lv=4 ;

COUNT(*) COUNT(DISTINCTFLG)
-------- ------------------
  300366              13259

用相同思路实现的python递归程序如下:

import re
def pai(n):
 c=[[i,i] for i in range(1,n+1)]
 return c

def peng(n):
 c=[[i,i,i] for i in range(1,n+1)]
 return c

def chow(n):
 c=[[i,i+1,i+2] for i in range(1,n-1)]
 return c


x=chow(9)+peng(9)
x.sort()
lenx=len(x)
xflg=[0]*lenx
for i in range(0,lenx):
 for j in x[i]:
  xflg[i]+=10**j


print(xflg) 
cnt=0

def mj(lv,flg,last_i):
 global cnt
 global set_a
 if last_i==lenx: return -1
 if lv==4 : cnt+=1;set_a.add(flg);return 0 
 for i in range(last_i,lenx):
  newflg= xflg[i]+flg
  m=re.search('[5-9]',str(newflg))
  if m==None: #each number occurs less than 5 times
   mj(lv+1,newflg,i)

set_a=set()
cnt=0
for i in pai(9):
 mj(0,(10**i[0])*2,0)


print(cnt) 
print(len(set_a))

def wt():
 fileObject = open('sampleList.txt', 'w')
 for ip in set_a:
     fileObject.write(str(ip))
     fileObject.write('\n')
 fileObject.close()
输出结果
>>> print(cnt)
14738
>>> print(len(set_a))
13259

修改上述代码中的参数,得到w(9,3,4)

x=chow(29)+peng(29)
x=[i for i in x if 10 not in i and 20 not in i]
x.sort()

set_a=set()
cnt=0
pa=[i for i in pai(29) if 10 not in i and 20 not in i]
for i in pa:
 mj(0,(10**i[0])*2,0)

>>>> print(cnt)
5296275
>>>> print(len(set_a))
5237550 

重写了mj2函数和计算总数的w函数

def mj2(lv,flg,last_i,s,t):
 global cnt,lenx
 #print(lenx)
 if last_i==lenx: return -1
 if lv==t : cnt+=1;set_a.add(flg);return 0 
 for i in range(last_i,lenx):
  newflg= xflg[i]+flg
  m=re.search('[5-9]',str(newflg))
  if m==None: 
   mj2(lv+1,newflg,i,s,t)

def w(n,s,t): 
 global x,xflg,pa,lenx,cnt
 x=[]
 for i in chow((n+1)*s-1)+peng((n+1)*s-1):
  for j in range(1,s+1):
   if  (n+1)*j  in i:break 
  else: x.append(i)
 x.sort()
 lenx=len(x)
 xflg=[0]*lenx
 for i in range(0,lenx):
  for j in x[i]:
   xflg[i]+=10**j
 cnt=0
 pa=[]
 for i in pai((n+1)*s-1):
  for j in range(1,s+1):
   if  (n+1)*j  in i:break 
  else: pa.append(i)
 for i in pa:
  mj2(0,(10**i[0])*2,0,s,t)
 return len(set_a)  

for i in range(0,5):set_a=set();x=[];xflag=[];pa=[];lenx=0;cnt=0;print(i+4,w(4+i,1,1))  
4 20
5 35
6 54
7 77
8 104

for i in range(0,5):set_a=set();x=[];xflag=[];pa=[];lenx=0;cnt=0;print(i+4,w(4+i,4,1))
4 368
5 620
6 936
7 1316
8 1760 

总结出w(n,s,1)=ns(2ns-2s-1)

测试w(n,s,2),计算差分,应该是3次方程,用4组值,mathematica解方程组得出w(n,s,2)=s(2s^2n^3-2s(2s+1)n^2+(s+1)(2s-1)n+6)

for i in range(0,20):set_a=set();x=[];xflag=[];pa=[];lenx=0;cnt=0;print(i+4,w(4+i,5,2))
4 8310
5 18880
6 35850
7 60720
8 94990
In[17]:= Solve[
     a*4^3+b*4^2+c*4+d==8310  &&
     a*5^3+b*5^2+c*5+d==18880 &&
     a*6^3+b*6^2+c*6+d==35850 &&
     a*7^3+b*7^2+c*7+d==60720,
     {a,b,c,d}]

Out[17]= {{a -> 250, b -> -550, c -> 270, d -> 30}

t=3的公式

In[115]:= Factor[w53[x]]
                                2         3        4
          5 (-192 - 6 x + 1525 x  - 1650 x  + 500 x )
Out[115]= -------------------------------------------
                               3
In[104]:= Factor[w43[x]]
                                 2        3        4
          4 (-156 + 141 x + 764 x  - 864 x  + 256 x )
Out[104]= -------------------------------------------
                               3
In[93]:= Factor[w33[x]]
                              2        3       4
Out[93]= 3 (-40 + 64 x + 101 x  - 126 x  + 36 x )
In[94]:= Factor[w23[x]]
                              2        3       4
         2 (-84 + 171 x + 70 x  - 120 x  + 32 x )
Out[94]= ----------------------------------------
                            3
In[95]:= Factor[w13[x]]
                          2       3      4
         -48 + 102 x - 7 x  - 18 x  + 4 x
Out[95]= ---------------------------------
                         3

容易观察到x^4的系数是4 s^3, 常数项的系数是-36 s -12,所以只要求x、x^2、x^3的系数,设x^3的系数是 a s^3+b s^2 +c s +d, 利用上述公式的s=1~4项联立方程组:

Solve[
a 1^3+b 1^2 +c 1 +d ==-18  &&
a 2^3+b 2^2 +c 2 +d ==-120 &&
a 3^3+b 3^2 +c 3 +d ==-126*3 &&
a 4^3+b 4^2 +c 4 +d ==-864  ,
{a,b,c,d}]

解得{{a -> -12, b -> -6, c -> 0, d -> 0}},所以系数是-12 s^2 -6 s,用s=5代入值为-1650,与w53的表达式比较,结果正确。同样方法求得x^2的系数12 s^3 +6 s^2 -25 s,x的系数-4 s^3 +97 s +9,所以得到wns3的通项公式为wns3[n_,s_]:=(-36 s -12 +(-4 s^3 +97 s +9)n +(12 s^3 +6 s^2 -25 s)n^2+ (-12 s^3-6 s^2 )n^3+ 4 s^3 n^4)s /3

当t=4时,和t=3的情况类似,前几个数值无法用公式描述,要从w(9,1,4)开始符合公式

wn14[n_]:=2/3 n^5-4 n^4-26/3 n^3+96 n^2-140 n-61

利用相同办法,求得:

wn24[n_]:=64/3 n^5-320/3 n^4+128/3 n^3+1712/3 n^2-822 n-46 

在笔记本电脑上运行程序算w(10,3,4)时,产生内存不足错误。考虑到列表的存储空间是集合的1/4(参考https://www.ituring.com.cn/article/509014),并且重复数占总数的比例很低,而计算不重复数只需要最后一步算,因此改写mj2函数为 而在w函数最后加上排序和不重复计数代码。

def mj4(lv,flg,last_i,s,t):
 global cnt,lenx
 #print(lenx)
 if last_i==lenx: return -1
 if lv==t : 
  cnt+=1
  sortlst_a.append(flg)
  return 0
 for i in range(last_i,lenx):
  newflg= xflg[i]+flg
  m=re.search('[5-9]',str(newflg))
  if m==None:
   mj4(lv+1,newflg,i,s,t)


def w3(n,s,t):
 global x,xflg,pa,lenx,cnt
 x=[]
 for i in chow((n+1)*s-1)+peng((n+1)*s-1):
  for j in range(1,s+1):
   if  (n+1)*j  in i:break
  else: x.append(i)
 x.sort()
 lenx=len(x)
 xflg=[0]*lenx
 for i in range(0,lenx):
  for j in x[i]:
   xflg[i]+=10**j
 cnt=0
 pa=[]
 for i in pai((n+1)*s-1):
  for j in range(1,s+1):
   if  (n+1)*j  in i:break
  else: pa.append(i)
 for i in pa:
  mj4(0,(10**i[0])*2,0,s,t)
 slen= len(sortlst_a)
 sortlst_a.sort()
 ucnt=0
 oldv=0
 for i in range(0,slen):
  if sortlst_a[i] !=oldv: ucnt+=1;oldv=sortlst_a[i]
 return ucnt

这样改造后,利用程序最多可以算到w(12,3,4),离计算5次多项式方程组所需的等式还差2个。
好在根据wnst在t=1,2时的结果,通项公式的最高次系数有规律,2/3,64/3,总结为2*s^t/3,所以最高项次数为2*3^5/3=162(实际上n、s的次数应该相同)。这样减少一个未知数,只需要额外的1个等式。
再次根据上述经验,wn34的实际值在n=9前比公式值略小,小的这个值一定是3的倍数,所以尝试在w(8,3,4)的实际值上加3,作为额外的等式,得到如下方程组:

In[15]:= Solve[
         162*8 ^5+x*8 ^4+a*8 ^3+b*8 ^2+c*8 +d==2661693 &&
         162*9 ^5+x*9 ^4+a*9 ^3+b*9 ^2+c*9 +d==5237550 &&
         162*10^5+x*10^4+a*10^3+b*10^2+c*10+d==9496287 &&
         162*11^5+x*11^4+a*11^3+b*11^2+c*11+d==16149624&&
         162*12^5+x*12^4+a*12^3+b*12^2+c*12+d==26085537,
         {x,a,b,c,d}]

Out[15]= {{x -> -756, a -> 738, b -> 1416, c -> -2343, d -> 117}}

得出

wn34[n_]:=162 n^5-756 n^4+738 n^3+1416 n^2-2343 n+117