连分数-ProjectEuler 946

mathe

试着用gp计算了一下，到$10^6$项大概花了一分半钟，得到中间结果挺可怕的
c=1000000
        s=6318511
        b=(2998812662783403143011355960986399961513893334537115366346166076040057022739460207773155459432929573293143717401653759958776970892108442485281267348348903399745107380121149019463318993639461427990715068959422829395215443413478943645464540634957621616287403834619745581224800579214496838645951779546937587341280914482849777862065562337016727143478454170778337875512064760572095882383755297441017080135706139205932840920432121819466682338636494351010027372184225789767617073900774958175842832533834024559359655925895793*t + 7850993477060526277261707522214466439833349348919927003878809714906149887159356208966001682468724602546044505004281681381920590746772472103445420813990251748565362068896281490713447714213409325423238432942698109121601758470308088527792075073008012685287158737652066391551423958684417797713825646853723822836360732345261599442373654093998305348760232716685354631640949453132140597845042022886597571549577466697954214502623500408891397101195437591951224753368006624100678023163749446389912370886510461106037089088275650)/(270402769118390527160910118249132699486133395617135256893889586010008140183798343778987439715538559919354198799909527445502332136452608208570628541415514546140414402906470200544091653754048063236896268240993763470208027976314175247921347110590982414699543838203639401834712263633630086171755100158430478745151773449260885005015960088774073878581259363308380776547905378223456518236391247699765815983707523812644140547499918391390115103702761396625767928445097853524368155108672775663605031492709986258509918404732941*t + 707923640204036839466635239496933072097232025154277176026789049612014660621411379066226067772801338626370423000294562969956327033002731780484504886235986498405187237186817884110300460229511510892901521088282071267126299873221458828261446971227917289137298302205404876978845021296351402197892395972759708966316892276274309574419496502912974980128194579479019885721575103403801783771307558568109742420624237366244224403086391720449382509208403428852312645814884299363270676823624101900610579103850799415275898547072335)

...
将近一小时，计算了8000000项
c=8000000
        s=49573674
        b=(362084373216698230555512285247751509930745478930759430356050873448246554226392137694023136992256721876951973069938427562719582628882428170368121117853252168378105871670661529911283317937019440669235604693154550835252322446241351755102717363171694574725596745762231747298037596474729985190527903109110060994857125324409048883274910345686164504748472386025608551365845472766964687569586811315087907954547626922066973724884863089889507292992804346081877094769231729583830353493215079473101757846645673979772137670969832824815798078220314652477919366243068590988566249909782938079387398910265623007761149346105110401961416986476355682458381729145300476998084594580263772922773761039566863675127831305370229907586622079772335006965345753947225853910679165711197*t + 947949195876518062824380169884129829886652242456048354196735201274522353892511233642233982634224554823626964900950573077615172850932874644286781433348332546493493028184937786020391383363168498991278442336497344627252091410761298731622525481178517204944132627889958554753360044989811455958391428842536193753369251609358175251488081395891156233877380514722400561410772825082290380572506421722726386224813961425968080290694483388964335767757936121868188590355583223145631205162367415280046096529531708851766817873761461716779699753810787160963151172308379506934170587891910592491636011863893879505887240302072822010718637949748752121823920115074075342104043124390154227674882873172027173921291452993785066682488857551673204258521400390715201812561452296086581)/(32649127565513316900241141199879520333569986180369176355014651545406564519458378188869132271785389717759174595553307072721290079192807006648332634649384543552088809206089041255273783634589766256329350671477562277980265252668920841777301684660556553287905444540780058892311470777567472315094403794250175833703628829473240774866572485535642607656761555279091639297872979440910964457903561517489127714950657142969743029969196776258084029883224110276540719777903631138361430400041043249907630641264782521796299646199481861293123752189845322602804100382005851819491690632374767426143873486016298428541260010997674828094768924551136737160638096481470509545115892407342090378355472345766641829795457999112344613818200721325666202057507152262072961968432326187344*t + 85476525669544972871198913529502110119577090258299556613216035307811936653061778814352280307165500262494064452209008341632432557453608436732957277430404251096456698016567922484366177041239794636378861207221935078242570590315838687731464146498561536260281526968758373016532107605972985964143342139522082532452876596671106086601611063351489944012399099316758365992318286329757323297078786870013232676889788241601291956964651328481135070551850513326991706865007798372110642858659433194583403825835047804672947392673759309480409116520001059536705513401416058925509926234906494586335029176459643972968733873620092011630190967115725831355931584421648299217663275878989590735896941146220029466943749460718983827044604704484733482216072011694274355569508763617177)

wayne

根据 https://mathworld.wolfram.com/RegularContinuedFraction.html 得知一个算法,叫做 Gosper's algorithm, 能根据x的连分数,计算$\frac{ax+b}{cx+d}$
https://perl.plover.com/yak/cftalk/INFO/gosper.txt
https://math.stackexchange.com/q ... fraction-arithmetic

wayne

我们可以借此题总结和推广一下一大类问题.就是 关于连分数的运算规则. 反正总体上就是 基于 Bill Gosper 的算法为核心展开
1) 一元运算:已知 $\alpha$的连分数展开,如何求得$f(\alpha)$的连分数展开.  比如本题的$f(x)=\frac{ax+b}{cx+d}$,但是我们还可以讨论 $ f(x)=x^2, f(x)=log(x)$ ...
2)二元运算, 已知 $\alpha, \beta$的连分数展开,如何求得$f(\alpha,\beta)$的连分数展开.  其中$f(x,y)=x+y, f(x,y)=xy,f(x,y)=\frac{axy+bx+cy+d}{exy+gx+gy+h},$ ...

mathe

所以我对多个[1,1,1,...,1,2]进行折叠计算反而失去了规律。
这个题目我们如果不进行折叠，那么虽然计算步骤更多，但是更加容易找到规律。
我们可以以计算过程中任何时候输出整数部分m,并且将余下关于t的矩阵b作为当前状态，而每次输入连分数对应项数字1或2作为一个转移参数。
由于开始部分没有规律，我们可以先简单输出并且跳过，经过一段时间后，可以发现会有一些有规律的状态，
比如
b=5t+8  
  输入 1 -> b=(8t+13)/(t+1)
  输入 2-> ...
(8t+13)/(t+1)
  输入 1 -> ...
  输入 2-> (13t+34)/(t+3), 输出 11 -> (t+3)/(2t+1)
(t+3)/(2t+1)
  输入 1->(3t+4)/(t+3), 输出 1 ->(t+3)/(2t+1)
等等。
可以看到每个状态矩阵都很简单，而且由于外部输入只有1和2两种，最后派生出来的转移状态应该不会很多，把整个转移矩阵建立好以后，就很简单了。

mathe

这个是状态转移图，可以看到连续1的最大转移长度为6，所以对于\(p\gt 7\)以后就可以适用这个转移图了。
而对于P=11之前，计算可得输出为
1,5,6,16,9,1,10,
这时对应状态5t+13, 所以从这里就可以开始使用状态图了。（实际上P=7时就可以使用状态图了）
我们这时有11个1可用，根据状态图，先使用3个1，输出16，转移到(2t+3)/(1-t).
由于还余下8个1，在绿色图上转一圈输出11，然后转移到绿色点(8t+13)/(t+1)
然后使用2，输出11并转移到(t+3)/(2t+1), 由于接下去13个1，所以在蓝圈上走13轮，输出13个1，然后使用2转移到(3t+7)/(t+4).
从这里也可以看出为什么输出结果里面含有大量11和1了, 而计算结果已经不重要了。

mathe

现在有一个问题，\(\beta\)连分数展开后，里面数字的平均值是多少？
也就是设\(\beta=[b_0;b_1,b_2,b_3,\cdots,b_n,\cdots]\)
请问\(\lim_{n\to\infty}\frac{b_0+b_1+\cdots+b_{n-1}}n\)是多少？