堆叠杯子的趣味难题 ProjectEuler#907

mathe

根据代码里面递推式
DP[u]=DP[u-1]+getUP2(u-1)+getUM2(u-1,DP)
根据getUM2(u-1,DP)，如果我们定义DP[0]=1,那么getUM2(u-1,DP)就是DP[u-3]
所以
DP[u]=DP[u-1]+DP[u-3]+getUP2(u-1)
而getUP2(u-1)在(u-3)除以4的余数为0,3时是1，也就是u除以4的余数是2,3时是1，我们得到递推式
DP[0]=1,DP[1]=1,DP[2]=1
\(DP[u]=DP[u-1]+DP[u-3]+\begin{cases}1&u\equiv2,3\pmod 4\\0&u\equiv0,1\pmod 4\end{cases}\)
于是DP[.]每隔四项就很规律了，属于线性递推，比如好像
DP[4k+4]=5DP[4k]-2DP[4k-4]+DP[4k-8]+3
推广了一下就是
\(DP[4k+h]=5DP[4k-4+h]-2DP[4k-8+h]+DP[4k-12+h]+\begin{cases}3&h\equiv0,1\pmod 4\\0&h\equiv2,3\pmod 4\end{cases}\)
于是可以DP[4k+h]的通项公式会很简单，形如
\(c_1(h) x_1^k+c_2(h) x_2^k+c_3(h)x_3^k+c_0(h)\),其中\(x_1,x_2,x_3\)是特征方程\(x^3 - 5*x^2 + 2*x - 1=0\)的根 , h=0,1,2,3

另外UP[.]的值总是0或1，同样以4为周期。最后s相当于将DP和UP都是从前向后每4个里面规律抽取两个，右从后向前规律的每4个抽取两个，最终应该可以用公式表达。
只是从计算角度，还是直接用DP的递推式和UP的公式即可。
另外上面$x_1=4.6134702675815553806829461901509264729...$, 由此我们可以推断$S_n =O(n x_1^{n/4})$

iseemu2009

我发现   s(k+3)+s(k+4)-s(k+5)=s(k) 对于50%的情况都是成立的。
上式不成立时，有的情况需在后面加2k，有的要减 2k 就成立了，
即 s(k+3)+s(k+4)-s(k+5)=s(k)+2k， 或 s(k+3)+s(k+4)-s(k+5)=s(k)-2k 。
上面三种类别加起来就是 100%的情况了

iseemu2009

数据验证，希望找到通项公式。

iseemu2009

找到规律了，如果 s(k)中的k刚好是4的倍数时，或 k±2是4的倍数时，
                   s(k+3)+ s(k+4)- s(k+5)=s(k) 成立
如果 s(k)中的k-1是4的倍数，s(k+3)+ s(k+4)- s(k+5)=s(k)+2k 成立
如果 s(k)中的k+1是4的倍数，s(k+3)+ s(k+4)- s(k+5)=s(k)-2k 成立
四种情况各占25%的数量

mathe

另外根据9#的代码，我们可以得出
S[u+4]-S[u]=2(DP[u-2]+UP[u-3]+UP[u-4]+DP[u-5])  （循环里面u+4比u多出的部分）
      +2(UP[u+4]+DP[u+4]-UP[u]-DP[u]) (双方初始值差值)
而其中UP[u]的取值很简单
UP[u]= 1 (1,2 == (u-3 mod 4))
       0 (0,3 == (u-3 mod4))
显然UP[u+4]=UP[u],所以简化为
S[u+4]-S[u]=2(DP[u-2]+DP[u-5]+DP[u+4]-DP[u]+UP[u-3]+UP[u-4])  
根据11#的结果，T[u]=S[u+4]-S[u]也可以根据h不同写成4个通向公式，或者说S[u+4]-S[u]也可以适用DP[u]的递推公式。
现在我们计算前几项，可以得到
\(T[u+12]-5T[u+8]+2T[u+4]-T[u]= S[u+16]-6S[u+12]+7S[u+8]-3S[u+4]+S[u]=\begin{cases}-12&u\equiv 0\pmod 4\\0&u\equiv 1,3\pmod 4\\12&u\equiv 2\pmod 4\end{cases}\)
对于$u\ge 2$都成立。所以我们只需要提前计算出S[1..17]，然后就可以用递推公式了。同样通向公式也容易计算，只是最好根据u模4的余数分成四类计算。
形如\(S[4k+h]=c_1(h)x_1^k+c_2(h)x_2^k+c_3(h)x_3^k+c_4(h)+c_5(h)k, 0\le h\lt 4\)。
而iseemu发现S[u+5]-S[u+4]-S[u+3]+S[u]比较特殊，同样我们可以检验前面17项，可以发现\(S[u+5]-S[u+4]-S[u+3]+S[u]=\begin{cases}0&u\equiv 0,2\pmod 4\\-2u&u>4,u\equiv 1\pmod 4\\2u&u\equiv 3\pmod 4  \end{cases}\)。这说明这几项放在一起正好可以将各自系数\(c_1,c_2,c_3\)部分抵消。
比较上面两个递推式，根据特征方程，可以S[u+3]-S[u+2]-S[u]也比较特殊，计算可以得到
\(S[u+3]-S[u+2]-S[u]=\begin{cases}2&u\equiv 0,2\pmod 4\\u+1&u>4,u\equiv 1\pmod 4\\-(u-3)&u\equiv 3\pmod 4  \end{cases}\)
进一步简化可以得到
\(S[u+7]-S[u+6]+2S[u+5]-3S[u+4]+S[u+3]-3S[u+2]-S[u]=8, u\ge 2\)

王守恩

1. RecurrenceTable[{a[2] == 2, a[3] == 6, a[4] == 12, a[5] == 16, a[6] == 22, a[n + 5] (n - 1) == a[n + 4] (n - 1) - a[n + 3] (n + 1) + a[n + 2] (2 n) + a[n] (n + 1) + 4 n}, a, {n, 50}]

复制代码

{2, 6, 12, 16, 22, 36, 58, 82, 114, 174, 266, 382, 548, 816, 1212, 1762, 2566, 3780, 5560, 8128, 11892, 17454, 25604, 37498, 54932, 80538, 118062, 172996, 253510, 371574, 544600, 798112,
1169658, 1714260, 2512406, 3682066, 5396294, 7908702, 11590806, 16987102, 24895768, 36486576, 53473720, 78369490, 114856026, 168329748, 246699284, 361555312, 529885016}

a(100)={105810375074960304}
S[100]=105810375074960304
a(1000)={26962901013997007773476449840406160338483621262752268591326543606162290720748003129661220722119075977873544416448534677592402973194126252297783679735507922650088522700}
S[1000]=26962901013997007773476449840406160338483621262752268591326543606162290720748003129661220722119075977873544416448534677592402973194126252297783679735507922650088522700
a(10000)={31127930776577878223912939652283917007309720910229733404980767096433599906671652710001938387364854655740186734590908649337356091922621078920216052890139240853978834431849749301441120982004337938233139987156772948822431045256171786679546693883324545906104252627630400702096311093294706718040576545408214926280868115722806877605932122196052494082623011237601970298979931170969289976068813021526033675831059992990730577890989631841080978643700990585349553415808578723546946916820500168052783156513103111742386023998708284150142291655601723334559334350809797047051453261439348573046522312830382200517486230343043897711784893663693727639779892990757401519165412641082883760609113920300639590955495404532854649163537655214920643816932518331623066331271805714815696682277895216235139841994069819511789713982835479855092133910403059302059116100776409354150790798813753273844790430066915461356063289968467220558523526513213228886433823912526799122628736929131922922546052623254928685742281834365448375330946022402419007606452114909224699039372244294796852125356356305365277108241592981671643800362193173209313411961857120583215000544762786535022200235646986267655715861136663687268833584319373370806212277764365165024210936127505563077780044446172967341034465955736157812639997021613611139824073413196090993637257782802715278260197652570802957241392696659896441870090217096364557746038056749198966414529129939068897927148170015396744856971653040759512206663129263689011235493838837509291180058136049825236927951726835961962015206904611846148078834006589220003421416282780956845528315528806900524316290074905895950832826283918972260149280538119730487802435228313647395336}
S[10000]=
{31127930776577878223912939652283917007309720910229733404980767096433599906671652710001938387364854655740186734590908649337356091922621078920216052890139240853978834431849749301441120982004337938233139987156772948822431045256171786679546693883324545906104252627630400702096311093294706718040576545408214926280868115722806877605932122196052494082623011237601970298979931170969289976068813021526033675831059992990730577890989631841080978643700990585349553415808578723546946916820500168052783156513103111742386023998708284150142291655601723334559334350809797047051453261439348573046522312830382200517486230343043897711784893663693727639779892990757401519165412641082883760609113920300639590955495404532854649163537655214920643816932518331623066331271805714815696682277895216235139841994069819511789713982835479855092133910403059302059116100776409354150790798813753273844790430066915461356063289968467220558523526513213228886433823912526799122628736929131922922546052623254928685742281834365448375330946022402419007606452114909224699039372244294796852125356356305365277108241592981671643800362193173209313411961857120583215000544762786535022200235646986267655715861136663687268833584319373370806212277764365165024210936127505563077780044446172967341034465955736157812639997021613611139824073413196090993637257782802715278260197652570802957241392696659896441870090217096364557746038056749198966414529129939068897927148170015396744856971653040759512206663129263689011235493838837509291180058136049825236927951726835961962015206904611846148078834006589220003421416282780956845528315528806900524316290074905895950832826283918972260149280538119730487802435228313647395336.

wayne

mathe 发表于 2025-2-22 22:28
S[3]=6
S[4]=12
S[5]=16

拿这个数据,跑一下,发现Mathematica返回了一个常系数的线性公式.

1. FindLinearRecurrence[{2, 6, 12, 16, 22, 36, 58, 82, 114, 174, 266,  382, 548, 816, 1212, 1762, 2566, 3780, 5560, 8128, 11892, 17454}]

复制代码

$a(n+8)=2 a(n+7)-3 a(n+6)+5 a(n+5)-4 a(n+4)+4 a(n+3)-3 a(n+2)+a(n+1)-a(n)$

1. RecurrenceTable[{a[n+8]==2a[n+7]-3a[n+6]+5a[n+5]-4a[n+4]+4a[n+3]-3a[n+2]+a[n+1]-a[n],a[2]==2,a[3]==6,a[4]==12,a[5]==16,a[6]==22,a[7]==36,a[8]==58,a[9]==82},a,{n,1,100}]

复制代码

生成函数是$GF(x) = -\frac{2 (x^7-4 x^6+6 x^5-8 x^4+7 x^3-7 x^2+3 x-2)}{(x^2+1)^2 (x^4-x^3+x^2-2 x+1)} = 2(\frac{1}{x-1}+\frac{2 x+1}{x^2+1}-\frac{x}{\left(x^2+1\right)^2}-\frac{4 x^2-x+2}{x^3+x-1})$

通项公式就是$a(n) = \frac{1}{62} \left(a_1 r_1^n+a_3 r_2^n+a_2 r_3^n)+(n-4) \cos \left(\frac{\pi  n}{2}\right)+2 \sin \left(\frac{\pi  n}{2}\right)-2$
其中 $-4981824 + 84816 a - 496 a^2 + a^3=0, a_1=164.52306508496125794..., a_2=165.73846745751937103 - 53.02036465022101892 I, a_3=165.73846745751937103 + 53.02036465022101892 I$
$ r^3-r^2-1=0,  r_1=1.4655712318767680267,r_2=-0.23278561593838401333 - 0.79255199251544784833 I, r_3=-0.23278561593838401333 + 0.79255199251544784833 I$
在$n>12$时,有一个简单的取整公式,$a(n)=[\frac{1}{62} a_1 r_1^n]+(n-4) \cos \left(\frac{\pi  n}{2}\right)+2 \sin \left(\frac{\pi  n}{2}\right)-2$,

王守恩

这样也可以。这2个是等价的。

1. LinearRecurrence[{2, -3, 5, -4, 4, -3, 1, -1}, {4, 2, 6, 12, 16, 22, 36, 58}, 100]

复制代码

RecurrenceTable[{2 a[n + 7] - 3 a[n + 6] + 5 a[n + 5] - 4 a[n + 4] + 4 a[n + 3] - 3 a[n + 2] + a[n + 1] - a[n] == a[n + 8] ,a[1] == 4, a[2] == 2, a[3] == 6, a[4] == 12, a[5] == 16, a[6] == 22, a[7] == 36, a[8] == 58,}, a, {n, 100}]

王守恩

王守恩 发表于 2025-2-23 16:46
这样也可以。这2个是等价的。

冤枉啊！！！看看我是怎么找到的！！！
Solve[{2 a + 6 b + 12 c + 16 d + 22 f + 36 g + 58 h + 82 j == 114,  6 a + 12 b + 16 c + 22 d + 36 f + 58 g + 82 h + 114 j == 174, 12 a + 16 b + 22 c + 36 d + 58 f + 82 g + 114 h + 174 j == 266,
16 a + 22 b + 36 c + 58 d + 82 f + 114 g + 174 h + 266 j == 382,   22 a + 36 b + 58 c + 82 d + 114 f + 174 g + 266 h + 382 j == 548, 36 a + 58 b + 82 c + 114 d + 174 f + 266 g + 382 h + 548 j == 816,
  58 a + 82 b + 114 c + 174 d + 266 f + 382 g + 548 h + 816 j == 1212, 82 a + 114 b + 174 c + 266 d + 382 f + 548 g + 816 h + 1212 j ==  1762}, {j, h, g, f, d, c, b, a}]
{{j -> 2, h -> -3, g -> 5, f -> -4, d -> 4, c -> -3, b -> 1, a -> -1}}

王守恩

LinearRecurrence[{1, 0, 1}, {4, 2, 10}, 100]
{4, 2, 10, 14, 16, 26, 40, 56, 82, 122, 178, 260, 382, 560, 820, 1202, 1762, 2582, 3784, 5546, 8128, 11912, 17458, 25586, 37498, 54956, 80542, 118040, 172996, 253538, 371578, 544574, 798112, 1169690, 1714264,
2512376, 3682066, 5396330, 7908706, 11590772, 16987102, 24895808, 36486580, 53473682, 78369490, 114856070, 168329752, 246699242, 361555312, 529885064, 776584306, 1138139618, 1668024682, 2444608988,
3582748606, 5250773288, 7695382276, 11278130882, 16528904170, 24224286446, 35502417328, 52031321498, 76255607944, 111758025272, 163789346770, 240044954714, 351802979986, 515592326756, 755637281470,
1107440261456, 1623032588212, 2378669869682, 3486110131138, 5109142719350, 7487812589032, 10973922720170, 16083065439520, 23570878028552, 34544800748722, 50627866188242, 74198744216794,
108743544965516, 159371411153758, 233570155370552, 342313700336068, 501685111489826, 735255266860378, 1077568967196446, 1579254078686272, 2314509345546650, 3392078312743096, 4971332391429368,
7285841736976018, 10677920049719114, 15649252441148482, 22935094178124500, 33613014227843614, 49262266668992096, 72197360847116596, 105810375074960210}