数学研发网设为首页收藏本站
切换到窄版

 找回密码
 欢迎注册
数学研发论坛»论坛 › 【数学研究】 › 难题征解 › 阿弥陀签-ProjectEuler 837
123
返回列表 发新帖
楼主: iseemu2009

[讨论] 阿弥陀签-ProjectEuler 837

[复制链接]
wayne
发表于 2025-5-19 20:23:36 | 显示全部楼层
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)$
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

wayne
发表于 2025-5-19 21:46:27 | 显示全部楼层
C++/GMP计算9位数的超大模， f(987654321,123456789)(mod 10^200+357)  ，接近5分钟
  1. 59323489710030590937498974041811361729503527938950884150633376951881783733141136038416761027788735538124911444839486131816924092080909443328747109357429495905930956053696986596855219597555140472081399
复制代码


计算10位数的超大模，花了40分钟跑完
  1. ./nm 9876543211 1123456789 mod(10^200+357)
  2. 14504998039335138058545206947322707754384149282070629521936669837364927124521582648623340291701330931085574869168541809125058605840423114271914268430205327952153272578962424373053385826970817718933164
复制代码

相比之下，PARI/Gp慢了2-3倍，花了一个半小时跑完(1987654321,1234567891) mod(10^200+357)

  1. 74279609833465405311847193715038711033744229842155861970904036861058714593564477277402764734345492171616002302526713318977511160946438725161458628036851008656618220429643566021101278292868436244638766
复制代码


  1. HELLO(n,m)=
  2. {
  3.     local(M,a0,a1,a2,a3,a4,a5,a6,b0,b1,b2,b3,b4,b5,b6,b7,t);
  4.     a0=0;a1=0;a2=(1/48)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n);
  5.     a3=(1/768)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2;
  6.     a4=((1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2*(7 + n)^2)/30720;
  7.     a5=((1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)*(7 + n)*(9 + n)*(26355 + n*(12083 + n*(1749 + 85*n))))/185794560;
  8.     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)));
  9.     t=a6;
  10.     print([n,m,a0,a1,a2,a3,a4,a5,a6]);
  11.     for(k=0, floor((m-1)/2-6),
  12.         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);
  13.         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);
  14.         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);
  15.         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);
  16.         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);
  17.         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);
  18.         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);
  19.         b7=256*(3 + k)*(4 + k)*(5 + k)*(6 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k)*(13 + 2*k);
  20.         t=b0*a0+b1*a1+b2*a2+b3*a3+b4*a4+b5*a5+b6*a6;
  21.         t=-t/b7;
  22.         a0=a1;a1=a2;a2=a3;a3=a4;a4=a5;a5=a6;
  23.         a6=t;print([k,t])
  24.     );
  25. }

  27. WORLD(n,m)=
  28. {
  29.     my(walltime = getwalltime());
  30.     my(M,a0,a1,a2,a3,a4,a5,a6,b0,b1,b2,b3,b4,b5,b6,b7,t);
  31.     M=1234567891;
  32.     M=10^200+357;
  33.     a0=0;a1=0;a2=Mod((1/48)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n),M);
  34.     a3=Mod((1/768)*(1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2,M);
  35.     a4=Mod(((1 - (-1)^n)*(-1 + n)*(1 + n)*(3 + n)*(5 + n)^2*(7 + n)^2)/30720,M);
  36.     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);
  37.     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);
  38.     t=a6;

  40.     for(k=0, floor((m-1)/2-6),
  41.         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);
  42.         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);
  43.         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);
  44.         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 );
  45.         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);
  46.         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);
  47.         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);
  48.         b7=Mod(256*(3 + k)*(4 + k)*(5 + k)*(6 + k)*(7 + 2*k)*(9 + 2*k)*(11 + 2*k)*(13 + 2*k),M);
  49.         t=b0*a0+b1*a1+b2*a2+b3*a3+b4*a4+b5*a5+b6*a6;
  50.         t=-t/b7;
  51.         a0=a1;a1=a2;a2=a3;a3=a4;a4=a5;a5=a6;
  52.         a6=t;
  53.         if(k%1000000==0,printf("k=%d, time took %s\n",k,strtime(getwalltime() - walltime)))
  54.     );
  55.     component(t,2)
  56. }
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

wayne
发表于 2025-5-20 11:09:47 | 显示全部楼层
这道题,源头是Project Euler 837 ,https://projecteuler.net/problem=837
根据源头,我们得知这道题又叫 鬼脚图，鬼脚签、爬梯游戏，在日本称作Amidakuji, 阿弥陀签（あみだくじ），是一种游戏，也可以当作一种简易决策方法,用来随机生成序列的一个排列.

阿弥陀签，日本常见的抽签游戏你知道吗？
https://zhuanlan.zhihu.com/p/150081147

点评

wayne
我印象特别深刻的是大学期间我连续做了223,224,225三题,两道题是 挂机跑了通宵解决的  发表于 2025-5-21 07:27
mathe
project euler还是很不错的，已经允许很长时间了，论坛里面好像11年时访问的人比较多  发表于 2025-5-21 06:32
mathe
这个是非均匀分布，奇偶置换问题导致横杠数目可以确定哪一半状态可达。另外数目不是很大时，状态间分布应该很不均衡  发表于 2025-5-20 12:13
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

wayne
发表于 2025-5-20 20:40:13 | 显示全部楼层
挺有意思的,抛硬币是两个人的随机游戏, 要是两人以上的随机抽签,这个游戏值得付诸实践.
我想到几个问题

1) 给定 n个序列作为输入, 以及 指定的排序 作为输出. 如何生成最少且有效的横档方案
2) 给定 n个序列作为输入, 均匀随机的生成有效的横档的方案, 那么输出序列的 概率分布 如何
3) 实践中 会遇到哪些问题
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

wayne
发表于 2025-5-20 21:24:19 | 显示全部楼层
下载到了一篇论文, Statistical analysis on Amida-kuji , 说是概率分布函数 服从Fokker–Planck equation in one-dimensional diffusion process

Statistical analysis on Amida-kuji.pdf (239.53 KB, 下载次数: 3)

毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

123
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2025-6-3 06:40 , Processed in 0.054753 second(s), 21 queries .

Powered by Discuz! X3.5

© 2001-2025 Discuz! Team.

快速回复 返回顶部 返回列表