阿弥陀签-ProjectEuler 837

wayne

mathe的计算结果 可以换成更加简洁的四阶差分方程, $a(n+4)=x a(n+3)+2 a(n+2)+(x^3-x) a(n+1)+(-x^4+2 x^2-1) a(n)$, 其中初始项是$a(1)=0,a(2)=x^2+1,a(3)=2 x^3,a(4)=x^4+4 x^2+1$

1. {1,0}
   
2. {2,1+x^2}
   
3. {3,2 x^3}
   
4. {4,1+4 x^2+x^4}
   
5. {5,2 x^3 (4+x^2)}
   
6. {6,1+9 x^2+9 x^4+3 x^6}
   
7. {7,2 x^3 (10+10 x^2+x^4)}
   
8. {8,1+16 x^2+36 x^4+30 x^6+3 x^8}
   
9. {9,2 x^3 (4+x^2) (5+10 x^2+2 x^4)}
   
10. {10,1+25 x^2+100 x^4+156 x^6+57 x^8+3 x^10}
    
11. {11,2 x^3 (35+140 x^2+126 x^4+38 x^6+2 x^8)}
    
12. {12,1+36 x^2+225 x^4+568 x^6+441 x^8+90 x^10+5 x^12}
    
13. {13,2 x^3 (4+x^2) (14+84 x^2+119 x^4+54 x^6+2 x^8)}
    
14. {14,1+49 x^2+441 x^4+1645 x^6+2185 x^8+981 x^10+155 x^12+5 x^14}
    
15. {15,2 x^3 (84+756 x^2+1890 x^4+1900 x^6+738 x^8+90 x^10+3 x^12)}
    
16. {16,1+64 x^2+784 x^4+4060 x^6+8200 x^8+6436 x^10+2066 x^12+230 x^14+5 x^16}
    
17. {17,2 x^3 (4+x^2) (30+360 x^2+1233 x^4+1708 x^6+909 x^8+126 x^10+3 x^12)}
    
18. {18,1+81 x^2+1296 x^4+8904 x^6+25380 x^8+30726 x^10+16766 x^12+3906 x^14+315 x^16+7 x^18}
    
19. {19,2 x^3 (165+2640 x^2+12936 x^4+27874 x^6+28105 x^8+12880 x^10+2586 x^12+192 x^14+3 x^16)}
    
20. {20,1+100 x^2+2025 x^4+17832 x^6+68124 x^8+117558 x^10+97790 x^12+38900 x^14+6741 x^16+448 x^18+7 x^20}
    
21. {21,2 x^3 (4+x^2) (55+1100 x^2+6853 x^4+18942 x^6+24717 x^8+14532 x^10+3466 x^12+236 x^14+4 x^16)}
    
22. {22,1+121 x^2+3025 x^4+33231 x^6+163812 x^8+381612 x^10+450044 x^12+272840 x^14+81593 x^16+11221 x^18+595 x^20+7 x^22}
    
23. {23,2 x^3 (286+7150 x^2+57915 x^4+217360 x^6+416416 x^8+416052 x^10+217726 x^12+57796 x^14+7056 x^16+340 x^18+4 x^20)}
    
24. {24,1+144 x^2+4356 x^4+58410 x^6+360855 x^8+1090584 x^10+1731744 x^12+1486836 x^14+680805 x^16+160140 x^18+17766 x^20+756 x^22+9 x^24}
    
25. {25,2 x^3 (4+x^2) (91+2730 x^2+26845 x^4+123240 x^6+291226 x^8+360672 x^10+231571 x^12+72396 x^14+9286 x^16+420 x^18+4 x^20)}
    
26. {26,1+169 x^2+6084 x^4+97812 x^6+740025 x^8+2814669 x^10+5791032 x^12+6687864 x^14+4346217 x^16+1561575 x^18+296346 x^20+26838 x^22+981 x^24+9 x^26}
    
27. {27,2 x^3 (455+16380 x^2+198198 x^4+1149720 x^6+3579147 x^8+6242184 x^10+6238260 x^12+3580038 x^14+1151010 x^16+197268 x^18+16416 x^20+540 x^22+5 x^24)}
    
28. {28,1+196 x^2+8281 x^4+157248 x^6+1429857 x^8+6682104 x^10+17292249 x^12+25854624 x^14+22623993 x^16+11529462 x^18+3339651 x^20+519996 x^22+39591 x^24+1224 x^26+9 x^28}
    
29. {29,2 x^3 (4+x^2) (140+5880 x^2+83706 x^4+574740 x^6+2132160 x^8+4458180 x^10+5347491 x^12+3653562 x^14+1362831 x^16+255276 x^18+21066 x^20+660 x^22+5 x^24)}
    
30. {30,1+225 x^2+11025 x^4+244153 x^6+2627001 x^8+14794065 x^10+47030685 x^12+88392465 x^14+100190601 x^16+68741841 x^18+28227717 x^20+6719625 x^22+876315 x^24+56727 x^26+1485 x^28+11 x^30}
    
31. 
    

复制代码

$F_n(x)=\frac{1}{6} (x (\frac{2^{1-n} (-\sqrt{4-3 x^2}-x)^n}{\sqrt{4-3 x^2}}-\frac{2^{1-n} (\sqrt{4-3 x^2}-x)^n}{\sqrt{4-3 x^2}}-(x-1)^n+(x+1)^n)+2^{1-n} ((-\sqrt{4-3 x^2}-x)^n+(\sqrt{4-3 x^2}-x)^n)+(x-1)^n+(x+1)^n)$

wayne

f(123,321) =8720560027328665154658721857141127550426817444170484901440812570613486430966452613389551773983777308244267241308

1. Block[{n=123,m},m=FromDigits[Reverse[IntegerDigits[n]]];
   
2. poly=First[MatrixPower[{{x,2,-x+x^3,-1+2 x^2-x^4},{1,0,0,0},{0,1,0,0},{0,0,1,0}},n-4].Transpose[{1+4 x^2+x^4,2 x^3,1+x^2,0}]];
   
3. Monitor[Sum[(k/2+n)!/((k/2)! n!) Coefficient[poly,x,m-k],{k,0,m,2}],k]]

复制代码

wayne

Mathematica 不支持多项式的有限域运算, 所以稍微大点数 计算比较慢,计算f(12345,54321)花了3分钟, 是一个13870位数,对1234567891取模是 232961944

1. 7799285511010828586160625392294329777217868025<<中间省略13777位>>44180209438908558808942353545069000162770988000

复制代码

wayne

自己实现了一个模运算,性能提升了好几倍了,但是要计算f(123456789,987654321)还差很远, 要用C语言才行. 里面的阶乘也应该换成模运算.

1. MatrixPowerMod[a_,n_,modulus_]:=Module[{b=a,d=IntegerDigits[n,2]},Do[b=PolynomialMod[b.b,modulus];If[d[[i]]==1,b=PolynomialMod[b.a,modulus]],{i,2,Length[d]}];b];
   
2. Block[{n=12345,m,modulus=1234567891},m=FromDigits[Reverse[IntegerDigits[n]]];
   
3. poly=First[MatrixPowerMod[{{x,2,-x+x^3,-1+2 x^2-x^4},{1,0,0,0},{0,1,0,0},{0,0,1,0}},n-4,modulus].Transpose[{1+4 x^2+x^4,2 x^3,1+x^2,0}]];
   
4. Monitor[Mod[Sum[(k/2+n)!/((k/2)! n!) Coefficient[poly,x,m-k],{k,0,m,2}],modulus],k]]
   

复制代码

wayne

搞到直接关于结果的递推公式了,是一个7阶差分方程,速度奇快. 前面的 计算f(12345,54321)需要3分钟,而这个只需要0.1秒钟.
$3 a(k+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) \left(92 k^2+44 k n+286 k+5 n^2+65 n+252\right)-64 (k+3) (k+4) (k+5) (2 k+7) (2 k+9) (2 k+11) a(k+6) \left(52 k^2+32 k n+504 k+5 n^2+136 n+1203\right)+4 (k+3) (k+4) (2 k+7) (2 k+9) a(k+5) (2 k+n+11) \left(2312 k^3+1660 k^2 n+27540 k^2+414 k n^2+12252 k n+109870 k+33 n^3+1437 n^2+23087 n+147123\right)-a(k+2) (2 k+n+5) (2 k+n+7) (2 k+n+9) (2 k+n+11) \left(3616 k^4+3392 k^3 n+29856 k^3+1104 k^2 n^2+20880 k^2 n+95288 k^2+152 k n^3+4392 k n^2+44624 k n+138600 k+7 n^4+294 n^3+4565 n^2+32958 n+77112\right)-4 (k+3) (2 k+7) a(k+4) (2 k+n+9) (2 k+n+11) \left(3560 k^4+2872 k^3 n+46852 k^3+870 k^2 n^2+27192 k^2 n+234142 k^2+112 k n^3+5253 k n^2+87604 k n+526583 k+5 n^4+329 n^3+8101 n^2+95911 n+449574\right)+a(k+3) (2 k+n+7) (2 k+n+9) (2 k+n+11) \left(13120 k^5+11552 k^4 n+176352 k^4+3776 k^3 n^2+121888 k^3 n+962192 k^3+560 k^2 n^3+29088 k^2 n^2+490952 k^2 n+2661864 k^2+34 k n^4+2844 k n^3+75846 k n^2+894164 k n+3730728 k+n^5+83 n^4+3653 n^3+66973 n^2+620706 n+2117304\right)-9 a(k) (2 k+n) (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)+256 (k+3) (k+4) (k+5) (k+6) (2 k+7) (2 k+9) (2 k+11) (2 k+13) a(k+7)=0,$
初始项是
$a(0)=0,a(1)=0,a(2)=\frac{1}{48} \left(1-(-1)^n\right) (n-1) (n+1) (n+3),a(3)=\frac{1}{768} \left(1-(-1)^n\right) (n-1) (n+1) (n+3) (n+5)^2,a(4)=\frac{\left(1-(-1)^n\right) (n-1) (n+1) (n+3) (n+5)^2 (n+7)^2}{30720},a(5)=\frac{\left(1-(-1)^n\right) (n-1) (n+1) (n+3) (n+5) (n+7) (n+9) (n (n (85 n+1749)+12083)+26355)}{185794560},a(6)=\frac{\left(1-(-1)^n\right) (n-1) (n+1) (n+3) (n+5) (n+7) (n+9) (n+11)^2 (n (n (31 n+615)+4361)+9345)}{7431782400}$

1. Block[{n=123,m},{m,n}=Sort[{n,FromDigits[Reverse[IntegerDigits[n]]]}];Last[RecurrenceTable[{-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) a[k]+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) a[1+k]-(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) a[2+k]+(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) a[3+k]-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) a[4+k]+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) a[5+k]-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) a[6+k]+256 (3+k) (4+k) (5+k) (6+k) (7+2 k) (9+2 k) (11+2 k) (13+2 k) a[7+k]==0,a[0]==0,a[1]==0,a[2]==1/48 (1-(-1)^n) (-1+n) (1+n) (3+n),a[3]==1/768 (1-(-1)^n) (-1+n) (1+n) (3+n) (5+n)^2,a[4]==((1-(-1)^n) (-1+n) (1+n) (3+n) (5+n)^2 (7+n)^2)/30720,a[5]==((1-(-1)^n) (-1+n) (1+n) (3+n) (5+n) (7+n) (9+n) (26355+n (12083+n (1749+85 n))))/185794560,a[6]==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)))},a[k],{k,1+Floor[1/2 (-1+m)]}]]]

复制代码

mathe

由于取模后每个整数小于1.3G,可以用32bit有符号整数表示，我们可以使用2^96进制的大整数来表示一个多项式，于是一次n次多项式乘加运算可以使用一次大96n比特大整数乘加表示，计算结果大整数每位肯定不超过1.69G^2*2n,在n小于0.13G时远远小于2^96, 然后对这个超大整数每位再模1234567891做一次规约，为下一次乘加做准备。
这样经过一次大整数运行我们可以进行一次有效的模1234567891上的多项式运算。
这样一个多项式最多需要12*0.13G<2G内存。
而我们需要计算(-x+y)^n,其中y^2+3x^2=4,
使用快速模幂算法，设其中任意一步结果形如u(x)+v(x)y,一次乘法会出现y^2项，替换为4-3x^2,可以回归上面形式。
于是可以快速算出，这样方案应该可以在16G内存以内机器上解决问题

wayne

PARI/Gp花了10分钟，计算f(987654321,123456789) 模 1234567891 的答案是 428074856.
如果是C/C++语言,可以更快.

mathe

wayne的公式很好，转成C代码可以4分钟算出10位数
time ./bn 1123456789 9876543211
898721522

real        4m5.255s
user        4m5.229s
sys        0m0.016s

wayne

C++调用GMP，挂机跑11位数，用了五个半小时。f(10987654321,12345678901)(mod 1234567891)  = 0， 计算结果为0 提醒我们模1234567891取小了。

1. +256 (k+3) (k+4) (k+5) (k+6) (2 k+7) (2 k+9) (2 k+11) (2 k+13)

复制代码

O(n)时间复杂度，O(1)空间复杂度，非常稳定的每3-4秒钟计算10^6项。

mathe

我们的计算模式是模一个大素数P=1234567891来进行。
wayne的公式在递推计算a[7+k], 需要在模P有限域除以256 (3+k) (4+k) (5+k) (6+k) (7+2 k) (9+2 k) (11+2 k) (13+2 k)
其中k的最大值为(m-1)/2-6， 正好在m=P时第一次出现表达式为0，所以会遇到除以0问题，导致无法推导下去，同样后面对m=P+2,P+4,P+6也会遇到问题，
这导致我们无法利用wayne公式递推计算这四项。而如果有方法绕过这个问题，那么后面会到2P+1,2P+3,2P+5,2P+7又会遇到问题。
我们现在查看\(F_P(x)\), 由于组合数\(C_P^t \pmod P\)仅在t=0,P时非0，我们可以得到
\(F_P(x) \pmod P= \frac 16 (x( -2^{2-n}(4-3x^2)^{\frac{P-1}2}+2)-2^{2-n}x^P+2x^P\)
由此，我们还是应该可以在O(n)时间内计算出每项a[P]的。
而类似对于a[P+2], 我们可以根据杨辉三角，对应组合数仅头尾三项非零，所以对应\(F_{P+2}(x)\)表达式也比较简单，类似对于a[P+4],a[P+6]也能计算。
同样计算a[2P+1], 我们需要二项展开\((1+x)^{2P+1}=(1+x^P)^2(1+x)\)同样非零项比较少。只是越后面越复杂，计算量会越大