wayne 发表于 2025-5-19 00:54
搞到直接关于结果的递推公式了,是一个7阶差分方程,速度奇快. 前面的 计算f(12345,54321)需要3分钟,而这个只 ...
给定n，该差分方程可以通过反复差分运算，降低系数里的k的次数，进而完全消元，最终变成关于k为常系数的n+1阶差分方程。于是继而可以利用矩阵的快速模幂运算，让时间复杂度变成$O(lgn)$，空间复杂度是$O(n^2)$，
最终是一个很简洁的递推公式 $a(n+1+k)=\sum _{i=1}^{n+1} (-1)^{i-1} C_{n+1}^{i} a(n+1+k-i)$
C++/GMP计算9位数的超大模， f(987654321,123456789)(mod 10^200+357)，接近5分钟
59323489710030590937498974041811361729503527938950884150633376951881783733141136038416761027788735538124911444839486131816924092080909443328747109357429495905930956053696986596855219597555140472081399
计算10位数的超大模，花了40分钟跑完
./nm 9876543211 1123456789 mod(10^200+357)
14504998039335138058545206947322707754384149282070629521936669837364927124521582648623340291701330931085574869168541809125058605840423114271914268430205327952153272578962424373053385826970817718933164
相比之下，PARI/Gp慢了2-3倍，花了一个半小时跑完(1987654321,1234567891) mod(10^200+357)
74279609833465405311847193715038711033744229842155861970904036861058714593564477277402764734345492171616002302526713318977511160946438725161458628036851008656618220429643566021101278292868436244638766
HELLO(n,m)=
{
local(M,a0,a1,a2,a3,a4,a5,a6,b0,b1,b2,b3,b4,b5,b6,b7,t);
a0=0;a1=0;a2=(1/48)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n);
a3=(1/768)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2;
a4=((1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2*(7 + n)^2)/30720;
a5=((1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)*(7 + n)*(9 + n)*(26355 + n*(12083 + n*(1749 + 85*n))))/185794560;
a6=(1/7431782400)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)*(7 + n)*(9 + n)*(11 + n)^2*(9345 + n*(4361 + n*(615 + 31*n)));
t=a6;
print();
for(k=0, floor((m-1)/2-6),
b0=-9*(2*k + n)*(1 + 2*k + n)*(2 + 2*k + n)*(3 + 2*k + n)*(5 + 2*k + n)*(7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n);
b1=3*(2 + 2*k + n)*(3 + 2*k + n)*(5 + 2*k + n)*(7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n)*(252 + 286*k + 92*k^2 + 65*n + 44*k*n + 5*n^2);
b2=(-(5 + 2*k + n))*(7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n)*(77112 + 138600*k + 95288*k^2 + 29856*k^3 + 3616*k^4 + 32958*n + 44624*k*n + 20880*k^2*n + 3392*k^3*n + 4565*n^2 + 4392*k*n^2 + 1104*k^2*n^2 + 294*n^3 + 152*k*n^3 + 7*n^4);
b3=(7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n)*(2117304 + 3730728*k + 2661864*k^2 + 962192*k^3 + 176352*k^4 + 13120*k^5 + 620706*n + 894164*k*n + 490952*k^2*n + 121888*k^3*n + 11552*k^4*n + 66973*n^2 + 75846*k*n^2 + 29088*k^2*n^2 + 3776*k^3*n^2 + 3653*n^3 + 2844*k*n^3 + 560*k^2*n^3 + 83*n^4 + 34*k*n^4 + n^5);
b4=-4*(3 + k)*(7 + 2*k)*(9 + 2*k + n)*(11 + 2*k + n)*(449574 + 526583*k + 234142*k^2 + 46852*k^3 + 3560*k^4 + 95911*n + 87604*k*n + 27192*k^2*n + 2872*k^3*n + 8101*n^2 + 5253*k*n^2 + 870*k^2*n^2 + 329*n^3 + 112*k*n^3 + 5*n^4);
b5=4*(3 + k)*(4 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k + n)*(147123 + 109870*k + 27540*k^2 + 2312*k^3 + 23087*n + 12252*k*n + 1660*k^2*n + 1437*n^2 + 414*k*n^2 + 33*n^3);
b6=-64*(3 + k)*(4 + k)*(5 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k)*(1203 + 504*k + 52*k^2 + 136*n + 32*k*n + 5*n^2);
b7=256*(3 + k)*(4 + k)*(5 + k)*(6 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k)*(13 + 2*k);
t=b0*a0+b1*a1+b2*a2+b3*a3+b4*a4+b5*a5+b6*a6;
t=-t/b7;
a0=a1;a1=a2;a2=a3;a3=a4;a4=a5;a5=a6;
a6=t;print()
);
}
WORLD(n,m)=
{
my(walltime = getwalltime());
my(M,a0,a1,a2,a3,a4,a5,a6,b0,b1,b2,b3,b4,b5,b6,b7,t);
M=1234567891;
M=10^200+357;
a0=0;a1=0;a2=Mod((1/48)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n),M);
a3=Mod((1/768)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2,M);
a4=Mod(((1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2*(7 + n)^2)/30720,M);
a5=Mod(((1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)*(7 + n)*(9 + n)*(26355 + n*(12083 + n*(1749 + 85*n))))/185794560,M);
a6=Mod((1/7431782400)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)*(7 + n)*(9 + n)*(11 + n)^2*(9345 + n*(4361 + n*(615 + 31*n))),M);
t=a6;
for(k=0, floor((m-1)/2-6),
b0=Mod(-9*(2*k + n)*(1 + 2*k + n)*(2 + 2*k + n)*(3 + 2*k + n)*(5 + 2*k + n)*(7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n),M);
b1=Mod(3*(2 + 2*k + n)*(3 + 2*k + n)*(5 + 2*k + n)*(7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n)*(252 + 286*k + 92*k^2 + 65*n + 44*k*n + 5*n^2),M);
b2=Mod((-(5 + 2*k + n))*(7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n)*(77112 + 138600*k + 95288*k^2 + 29856*k^3 + 3616*k^4 + 32958*n + 44624*k*n + 20880*k^2*n + 3392*k^3*n + 4565*n^2 + 4392*k*n^2 + 1104*k^2*n^2 + 294*n^3 + 152*k*n^3 + 7*n^4),M);
b3=Mod((7 + 2*k + n)*(9 + 2*k + n)*(11 + 2*k + n)*(2117304 + 3730728*k + 2661864*k^2 + 962192*k^3 + 176352*k^4 + 13120*k^5 + 620706*n + 894164*k*n + 490952*k^2*n + 121888*k^3*n + 11552*k^4*n + 66973*n^2 + 75846*k*n^2 + 29088*k^2*n^2 + 3776*k^3*n^2 + 3653*n^3 + 2844*k*n^3 + 560*k^2*n^3 + 83*n^4 + 34*k*n^4 + n^5),M );
b4=Mod(-4*(3 + k)*(7 + 2*k)*(9 + 2*k + n)*(11 + 2*k + n)*(449574 + 526583*k + 234142*k^2 + 46852*k^3 + 3560*k^4 + 95911*n + 87604*k*n + 27192*k^2*n + 2872*k^3*n + 8101*n^2 + 5253*k*n^2 + 870*k^2*n^2 + 329*n^3 + 112*k*n^3 + 5*n^4), M);
b5=Mod(4*(3 + k)*(4 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k + n)*(147123 + 109870*k + 27540*k^2 + 2312*k^3 + 23087*n + 12252*k*n + 1660*k^2*n + 1437*n^2 + 414*k*n^2 + 33*n^3),M);
b6=Mod(-64*(3 + k)*(4 + k)*(5 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k)*(1203 + 504*k + 52*k^2 + 136*n + 32*k*n + 5*n^2),M);
b7=Mod(256*(3 + k)*(4 + k)*(5 + k)*(6 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k)*(13 + 2*k),M);
t=b0*a0+b1*a1+b2*a2+b3*a3+b4*a4+b5*a5+b6*a6;
t=-t/b7;
a0=a1;a1=a2;a2=a3;a3=a4;a4=a5;a5=a6;
a6=t;
if(k%1000000==0,printf("k=%d, time took %s\n",k,strtime(getwalltime() - walltime)))
);
component(t,2)
}
这道题,源头是Project Euler 837 ,https://projecteuler.net/problem=837
根据源头,我们得知这道题又叫 鬼脚图，鬼脚签、爬梯游戏，在日本称作Amidakuji, 阿弥陀签（あみだくじ），是一种游戏，也可以当作一种简易决策方法,用来随机生成序列的一个排列.
阿弥陀签，日本常见的抽签游戏你知道吗？
https://zhuanlan.zhihu.com/p/150081147
挺有意思的,抛硬币是两个人的随机游戏, 要是两人以上的随机抽签,这个游戏值得付诸实践.
我想到几个问题
1) 给定 n个序列作为输入, 以及 指定的排序 作为输出. 如何生成最少且有效的横档方案
2) 给定 n个序列作为输入, 均匀随机的生成有效的横档的方案, 那么输出序列的 概率分布 如何
3) 实践中 会遇到哪些问题
下载到了一篇论文, Statistical analysis on Amida-kuji , 说是概率分布函数 服从Fokker–Planck equation in one-dimensional diffusion process