From fbb62f01c58a233fd3b25e5fcf985f2c09abbf22 Mon Sep 17 00:00:00 2001 From: Nick Date: Sun, 28 Jun 2020 14:20:52 -0400 Subject: [PATCH] Centered continued fractions (#379) * Centered continued fraction [CI SKIP] * Document centered cfrac. [CI SKIP] * Unit tests for centered continued fraction [CI SKIP] * Kick off build. * Fix syntax error in docs [CI SKIP] * Fix ADL. --- doc/images/maple_cfrac.png | Bin 0 -> 32918 bytes doc/internals/centered_continued_fraction.qbk | 55 +++++++ doc/math.qbk | 1 + example/centered_continued_fraction.cpp | 46 ++++++ .../tools/centered_continued_fraction.hpp | 154 ++++++++++++++++++ .../math/tools/simple_continued_fraction.hpp | 2 +- test/Jamfile.v2 | 1 + test/centered_continued_fraction_test.cpp | 148 +++++++++++++++++ 8 files changed, 406 insertions(+), 1 deletion(-) create mode 100644 doc/images/maple_cfrac.png create mode 100644 doc/internals/centered_continued_fraction.qbk create mode 100644 example/centered_continued_fraction.cpp create mode 100644 include/boost/math/tools/centered_continued_fraction.hpp create mode 100644 test/centered_continued_fraction_test.cpp diff --git a/doc/images/maple_cfrac.png b/doc/images/maple_cfrac.png new file mode 100644 index 0000000000000000000000000000000000000000..4d362065ff3355ffb61ce3c5153f42dd2c357c44 GIT binary patch literal 32918 zcmeAS@N?(olHy`uVBq!ia0y~yV69+aV4T6h#=yYvgQ05~0|QS=rn7T^r?ay{Kv8~L zW=<*tgGcAo>Fg1~C&e0`PYV%>QeX;l%3tKctDxE<=IF3WK}(TSWr~DQQ`g+mNY?9m zN(s*4t|H~7tSP!3JFdBI{C>2lYm=T>&Wb&poNHe1`?UAxzMY?+pP6a=e9oCW*31ol z>hWr)T}v5HOcp!w_)v4?}W&$Jvu-l~dE zU0eU{|Mf-sBAP4=3Yim(TeuwMo`^AgNNV(nV_-03Ikh%jk917e z2Nbr=+Rpn*vch){<3-6KJKO&FDTO?jWtgHB=BBdhYhZ9=6ZdEK z853e#|4i^Wv24QOB|9$$N>#coQ|V!sTqjntY)QlP1(LfXmT>f5zp-YTe}A{+BOlA@ zizTDaPv<|k`ci)IGnFEfZ|l2V_bC2e-+hGhUBgYZcCds*AY&^SAPYEWv1HB>MQf#_OI9q5XjwENp!7||ZjV_jI>mQuYjgBn?s+@)>9MRYg_|_YV^p^{ z|J~&uxR-aO2TDLkW{qJmxcI6*R&!ud; zv{{$-EMN(g*(bO1?UJ{aIj&PL>VIDR%{19!^__?F)<1txrg~yS5)xY%EUN3Z$8TVYclg&kAH^p@y5%wCUaf4x}JBvf0kW< z!mP*zw-WBW+!R^lIc-P7(Z+a2!ESZY9=W4M#ihlp4qtr(_b-bxWZy5!@%c>rraxi~ zHC_u1Uh7NVmuhJ4XGqxd(>v|2P#*t}b@nxm3#N##O5O3wu_2#XcE{qLM#~Lsdp7Ppt;E9 znZkXA)Dv1ugl;xnQLW}AH ztPW0j;PXNC1IJ{J$QH>KAt%d3>xT>vS?>78wwSvpP4wKNUUm5H45u?{=hXC6^Az5x zZWFNX;P0qa5jr{JiO(nXqZ58CVG82Rw9Pbssqs>2W%rl9HLlwfc27KfVr7cv&KWU& z_gw9q_IcF#yj%VvC`WVO6p5}0M%|m{Z0fz_W9k2L(U;(_E53!)XEB-1-Zdf0aqaSX z3+^r3cj?llO_x4hoV7yrV$_1rgP%nrd=$%5S3w`vHnHL zE?L`-d0zSp^cP>h)PMQ?rTdx2TimlhZ)v^tdF$#e(p&Vi?5{6;apP4%$+A}>FQ>eI zbKP?N<;C6?w_bkw>Qm|SmtimEuIycQd-?GT(yPT6PhY-0B!Bw7%)RQrRBEmN2GlP6 zr}TI2uj60k4Gh{$n@`Tl^o|8}_DZu-Nw}R~^TPj;M&n@06-mkp2Jd=6SrNSiT z%-*CV`B5w{G56y&b4dne|l$EU=n*WYe@C-7pz`+&yuGzZ)Sha8ApFMo z{__dvtK;?WGu^MaKzmMfOAygBFYoaYDE9sKlw`C$IeeRp==nI@p4Vzb0#iEW6ziQyKb zFL^UUBf=tLBl35o|9JCZ8i$ih%8^qY**3=W8p2P%mVIwq_wL6wkEoQIl+-7GHI{2s zhJ0ozIrM{bwfddDZN1an%#Srl#GQFH!+PVljXyW;_fhfdnHypHG4Rff9hQaak0n1| zepG&RzbA)Z#JmFSi;96WZ_HXUdCIhq*GNhkHOW4?=D2I3tNoc{>D+1C z&uyHux!~O1opW|ZKU4f{^la<7?eq3$)Ln1-z-lcXB>qh2&)OSXD?&fM&Y!eDGcILY z#=%AG_uJ~D!xzOW>F@kD_kQ_ovB%x2-M6Rj(2v~zq(@Ag}7)o!2OR)0(J?!}$PdFjz>%X8oAU5-8d z`hrzz?DuzD?uvhFd6T)SeC7R3(|2xv6fF1akw0$ESZbyEF2HnHUQ8R0!>kp)GuG_;%q!+n z>o2SFVSYrdzKl*B!y67Jg+fM#zz;o}1jC|FnA}*^_lnbzp|o6n-K<3?O3H(8FFGjA z|M2dC1KP$69%mV9Ch_xp$(NA$!Sa5^4PJ(*Pgl!*LlqD0KA5I1I#)h@4?|7pmG2ky zlm6^t)#A4~aQ*FmhU=Hx7qi%2-NV4Z(3KkDnda-upvAzzz`?-6D8<0az`(%Bz`$U~ zC=F-3F={X{gT6EQV5(Uw7BIuu zj59zQ9fO@cJsB7nf|H95@=I8s@6ljjVBjq9h%9Dc5K{$VM$aIX4-5>9*E2&R zN+NuHtdjF{^%7I^lT!66atlD_FxXUBRpb`rrj{fsROII561Ai&D~Tl`=|73as??%gf94 z%8m8%i_-NCEiEne4UF`SjC6}q(sYX}^GXscbn}XpA%?)raY-#sF3Kz@$;{7F0GXSZ zlwVq6tE2?7NC5^Q?o6%7MA(#94E0uWey%=9M&D4+Kp$>4$as*bRX}D%YEFbpW^QU; zab|v=ouP?=fely#h7`g|8?e^U;?yGN{M_8syb?P^s6#a)OQLH8+Zl|?L93!CEvE#aQCJYpjGMO&^Aq)%*Odt*e1H%~x1_n@i z0MTa{7~Cojn@kO6U~pjYba4!+V0?3z^Nx(`++!cr*&7vfPAfhX5=d%b5o>HXI8}hT zgUvNnarOy8K^76Q3$9_|hA~ANO)`rXZQ#^q>19&lYLa=#C2?qJ=elDoE{T!NoScm< zEoaVM&OLkE_`S^C&*y&cEB^E9^GqNA^3u|)+tcoO{$IYU%0X|wixUGQi+}?bBEzzg zk)qi&l^Vfd@}*K{PhAM>S&! z85#zd!+zF^{h#N)f6;FLXW_p;&+T8^zOU?Gx9?Zh+THJV`JUZT@iFOc<@32!|9(E- zdYI2TX?AYlr#m*L=Me?f<^+`)<>$D=YGTw~6ZS z`OqX+@qke`YKzD1{Qb34>UY20#yvaRZhgc?r}H+Sd*)TW(%k#;n6%Z~Ethw_oxJql z*DuM_KE)ode!n+-|L?o+r<{Je<+7jk%LUDRm#ki|xxDsX+{B}wZ>GnZ-9hlBawx&&;=v z|NQLiVyPse+u z!>??+aO_guN%qF#>UTTS_2eVY$=3aNxU<*%p2e;A!SjEGpAC6>@@LN0erT3IBU|-i;l)qZ_dfQ`**SInRFJ^I?WwoQNMEc=sEr+cmD^?JPh^#2#NrvKk)H2r^}>i4Jm^_C(>izjpFL~c@< z{$J5@`}MeL*LvssU)Pn3W|?JPS`xGL{beULF_f~eQBCXv$6-6)<1)oQc^4Kq?)=kn z?!^wi336pO5_kT2an5e7`^$sx_Wgc$=bZKXIn^JJihG+s{rBto{%4<*Ba45(-(UaE zJN{RYo_*3*ft1glwRJ5i+iSJf793Eq-=#fWo-sy4Ik%zw>@3rrhRMf7zG;^`Mfv}H z*DavnW0$}6YM85I^n8Wg_h#qst9-Zf`MmB~KVGd~@1}R}{f@_d+F9Rfzu%o+eBk2W z+}(dPu4P!tJf403Pg>x=C*6^o(_D`i7JOgSt#@i|>~6E2J2NbE|Gs{Ee9P0@&vWyu zqxI}pS$&&Jf>{_?SMm?8gCM$|C%)Rc_7me1!1cRktJ$#N(B-|v}k!gpD6#NMBy;wo%#R1EKlt(J0j@r(s=iYTmLOC z$8FW`cB=1O`DzpM)Ve?azVFw6GvDmMYvG{&I^N^D;t$p?|MN%nXyAmq$~=?#`&s;! z>DPVizPM3;&xb=#Ud?&^J}{f_*}1va+c)Iy|NHH<<%1vf|G(F_mj29?h$;r~R)BPJJ%exHSlxuIC#jyFDqs66jvKdK2@VciaD@-f7K^wEcb~xqptYtWJr+ zS;nOwcD!6Rd(o%$yI!r@b^1K4Y}tI_!+!lGKlfN<-vTG-UrTiDpYWgG`fC1F@1&X+yL=QRFg+&=U6YyZaC zIh#DcY`Yem=Z|5ugr@b6`oidvrjw7t(} zWpn+W%Kf43@iN`}wGXWzxx5l^5rx&#x`JxUEQk)0O@6>dplHV-=4PxcYnV``Y~NtLk1j%Xe+N zTlf2IYi^{C?V+6ad%yR^aI>n2*=U_@RP#S65H#)6UVpLqRlm%=f3EWn^M5ZKbnCqL zt?$vcZ%bn@O_P5mJMZ!p`_<{!{zU%WJZ1jX&HDHLZ#(~~(*E-q`d`$`n96hEx`;y3LHL+f(64>5<9s&{3c`slGp z;K>@%2WRgeP&ZFH(sAd>$;pcw0&lF3YB10JsVMNUGTHX{+UJUT(`~rSWv)**n67sI z;On>6?{*|_Ir;fx*sA%>;pbjAERmi3^78zDC0la8-!1ny{pxz_mwNxim~~hCBTw-w zU%9(2d;Q*R1rq%N-&Yu(50}w>p}4^VUAdS!yIB!l$Qmn{q1KD{+sXn z*qg!km{VBIW%2s@Kab_x4?mN7`BI^~|D`lX>y;IsR{v~Co@y8V>i(sgy4baUbmpDh zbh$2W)4uzYj(g~ZXXWpGx_oW@Ccea-v(isYu30y2f875qPt`xGPS>tf`>Zlu`s=^* zpDMrKSF3;a=g0Z~f6k|rE_zz9wDrGSn|$f6)T{$tCie>;8t|?9f4%;1__=W5AMCdN z^B=F)oFiBLX5*dh`@XK7Q@`SK@1sL!o#m@uEWGpk{{Me3)TW<&nVNQfo-Au=^TFHq z|9yK=eK(k2HF#^X)%@R9r~fl&O`UGDZF$s(5U+hRXElG!TRlDYx%}T3?kCr@KK<~& z(aXL>%xDv8Bh_T4LqF4=eL2_6`#zmL`EJMKz85an9ChbkX%n_R|M-Ql1^b?wq|!^C z>aOc$4{$HG+k1~OjJ?)&`s#S~_XW9|e|~=Me$aUSil0x;EIxhj*R}0#m3#ePK3$mk z>RFNPrBmidy=Kc;{C)rb-+PUO=cgaNmfZU~c`4@stNHpbx6i%(`pOB0_4)rw@0MPF zY9re8yrXpg&3pX*%O~=z`n6+WzKzyKleyLJcDB}LpS6BzSyS6{$@=X4f|2ytt;^r&+Gmii1Hu)Rj|8w5;kMjQ>%5%N_&|59@Pr{s&wpP4n2-g@eVzq+3?Z_Cr=UwKl~f9^@C-QSu##V$N6-}Y(x%KA-mnYM@R zGC!w(Js(wfVt&oKCo|7Xj>)_7+y3YeOVZrThZfiV# z%RW!8*{`KK7zy6k!=AB=@-?4hVMmYB2qp~0EHzwG~ zw(mdOd{BJfht@gyvtG(8IWO>;zdzN@*QlDy+@|>Vnfd;cm39|n%Wf{s_}|QD zrE%)@&*I6N-gj-RPEJo*=4mJR%`@3z<=4&b>3VT{EY7}jt;`WUwbG>N`t+N(KYyP8 zzb5g%8(Rf;>%ZEA#jimHcWc$(+x!1U2Tt4j?UwefX9?Lw>%!By)L-t~8L;)q#l7FY zg~%>7-^{r- zeHx!%d#G-!&E%hlit_jTcfNVm?!DYiDZBfZm-%w<7MXY0BjEYTBH?pRjHY6nYK~j( zc)8lwX8X@eOTD>O-)}x|7bq?Bq4bj|+q(C8yWfhP-MsH++HBWiq2KP&>%X74pTea) z&-_-#;=Ev~CAP~>e=gsDdi}Tm4LVONOSaZ?eEaV6$$h!W*tH&aM6Q#34Sj;8L(jh028I4ZJ6(xBgmk^7AUS$;}mU?=_{j z-Ok&6=l`Gc|DRl26B%5xSv~&VtG_n8zI3Y3^LRH;+BnT;OMQ__#Tw4{eNu?`}r=szhhf0$A_KwW&Tm))6bU%zTfWkcJ1;%I$zgczx?~_wE7jn(>@)xd$@B_ z`i;0>cckC{e`2Urd(YO@@yG1&yi-!~y+Yf}e--`LxnKCt@oIhhpT#rlHF!^kAKX#* z@1Rjlc>B*qvUh5ZTi*-(kvHEc)oaT2%gU<3KmN_LC`>9U{i2+6E4k`;fZx*p>?xVe z`{n=seY!esUgfiyi>g{zF0zT*8KpJ%{GO<|ZHg%IbXF-Kyy5#{`F{l^ho*nydkE@w zOxbTz@Oi{T`X71u&+nFIVwdZD{#C9tWLDPH_*wQpZ%yBB zRq_7u523dyzkZ&tpZ9si+W_|ddtPz3qTkv+WN3DuT@}@H>51UMg=apWx8Gkm_CI~~taR%s|JO6We!kT4J$mc6G+q1bzpuaQ++82|Z|T#v>!5D- zpO~}8Gxc}gx&6Lge)XN_wHFt;wx0i5zW?{#l&OnS*K^mex0&js96as)4p2+I|5f3Y zKzG;B59O`v&lKF1xwEy>EB0S)MkGjni0u z=|A(PDbN3_?^ws@r9Qd7q~O7|+ctaMwx+aS`uTA6({yRn<_izgOy)i3?SAJZE?N3~ zwfE$jN6&<{`?qnuJSrZq)5G?aGkfQgz=C;?um7Ibcdzupsn@op&Jv9w{&!>CY#HOO zo8RF%ky;TFIxR$-;kESL9d}InpPbEilFgbb|0%TU-}L&Q(_h?Js{U2+|Djgy;&=7( z>OSWmT->{_WJgR9!zJ}|de?tXQDcAO9B^~-53kAYzU#f;a|E&3{|6Jl)=qH}r-t)Deo_`b{mG-H0>8Ht+ zS(R#^+n(;L`v1FP>HPBhZ?t317JNJ^-mWD&Z~w2W>)qy8KKJzN4mc6;+-cf-i~AbJ zXV~A#&*HOwvtiEi*bny;&r82qknVAEwo#Scw+Vl?U#q|T$9cv6e&aTi-X`gIv-3s^n;Wgn>W&hvBi*7$$P_{Nj-rfUc(9uDgGllQNzYp#9J$_Sa znm3zhvT)X1db;zY)2EKu`2{~OU3|isbmQG~1M~Usf6Vonx9Py^J^S-F2nSu;RcSvx zFu&pY^yXhJmy|WU-P78vLaWw^Hhg`$YTmqsUj*HFci&TX>*?4weK*UFu%A<3iC&vu z^XX*Z?n#af8+3oZkzUn$;LSyM`MGyOBDVGF&$fNgz`V-p?cVSA0_$b!|9ot%J$z{I z_j}U0#k;fb9y$Aa$=dy&uU+ncdg}X$OR1pAk6Y<=>refM_j6G6I5Kl1#x-thi$ z<(KqZQxAbz1@ zCoJ7_T7SRIs^FLZ#s55B{on6b;VI4KDyyT*W-{bQiiIg56>!0Vfj9YKVf`UWe{=To@zd4ex@7yTrN-whlvnjX z&xzTV6S@7(y7vd!<>!3ntv$usEno4V(dokXBf|bYr-~jdG5P7{nz#LVna|`)*Y2lY z1qXjW@4CmX;;~oOcI|s|^R~_TtgRuZjd!ls-z{fTQLw6^{D!{w&&2szc9Fl#r-c8C zS^3AUBJB5`E1&C5U%FRyN`LFhu>Y6#JS`dUsQj~CB9k^>+#U#!=2*uysH;A z%zOT`SaPbHh-EQ{+s{u=C&*}gev$b8^+eBoI+vC0YCH3<9#5}V2!ma*GC2zN0 zcl)YQe%Rpr!50@j`SPZC%eBhLOfKe_etLGybh8JOm;KpU{QOzX=DC-wcHMJs<2h;j z^U36s&gIeHIt;cSEdDfYpZJ`c)35i*+uQ9*?9VrN#;{&~p7Zr!%O97XB_3`QtrMtt zxAQsIbdl;WWk#*iA7e_dhPGafd_8}w{N`!zE6d+LT;Q6?RlewZ#}2QR)6$FV`>VeE zIsU%-{_jbvSLLa>)y!I`qwM3B9EfX{PME*d;PWZ zzbC1B7lDUKW`5^&a<2Mv(S32LbpD>gi=o={)Q%QU`OTJ-7`1SN|E=db_Y3zirv2#K z`RSB)H{ZidWi|hKGJko?K3or9WRd^rSjxG_tAqF4yZ`oIwD2=~{$TZ=yHDRgnHcwV z%kh%bJMEB+V~}O_^FQnP|Nk%c*O~m<-*o$hwuDWr>!++mN;Cd&=WPnpI+Tq(I+MZL zc!n{K{neqfXKlPrIx4d>$sYK0LcD(q&&kCe&IV?d`B^in62IRq=Z`#X_ggK=Zs(UT zeESY>cYjvFlD?y{*ZuylG_8G{-NDm#wKDJf|M&fRtNZi(Gi>`8GS53MTOM;I?0(5{ z^)q&sfBXOc=nn{aXnAMhozm|`JEzR~yy1(Q>C;c`rB|8k6#M62p7Q=@cdUEALovg@ zY0-Hne+CuCFNKa;nFpVEzuyU~*K_)0qAK_CoxJd=-Xg!plBc-NKivE< zp=RxrWq;P#MgBI={+SoNcb;c|6nIo(&0F45kLbD+@xNkLf=3w3Z8!Z$*WDj?{=wxJ z-}j$B`>XQ1-S@lY@gg$3m(TympSFDQ$MS1ElE!nEPWTfm{Eu&%^wang5vhu&%jev; z{g0^?#0&X*FF@q+qdQZ@edx0K*I;-shpoTTD^E>-DEvMCg|3#)<1HUPbwR>sD1N?4yG>sEY4768R`VOGo^1MFDzM!o zRr7vu?2fk!<>U&VX|hZz5A&b)-ja{s>cxT;`)fDPTD^YXE`z;JM)Rj1o_ovr=G1xe z)o(UliY=D2e9_SI`q@g`nxwi*p6aF2nRB4?E6&O?FFtU@eV_RJ{uh~3dAqM$T1@td zw3+<7@}~0TnwMLdv!@14`?MJ{h~I9O`jO|i%KYrwxH;x$f9qyfPcHQct~07H+8OkB zb5+_|`BxLJ*=PT>{aA8W>hHhunO|$}?Va!bp-&#wqLIVw;pbvo$^boGsUa!d|buD)>dD>bGv1xn-@Jf zk#|2jv+w1hOaJT@ax3>)srVd!Ag?L=a};UYWNy_2sOTdi8pZoQuQx(af}bXYKz zg>A>CpTe`wXx2Yk^5+ZJsrfgw?bBC%ay*}x`sDqUwSP>wPxZeGty_0D@Al<&GmO*w zKB-?g#r|tiY|PU?EBv0iD{9(1?T9wad?i-4>+7}Xsry@ke@buEs8^B+6SsY?zdW># zi%BH^wCs&^qSj5>zWnuPuvweSo$Z0_4D-k z+VTzG`894;KA$^PK6l<78Qt$kjQ<~QjM-juU~lx!MAvJdUVr2~?)CZ3Y5gYO%s|PhuzUXaaT*}YB;<>W6l)d8e7Nz=?t3Wf9?u`%rbj&i@{iSl=pN{a{olmFDvMOEWXaDz$oxL49 zcnJ4}3?sYi=gZeM><{fLmb$B?`}fXsyXCX8|>*IH}vDH4CL()6QvPa$xp; z7F(&6Gv@1u{`7u7b$*wERphDr8EaQY+JByT{>8#}xk>KY6Cv|P8Cx3M_dHwf9Psa1 z^(p`Ni~sQ6lB(r6^}6xEEYs{otAjrUTWi>t2J?F*g`+03_pDN455B5S`1!MF>VCml z`~JSo*RDVQ_Jz>#Fo~s8_DAn1yO}EfH2>bRKYv#FYdJB$AC z`oGriuoXleAZn}=%7|#V{`I?WSgh zF1D6vSq6%vzZwn=+vPu6nxJ%q8yJ|@su`58Gs&E^rEKk%g2b3=rafz~bzc!@WVz+p zz>shKRV^2FCQra&!F1N5oj(}3#XOYVj+|zxsFs}jNYH_yQHzPCV)cU!oTx#|z{qln zdrgccudGzQ@`31k9a}jR8XV#|1RlKQPVH+%sbLiw9AJSo#!i6XdU^k8FV;rL+2Fi}asfuWI& zL*PL!gOV%y!kY`3Osm7ao!gGAoo5>@WxVV%ONKrZhr$IJkOhn7SFoZ7hpwtY{XUbd zkeH{LmzQV;N3Hr5A_fXRX-1|!xtFJ!u)tzEBfjxy_fqTkQ|8sj-AZK7d=rs*TKD$$ zy<2!d%F`M{S1&1iyrlB;lBcI1&#F>+e`4aLGczxJcz8CZ`c)eP6UUbM42<`_y$+6Q zf%$98)B~~Qp2F&0>-Su`x~A7`Qr@1xSq2NUWBsRvt?T(7m%n%3ltcTF0`$vX2F82) z%%684#57rLipsa&T{5?9(zo0GyS~T0oap{&?)KY$Q9A;pEDTZ;LpQaa?pp4@ep06| zLKlb0LI%crYHvhV!5yO@z+LfVUGwhebFZrw{rh<-GX3cK_dDIC?Du6y@AYP>cs(`Q z9+WM?!O{3nOyR)o`?d!@B1%9lt{>~({C*d`=jT!JrMcTCURyV{?yGj@`@P-UuYc}h z$hXaU(+u;u!UY~sd@K}S!3>WNhb+eX7t$a1E?v8Q+PvB@DFcPC>qL!(9<22^%!aU7 z4OGTGk#+eBEz}qoS-M!)?e??!8uIer@1++PcW=M?yNe-zZ_%4}gl+ekSSo(+X1B0I z*rwGO~ z;-Yn5h1EQ+$Lr_r%Dld-YPyu@gRd)+kjH%HDjMX?dE#EL#Pusm&|To3G_&s?&f=Og2p6QctkO>{Bmhv$e(Xe zEC3g4bmOX6)!J(|>Ds!fVQVH_i`|~SRQSQ$newn2mVt@m3n-WT+x+|lA`VJHx_Hed zxb<`_&swp0#-Tix3T_b{M9$v=D(b{u8?1tt7$uVqWUoKEHn#rnqpa55Qz;}!L3Qg3wGZ3(d*;{g-cq!*Fu@aL|H0cCnc*RJf|k0aNst}TnABv$-B7sd~7pG>-q9za_sGE zk?vBaUcY`_U;6R!?rYymIvMi)S#QdN3Xe4hHl<$r^77J=&Py8-dv#}D0+n!O!U_j& z=h=Ey!^~?m6RP;^`||q!X=3`*w%u^Lm8pKx@t)gq?R6dAHeau}P4}q;2U;Y{G@Xfh z`)%rUVztiBRPL3VyOiYyaI$iIFiTt>pbGaVIR7)USotbU?~dM+&}fq6vd5=k`{jGATf`YzUI{83 zi0+sA4KWu~{zF*~>p6Zro8{bo%vJcoThSygSB6G2kO$Lby^3KTgc>c7k;>?Q(I6w> z4inR!wRa|8k!580<<-EDAMaT#gYFK26`GLxakJoqwI^LcC4~b^MWry*T6m&uaJUSr zF=nk>pv0kYK?zjbE?K-n4PG8V70y*~XxJ`XHPcL~!9f}1_n^s}(4FV7z@3SuqWZ}3 za~uj68ktxstT$ZFG(tD&td>K=_Q>kCD#Zqe&0GQxa<9!>7zd6YB%)*i1LHluHzl)r z7?^B@6b{^0wkyME+%>v^v@cOhRcvrj2Q?{L;+FLy83Oj|mJZNL)t4HpIvJRJK{b5l zuBI-q3^H-x)zwQA6gM7zbMw*+!x>jnY>9H#IlR+x|L?X6&Ntyf2GWU<# zXCwB(89LSG`_(+FuWz3b3_#XT?02mAypRe@)!Ujmb+RH?MlWa;s^~F`L?| zwR3XZLD?Ie(KEyunfAQBGx#XKnAKe^A+ZP#t5Q zSCbLu& z{Z-876POTOlw`uimAXXxe{r=OW{ZH?yq-?OqGeN~?mQ1|=6!Ann0 zKHaW*z3R)mHIXMxbwFtkoZ#nbZutMN=fT@E^H5sTA%YQGwdT!BoL{?S_xp8|3?DB8 z)o&j+q#%;d3wK7h58ssS+7bEng%6|uiSpWUwPsf(i$2YuUOYby>l|&m*n2_jT+g{kn@c z)C3>ou2EkKZ=9?J6ZcfErafnGIA7spWVys8@W5Bw)ecjwG^77jgMx%RGNA6nS`(JWJ`sfj zwgK4gd zD7rxo3pRo}W3N3kITS9`gPOkQ=PARJJS15tG&tCE3OsnLUcuMRz{Cq`RTl(tfuaX$ z5u}u8OasZhv32{ri{?k z=)OM`g`8|Ol#$0O<|z6$d~al5qTDy(-#7D@e?BiYK0k$xN28b_3=tYBLJ9|FcXppa za!|>(hWA$Qw?NuwG1>=Szm)iP)BO`~@4Fq&Rrh~C`ujs18s-8H3!a18GkG4I$ol6T zI5qXsg@uKiSF={c*W1q9chKC}>i3(9uC~6ziw?q6I_QJOB3xs-5T&5Qf>>4?*TV}O z3tz8dt?;*t-1RFf^ZVVU>hmVB@#TaI-fD%mwK%@?Gcex!SNMDd%rsCyZfRJa#c}YK z7vnv4`6%3k1;ozm%a?eD*c>nhg7Q8P!zBO$3N?llk<3bLoAG+3k1wTxrfq_X^ zVS?I$-EX&D`h4E|^;-40mCH&FabAl5uQj)7)su^R{kGpxs(qsTv8euk_5Z(<A3*(U18w%b|nU(P--x7hD(%@b(p*7!{@qkZ-Pp$ED1 zzboq^4O%TY&yvMj5nttbE2(>j)${rKwI92WhFU)nc=_DE|N7oPo3@-(J9^tY{_G4z z%ZwO*DH~AF@7p}Nz1G=bD8}EYuXp>f%|HK=5K@Q=c<2UbANc)j_N6xIV~g|le=~jA zW8Brf??KbL1@DFZJ*wXX?wK2o@bF6(iG#Dxb;ck>I20xbtPp)rdad@`#AjzF_Swh% zdRH!N{apP0Zht9zySV3-mFBnCR9LS^WR6Rm0uOwjyC52Sj4T2k97{Mq%qyAnW$t^A z>bI)4|0c5LS*4%eE`2Ga3y(_L;5bJ2PnB(_gP2=PJ8-^xJ9W?%DXnvi$zq`&9iefdFA7#*j5$ZOuaOHpUT=; z>H9u!zwe*kXW_K#lUL@42k(AHf(BPgW;LuYnE@^1)?AMAuNHJzupE@TwNCs)3Ic_I z;0If;FS-B6P0GaMR&6=cZqtJ&8W@;l6*r`Gmi&DBG5O`CAIGLFb}x6ioTWXyn;MTIOC(#O0UIgrN`5vGD!XH_0Z2UOJ+ zsWa`_$oyZa!NH7?Y0u3)C$MbbGx0ePY|0}eA@1I=J#mdli-S9;LR|4~*+#HQ$ixL7 zP%*L6d<8etT1|(B?WVtbqcAH}X;2Td*>4Fa$Cu3vjQ8r^G4O%v3s4n?Fl?^2L&Nss z-}9y^U(jG=+S5B}31*+xL18fiK*8R?)J z_0H}sycnuxMKv(w2bc3tbzY##!tr67i+P|Ya_R+_D;e$}{Wn%`;cVO|s&L@;_Iomz zCJAV8Kq}|`bHF338#-B81RN0U5e%!P^r0iG@3v%Knxffj&!(!>;7|=3a(@4BX(=K| zI20}@G>XSfkgITr+TwAmuKMNO_tUoBURU%mf1jr9-z{H09KKZb`shz*f#sj#Vx!;f z39fo$_;O0{5?}MTQRju2IJUGOSQUC{N8zPCm6tv}Eu3vLsr-pY<};tnA0NJDsIPTk zXvESWU}T8`cNEXpg+*;#v|5qxNX<*{OOffNYq`Dk_p12LP_X?z=gYnCUf=I1>qaWo zz6!o%zVBkGW!Z`NUw@9t1)HUD49o;UT=^`?HirOkKl+)U-%-MKOnc25K{ zuj#zp@p#Gm+V$G?FDIS;S}0d*a?hz@dq~H1XxL2$GZxL)?`nbMlV*557W9nNTEdWkKGN=kKU=P&U|<0hb3q6$D^@~ zchcv-oz-8y`RtPE@zdVzi?;f3;ECL%=A-;}8nR_GzC4x}&;P8=UT3@g{nc=9^*I9h z>IZJ$?KFk@=7Oq2!}ilns6$N)JPs_*sQ>@_Q}F{UiPwTVnc0tBHpS7u<&9C-~C_Hwmj~;)Wmx9c3idZt>pfSs`l*lQ?JF_-}7tOexnmsEd+pO zA#OCoxZqBBgM$|1eT%coCAG#|)aUn<&$N%1GS!OP_mJ&nvb?t5tJm?<&siw*aak|e z{ccz6ZGS7zU61)*zFhwJD@(=hoCs+5xA7DcONI9B4Xh{uwL+udqU$B`eG0WdPhYZp ze&l9ceQed~=p`|F`##RfKQPVb`CRem->v$;eJ$bhv;RHw%R_!|KB@Qn_dhA(_;AeZ z4{TKPDl5l_uc0LTA9!1s@hj`wey<2IAv5-Vw3*gLo7`t>w$E_(b+h7z@Be|+uJ2Rhl~TD|vDhRxY7Qvh z-P==fQ@Q;ApY)d(-ItbL>%4wd`ho423l1} z_N=?K-#E=9`&!7YsM^2av4N|fD^|O|KQ}AXS{~Xu&tPNpuRX$XFYN{AD@Y1~*4Z2; zjSb6#Z=}m-uHU<4rtwln=2hXRf=%}PTn0*fC)JN%Hm;b@%DqIkY{IwO{#7p)zDjLb zYyY$JiQq0+L(xG&^8m-aC2Xgm2_Is<11sbC?bQbuKWaqX7seGXLJ8w;N_lnL{{rgD#WAF7zw@>Q|KhRx%$O@XwG(bJ` zIb5%i5}gA>nhDN5P7coQEZi-3kc^8EBqn>DTsjo(0B zUb|XPZ)6uQ;ADBX%UjCAU|#en40WB@q&zOOjA%ms&E&UD4 zGbBtj^7t9=eR~-m)#LygeY{$7qYb74l7&I3c*}%_qpHIL7y<+hI4Tvi~l( ztQPDmWa0v7zSw%%>P!tLRz{{htX!y@n88W85iz>D!o2})1`%W(*kZ&~rh`5x=2H{8 zkP|-GB;*N%SsOFx%5R^`;KqFA@HE((xyg_5`C-W9xLUuI7PKRAm+xF7@_kq?~<9+6jF<*O>_uz(%ZzHd0b z4yp>`NCu{(Ox@z%<@Ypg|L@synD5e!jhCiHi#>lR%EVz32dXpfSY(3x7&~hEpU7GU z&9V%ZvfK0V_!1Wufeb!Irafmb@?Sx1b!JpH?tV8(&2LIv<?$`7P% zbN6iA^7Y!K&F8(C*`}O4q6M|>#Y6_idvYH)8l%|O;Go+0dfj9{+t6MA>oUK*@Y)#x z%h|S|(n8qw8)EDol-dMds7*NSX?{mRwt7v;kB66f&5f$IOABd39FQ@Af$<)jV-OFj zB@HFN-+o+i*R^|>iJH%Zb^GIPm0Vs|`giY>z@v2&z&fL%!9(zVek>K%JHKaIqv*eo z&=?*+bzX&2)eFU!m%NuwkMH9z`1SHq>hz`7H^2k(^YrpF*Fol-V{iZcdNbYo`JBtO zw%HJ~y${r%U+$N?*2;X>HkJx$T`$O-BO<~lHZ<>keBArIjq=+ai?<{=TFsC&USgVk z>8b8n(~1oyl|?Uay<2ZvReo*FyVu}8tHOmSrr+<+bKLvS&cN{hzk>BDT?PiN=b$>| z?cM4%v8glJAH;L#YKCo5YTB+~y!B|XL8z7{)6$|&7A2M@US7o`8LwEl`&@+gBrvA4 zEH+?rk*`S{P(lqn@z9 zhqB1Sxzmx=uNAqmQF*3;;!GRMpJ(%D8c-Ml>fCw=(t`n>Hc=Sn!@}#xVlezVM zo}Kl+_i38ZmA!vHefquidYS3Ovk(U^Ki=$%oEz44{uaR{`sRaD#c=Y4r)1xeVbe+m2nVLle zKJW>hGC>XzaO7|266w#lk-xv_JJY?qzuFt#mx&5|*p{7s1ugk-2xJI1J__Vr>^}V{ zN5$w&k8O7t}aKKtaZ8AqsM?*uo%Z=D6$ofHc2FNlp=9gA2z4T64;6qvb z`L{FBRT!u+GQOX%!$<8x5F=y0_s-ZJbQJ-%jEwnTtCl@_Cn)ft%zyn`3|kCz7#ZJh zxG~3NfhrTzy=ex!IY6!h1s%++1-qEm@0rvi8+~nmPNY$?Tj}>jIMNO{9uG3D4p$Fe zKCN3XVnR$wV6%RLB@;`VuEGOZv$?`30kt5S<=uY!$bB|5&GjvRJjlDbwN}VuTdsH7 zS+B6QUS_#ku3grZ)m4*RyC->0)(})MxWvSCZ{oX#*C@%sA)5Qo$NrOBGEaKztM31s zTXZEbH+J14uv<>Eu;39}o=!P4qv-Qt-iq^9tM~jaw`h5MT)q7M-uK4q_ZaPbwW{{KaWexW zQ>duMzwdYc{0v>Y&$N`gBKq+gEl^Cropa$C)Bm6IPu|&i()PQGyWGv_9KjD_zPlH@tK|OYGYSq3 z4xy}F?oW3#mp^#NS;79T5H%m{a(Qq>c(S{k$1AQB$01VweqGGX z{Zk(4-MRl~naq3}%b5m>J0G?^X%YT)-qZ;k8d16m4|2Eec#RT!4og|ue)HV($_?5G z@@YUMBV&GQ^zs6*bufQtv^F$s7d^EOT&5lt75Kmxey$Rv9Hv0wf+!Qyy>%IzeL>N< ziHYgnscn7_QB`m?u8p3|%m$s}#A#x1V4IbFS7?$*!WTb-HC z3JH8L^S-tVMLXE0zbx0*>({T?@MGg4Ug=4q+LK;h{_Sh*&~QM3m1R$Ji%~C1?1PfQ z0=vek+LJ$@_b+?rb8VgOb)gT}(qH9(icw}v<8=i!d1 z#dI&MW4D_+H%Eo_&fVYKU}k}R$HF0y;S4SETRAIkZ*cTlHKn0pd(hh;9Z+^cDFc5o zG2MIh{Q}=DcZUPF?e+$6VyMstl{FXBX6H<7XxQ$0LC+c8GKByxj*7pk0vU4}8nzcs z3xhNtkb^~{rJRk8^QZ1!UWOik4y>$km8z1~+jMGf zZb~`U!oT$FNY zNm|_9TP#+=4hL?}tlg@O>dp<(51i$_rs=%wf4J`CKQ%WF0fhipM#g;6-OJKYRT^+K zw#)maoqJ>dF*o1)+S=E(Y+3&$K^hzuXfZL}TUFvU7d4a;I+&KvohDNfP;}Y1)?dcn zZl=Y?_mGxbfG;CszUl8}S!i}FIIy^1?YR8>KfkV@tiEskAcP%M(l0o~#B`6#-2Ww- z`T0!C=g<4|W%t5g4Gr68)driQM$V1Z2WFXieq_$B{}ue{g7eeobFELm3|_9b|5s?y zgT`vxyI|!jCNwl`cRCxS1?GX02$&(!cYv8+MLcGL&F`GbZ>rN(%UaSVz1=qf!=Ya8f75TpzpSt5!l8>p(HdB3k?CtJ9o6mV| z%enb|Pn7Yoe*bs-c3+#gFF8DJW=v`5eZ%`40x$R)AJzYRe*S&mUbF1!O!v}GMxnIA zSva=LJaE}xeX`or#<$7+s>fwEzh}zdZ+mlN^8L?QXK&R1PdPUy%KcZLX!(1+*jwdi z%$ZnrIX_r3`Od-U{N=2BzH(X32WKRB?hf&Iu;=q6>-Rd$S7#b)hwKOf#~Ua>)w&%) zfeZJN!iAYkpP#92zcVT4#sr@k7iLy(0j0^3$qfzLXYa~gj}pfhzA_z`4HwdhzV=mw zfkQyyf*=#qy}pdazHp1djK;Xe-(O1`z%1MBD1$!1OE$w zDzB-mEPF!V&ASD18X`wI@HU$7x2$~F`sD6=KfOH$mcL#++4p~CL zNfG3M6-yc#wol3pHbIS{8@msfURMc^pPO?x>fB5VL(7jxo~(}dUFLW9cV(ZziB+6!lW%UG?9Bi6$=qAf+c`NZcCUT>%^Tqa21ce>#SOL9pN?vHe0*xjbZ_C!nRov&G2Q$2ej*oYaDw!_63~#lxZ!@?QzoW+Rtj@( zO=DuZw=Ji85oS3nr5a$fD)GyUB2&iv#UX1i>oGCiD|^{ug{g9ulEQ;ruGOx$ah0GA z4#g}id(IYa&8p_$sJOkDe-$FhK?-#afs9EF4ck38Ov~kFWXzA=;wXrc*ubh__4r!_ zg$KFqdMg@19zt=>f=ni+d&@3H?cil$*^_&LMHQR`Q3BV*-{HV*n`rI+Jq-=zPpTI% zff5jka~jpSI4Wj$WWN*?7x>Uu!7T!|9K|^od>X^!JdgG1uHOLeF&M8;%VkuU@^uz7EjRyHz4vWf%H3V}(l>saZ@6RP^?&z&ckk}p%OQZHSHr-_R4saA>uQ@V zkPdNWJmdS!<*MwmQ|^=m-%mZ>@2}?f#@cL7KGVI2D=ga?8q+{S0;<`uXz{fmk?Cl% zU{}Mo3SN;9cE2);{(L;SzJBhW?|YyAar+;B;C4pfYe9v8{)UF_a<FLRm#=pcrclxW`|C9D<*Y&UVW^)7?-xpk&puz>^WLo)h!?K`||ur?Q&Z7;^qIu)~^hn9RF9#^232o)wcgdKU`b9xp{ z({#5!UmF+O_x;NiySP6RckY&7x7+so#nf=M@Yt=-JI=FM`7iikJ>P!+=Wkd2tc5;o z3pT$3=`J8#d__be??AN9w=b7ZF7ZrF_~VcFBWd5aKw zOshbP4RL*^Taw1V;$M6`KDqSuq|N93x}HDS_wCoY@=4wLe$({*m-$xbmi+k7RneaR z+=q!p3N+fXXQe25;$!1p6F1jqZj|1x1jkQNCw4tgz4zzYlUDIb<@>F7emL~x#p21o z-+Ax*|M$u5`|E1Y&3n%O?Q{0^xE%!>(r36gFf={`<*D@SRP;RM*O+^IQsm}Ip{pnP z&OUV%X_TT$zu>It$*bY&(c7mb3#+TC&zn*6b@kqd1@GTHtQY#gcYodFb9Z(aYC8Rc zRI5wo-t@l8!m=my^^7WX*YD-Zxj)b5qsvd{ELRSjtSfQL`F52=j~WnBW{0ddy=bBZB?uJk}s#YIV!a8J4vII+2DZAXhj*23m5sY&FP*gn~=Z< zz9o+(z)DawHU1aBv32#I4UhsT^1$o!$xB#S_9)+FypCCaubATSanAFPk55lx-7_^S z)?G~d!P^ttZi!$R>#&)n%iZM7_V}H<8n#ECx+c4(p`m<3&@;@gj=+!E4j)&vU0m$V z#kk(SbE25QhrSZ_81%HP;IJTwiRqr#*42-@ML(3)bmvY$t1}cF7EEMfy62V^wK|TG zF@JB72d5J6@Pg#D^3B4riO-cj)-Mu zu?`1rmn4a!ge5o_UZ^lK=Bt-({#7mbA#8tp8F~g2P`Kd4$e7Q)F*a*H6Vts~p_!N} z92SH>m>NE5WAf?ur>ChJr+tz8-085Im1WObv2xV0D3JY)UeFwnmw2%0$)awpy|PEc zV;l}de`M1_562fijnC)%zuRwL`D*2p`TsO$S{hfrjf8Xmt&By)!ww7DxS2%xamRCui zbJ6;lFe4MF_sWJS@}S|Fs}gW$XZ4qhpl(s8x*EUT44YpWkq;iWt7l(VTkO6)@SvBX zLxKw^6`p$k2(4^pXxztibnb4~9lxip|I+vW_I}myn83RDJrH!&UG;3{6WRA3W$%3IzOH`4EdQeV4yXeW?r`At-o4w4 z(Y>g|bds?+0U=|oQQ znICsQ3RHH2>b<#F{dX&(#kQ44!u54e_W!s3dE@w`&*!}BJ}{TRFAJ#7ZFO*FVcEl* zxCL{lb|cgGJHhw1?#_#Szbk#t#6Q>TRoUh5Sesp%-?h7|dwX~}*O^&UMe^5A6c+f< zH)pRldPuN*)B149`{d((^?sX4G4-{d4sk!dQvdh*6Z?P3>z9AEEC5$=wQX#FQ~&&U z@^1Hc+rOu@C$C<2>${DyRq%or|Ns1)yL0=hH(m##r~0j&jOPCrQH@G>Y=1s^V!q$= z-_QAzO|#kD3xB^o*(cjP`*;85Z};tI8b6nJ|6Q*S8oAr@f5*gk0^oYAuVgtZ_JW7U zDdFRzCs%?`cS2@;nk`O>FtJ$0IUKnCQFj|!|IJ~6FiV%a?w*Idpa$BrvnO{H&fTD0 zt>CZ#SGlmoVL`=}LoJ*qpPYQUSSqT(bGm}V0wz$O-o&#LE&Vxs=h*Y7Xyeq+xINxDb?S<=TF_+ zt+`=aiZ}!$0vQ?KcYNAy0WL^DSsTh&5T5Y+TgtJXC!Ok5$B*@RF7xwUw_6XC;tS_n z1t&Z{mh$jW%C|R9E;w5k_m^FbU}Vfswq1)c5Vpb_+pN#499S<-+Dt_|Oy%ZE4p{4)T4&>V;P%wq?U`t~r-7ld zOz_6m)oW7D&Pw_8#=@)bFBw&o7&7X zx|!~=IYzaiB{EQd*_P?3vedWR>th@bME9qZ&((hLc3s_VU1*pgvZcdpE=coK$M3-H zWxiR5H?i!=y&f_UL3Oar z8TyQj`PE++YhC52Sba&e6uq_0!cnrNp<(;t+qbuz<;0lZgBFp^;sPJ|?wvXVO0nx( z{n64OC@kE391d6~eVMv%L&Nq8i8MzlTM#OppH+RCH8y zZ}QigTf4T}{BB@=6?t?g%PMV!2W#`*hxj349~4*|0vXYbs@{`={XE@zz5FZ|)@*yz z%&)e7U(~dk&u32_6~7+;FT39`V#kJmTh}-&=wf2Jr&Z>A7TNcpplvi0pJBO~r~Lg6 zWy>cMHr;TS_qF@;@X6BYlbZSc=GE`}=gU#?+I)-CR#uihSKrNlh2e=@^#B|9Z58+4 zotpY|f3?_;t?N~%$H(>A{7m_K{;}wXcPn2Pn{#ng?B4kLn-7K$7J+7d<^N5PDGmAG zy(VI!Om&V?RmrAA=TE<`f1htWMUe4*z?L&_T|mu`va9pJSry`WklPuVa5+WpS!oMrllXzSN&Qch31^*MaK-tp?HEO~i*yPdz^-8%kFX78s{ z$L<8~c|L3Ihn)p8?e`0pUD{bK@S!X(|1HMM=mj&T)nTf~WvAEtx_mNxpV5|0d7x3m z#osKq-;3%yKR5dKkAJ__i=SPw+x9%8@cX+M-a8ZW4NdbbgdK!azUJEmcjLnW&ZJUi?8@o|Q~ZD`wWVfAT3x3{ZHo2iuFnHW=FTUm2v#zSpR zP|xbjIRSwWWp46so6u9h1v94EW~%A)W{Q-*;QKPmTz}`wWsPZXe$Qwqzwee-vaO+E z`{v)-e)L@g$H+o52s`3XJl#?h}dIc`SnW5tu0SF)!7OU*)8L!_|Ev` zLmX&!cWEJpa~&ESa#8dCw>b-YxGLHgrM{cT#B|Trv<4%&2`Fsf2lbbg`+d6VyWK5( zW^CveW-E+1V9lWtoW;dak-bR*v$9tRh<+e_-^1?vPDmSSrS0--FCx=bqqkkLnzRNw zz+zWdQ}q4r)_%#id5nzt<=^kW!zfo-+7vcqg%;f@Jh|L{Zp~9|XN!;d!QlAe? zFGMxwZLdD1d6JP`<@>$aJ<`{!*9gvM+0z^F=A0?Zp1U{8cVmW_#1e;(ckVoj-2M8x zdbjxZ|AAZE<@MgbS*!j<^aI}}$z@-K1wMRJHbgFCK{cy@g2Mul8+$A@>tx=pTz;yc z{Lzu3n&0124!1q&v!3*Ny?%ecq;bvGD&~8j{O2apKjFq_(<>is&wt%11uD9~cuvL4 zB)VcZwrE;@KC@(pbNPbndn!-*n)8|@-#`9}vqHXpPmNwd;nYW;f8Ps_%U|*M_qt0x zWg-F}wk3yOaluGD29KDg$64)oux*F;H~IH-WV1egW!ba;cWlv>z^9_&>;0DbXxj6a z)~qm{eP0eVK>wmy3p2y*iU4K)v*HUs6~y`5McOTT_~G#HdA3(R?o0n2#tNGL+qJ|E zOHyayxFR6%;hKXXrn?*(7#g)WI4V|OP&LI8g9;7{+8P?Rvvq;zQXmy3QZR)$Ivj{@ z)mu3cy^I9MjwuVv9@79G%=RQGuDC#ri4xCVB-0R$s}+h057ss&bzn6%R6yWE*h23t zEcwQv!GV{NF`rk&8N0Df4Gr75VuO^S2?jNmvN$*@vKL5T&P;ebodKs57F0j$QZQc!rX)@k|`wAo-t>^6eZo|Hx#goma-fRZQ;(Rxt(^Xc@>85b9= ziQ1~QI(+@Tzjbdm9{;q=eDBJ<*Jo#&KR>BH|IOxe)%QNn>yfkFbw1+b(jDK{#qR#q zv;Xtl_h%g0Wv`^KxfwsN`rXbug~w%o-hE%U{rBJJ_Wvul`>vfeJJh@I?EZgW*Eh%g z>&>r8R`Z*4V~ZW%E&4{?2$Itdp@r^&+hll^UwaiudkoK^X$J* zXN=E3IcNR;%=Ud>*WP7(zjFD!S@ZwCNk6+UXa3ja{`2?k7O!Pv|6x@A?oQ3$*YWFD z_cQC*KWcw=X69tj-A0zLRxEbgBgTqIqOiPBvJd37$!fl*K=(M*e40G}jPCY3Mj9V0 zKRt-=EIrG=S~B;ccxkHErakvDYTRaDiscL{*gbu#AVMp#V~o$^@^UZelF{r G5}E+aMQacM literal 0 HcmV?d00001 diff --git a/doc/internals/centered_continued_fraction.qbk b/doc/internals/centered_continued_fraction.qbk new file mode 100644 index 000000000..0a012cce4 --- /dev/null +++ b/doc/internals/centered_continued_fraction.qbk @@ -0,0 +1,55 @@ +[/ + Copyright Nick Thompson, 2020 + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + +[section:centered_continued_fraction Centered Continued Fractions] + + #include + namespace boost::math::tools { + + template + class centered_continued_fraction { + public: + centered_continued_fraction(Real x); + + Real khinchin_geometric_mean() const; + + template + friend std::ostream& operator<<(std::ostream& out, centered_continued_fraction& scf); + }; + } + + +The `centered_continued_fraction` class provided by Boost affords the ability to convert a floating point number into a centered continued fraction. +Unlike the simple continued fraction, a centered continued fraction allows partial denominators to be negative. + +Here's a minimal working example: + + using boost::math::constants::pi; + using boost::math::tools::centered_continued_fraction; + auto cfrac = centered_continued_fraction(pi()); + std::cout << "π ≈ " << cfrac << "\n"; + // Prints: + // π ≈ [3; 7, 16, -294, 3, -4, 5, -15, -3, 2, 2] + +This function achieves the same end as the Maple syntax + +[$../images/maple_cfrac.png] + +You should convince yourself that the Maple output and the notation we use are in fact the same thing. + +The class computes partial denominators which simultaneously computing convergents with the modified Lentz's algorithm. +Once a convergent is within a few ulps of the input value, the computation stops. + +Just like simple continued fractions, centered continued fractions admit a "Khinchin constant". +The value of this constant is ~5.454517 (see [@http://jeremiebourdon.free.fr/data/Khintchine.pdf here]). +We can evaluate the "empirical Khinchin constant" for a particular number via + + cfrac.khinchin_geometric_mean(); + +If this is close to 5.454517, then the number is probably uninteresting with respect to its characterization as a centered continued fraction. + +[endsect] diff --git a/doc/math.qbk b/doc/math.qbk index 09d8ad9b8..df0770f6f 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -752,6 +752,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include internals/agm.qbk] [include internals/fraction.qbk] [include internals/simple_continued_fraction.qbk] +[include internals/centered_continued_fraction.qbk] [include internals/recurrence.qbk] [/include internals/rational.qbk] [/moved to tools] [include internals/tuple.qbk] diff --git a/example/centered_continued_fraction.cpp b/example/centered_continued_fraction.cpp new file mode 100644 index 000000000..70663af83 --- /dev/null +++ b/example/centered_continued_fraction.cpp @@ -0,0 +1,46 @@ +// (C) Copyright Nick Thompson 2020. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +#include +#include +#include +#include + +using boost::math::constants::root_two; +using boost::math::constants::phi; +using boost::math::constants::pi; +using boost::math::constants::e; +using boost::math::constants::zeta_three; +using boost::math::tools::centered_continued_fraction; +using boost::multiprecision::mpfr_float; + +int main() +{ + using Real = mpfr_float; + int p = 10000; + mpfr_float::default_precision(p); + auto phi_cfrac = centered_continued_fraction(phi()); + std::cout << "φ ≈ " << phi_cfrac << "\n"; + std::cout << "Khinchin mean: " << std::setprecision(10) << phi_cfrac.khinchin_geometric_mean() << "\n\n\n"; + + auto pi_cfrac = centered_continued_fraction(pi()); + std::cout << "π ≈ " << pi_cfrac << "\n"; + std::cout << "Khinchin mean: " << std::setprecision(10) << pi_cfrac.khinchin_geometric_mean() << "\n\n\n"; + + auto rt_cfrac = centered_continued_fraction(root_two()); + std::cout << "√2 ≈ " << rt_cfrac << "\n"; + std::cout << "Khinchin mean: " << std::setprecision(10) << rt_cfrac.khinchin_geometric_mean() << "\n\n\n"; + + auto e_cfrac = centered_continued_fraction(e()); + std::cout << "e ≈ " << e_cfrac << "\n"; + std::cout << "Khinchin mean: " << std::setprecision(10) << e_cfrac.khinchin_geometric_mean() << "\n\n\n"; + + auto z_cfrac = centered_continued_fraction(zeta_three()); + std::cout << "ζ(3) ≈ " << z_cfrac << "\n"; + std::cout << "Khinchin mean: " << std::setprecision(10) << z_cfrac.khinchin_geometric_mean() << "\n\n\n"; + + + // http://jeremiebourdon.free.fr/data/Khintchine.pdf + std::cout << "The expected Khinchin mean for a random centered continued fraction is 5.45451724454\n"; +} diff --git a/include/boost/math/tools/centered_continued_fraction.hpp b/include/boost/math/tools/centered_continued_fraction.hpp new file mode 100644 index 000000000..2016ffd20 --- /dev/null +++ b/include/boost/math/tools/centered_continued_fraction.hpp @@ -0,0 +1,154 @@ +// (C) Copyright Nick Thompson 2020. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_TOOLS_CENTERED_CONTINUED_FRACTION_HPP +#define BOOST_MATH_TOOLS_CENTERED_CONTINUED_FRACTION_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace boost::math::tools { + +template +class centered_continued_fraction { +public: + centered_continued_fraction(Real x) : x_{x} { + static_assert(std::is_integral_v && std::is_signed_v, + "Centered continued fractions require signed integer types."); + using std::round; + using std::abs; + using std::sqrt; + using std::isfinite; + if (!isfinite(x)) + { + throw std::domain_error("Cannot convert non-finites into continued fractions."); + } + b_.reserve(50); + Real bj = round(x); + b_.push_back(static_cast(bj)); + if (bj == x) + { + b_.shrink_to_fit(); + return; + } + x = 1/(x-bj); + Real f = bj; + if (bj == 0) + { + f = 16*std::numeric_limits::min(); + } + Real C = f; + Real D = 0; + int i = 0; + while (abs(f - x_) >= (1 + i++)*std::numeric_limits::epsilon()*abs(x_)) + { + bj = round(x); + b_.push_back(static_cast(bj)); + x = 1/(x-bj); + D += bj; + if (D == 0) { + D = 16*std::numeric_limits::min(); + } + C = bj + 1/C; + if (C==0) + { + C = 16*std::numeric_limits::min(); + } + D = 1/D; + f *= (C*D); + } + // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1]. + if (b_.size() > 2 && b_.back() == 1) + { + b_[b_.size() - 2] += 1; + b_.resize(b_.size() - 1); + } + b_.shrink_to_fit(); + + for (size_t i = 1; i < b_.size(); ++i) + { + if (b_[i] == 0) { + std::ostringstream oss; + oss << "Found a zero partial denominator: b[" << i << "] = " << b_[i] << "." + << " This means the integer type '" << boost::core::demangle(typeid(Z).name()) + << "' has overflowed and you need to use a wider type," + << " or there is a bug."; + throw std::overflow_error(oss.str()); + } + } + } + + Real khinchin_geometric_mean() const { + if (b_.size() == 1) + { + return std::numeric_limits::quiet_NaN(); + } + using std::log; + using std::exp; + using std::abs; + const std::array logs{std::numeric_limits::quiet_NaN(), Real(0), log(static_cast(2)), log(static_cast(3)), log(static_cast(4)), log(static_cast(5)), log(static_cast(6))}; + Real log_prod = 0; + for (size_t i = 1; i < b_.size(); ++i) + { + if (abs(b_[i]) < static_cast(logs.size())) + { + log_prod += logs[abs(b_[i])]; + } + else + { + log_prod += log(static_cast(abs(b_[i]))); + } + } + log_prod /= (b_.size()-1); + return exp(log_prod); + } + + const std::vector& partial_denominators() const { + return b_; + } + + template + friend std::ostream& operator<<(std::ostream& out, centered_continued_fraction& ccf); + +private: + const Real x_; + std::vector b_; +}; + + +template +std::ostream& operator<<(std::ostream& out, centered_continued_fraction& scf) { + constexpr const int p = std::numeric_limits::max_digits10; + if constexpr (p == 2147483647) + { + out << std::setprecision(scf.x_.backend().precision()); + } + else + { + out << std::setprecision(p); + } + + out << "[" << scf.b_.front(); + if (scf.b_.size() > 1) + { + out << "; "; + for (size_t i = 1; i < scf.b_.size() -1; ++i) + { + out << scf.b_[i] << ", "; + } + out << scf.b_.back(); + } + out << "]"; + return out; +} + + +} +#endif diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp index 4663e5952..e01553a9f 100644 --- a/include/boost/math/tools/simple_continued_fraction.hpp +++ b/include/boost/math/tools/simple_continued_fraction.hpp @@ -140,7 +140,7 @@ std::ostream& operator<<(std::ostream& out, simple_continued_fraction& out << std::setprecision(p); } - out << scf.x_ << " ≈ [" << scf.b_.front(); + out << "[" << scf.b_.front(); if (scf.b_.size() > 1) { out << "; "; diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 03c2fc7c5..4833b6de3 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -895,6 +895,7 @@ test-suite misc : ] [ run test_constants.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run simple_continued_fraction_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] + [ run centered_continued_fraction_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run test_classify.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_error_handling.cpp ../../test/build//boost_unit_test_framework ] [ run legendre_stieltjes_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] diff --git a/test/centered_continued_fraction_test.cpp b/test/centered_continued_fraction_test.cpp new file mode 100644 index 000000000..e9a0e9a9d --- /dev/null +++ b/test/centered_continued_fraction_test.cpp @@ -0,0 +1,148 @@ +/* + * Copyright Nick Thompson, 2020 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#include "math_unit_test.hpp" +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif +#include + +using boost::math::tools::centered_continued_fraction; +using boost::multiprecision::cpp_bin_float_100; +using boost::math::constants::pi; + +template +void test_integral() +{ + for (int64_t i = -20; i < 20; ++i) { + Real ii = i; + auto cfrac = centered_continued_fraction(ii); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(1), a.size()); + CHECK_EQUAL(i, a.front()); + } +} + +template +void test_halves() +{ + for (int64_t i = 0; i < 20; ++i) { + // std::round rounds 0.5 up if i > 0, and rounds down if i < 0. + // Should we care? These representations are not unique anyway; + // In any case, this behavior agrees with Maple. + Real x = i + Real(1)/Real(2); + auto cfrac = centered_continued_fraction(x); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i + 1, a.front()); + CHECK_EQUAL(int64_t(-2), a.back()); + } + + // We'll also test quarters; why not? + for (int64_t i = -20; i < 20; ++i) { + Real x = i + Real(1)/Real(4); + auto cfrac = centered_continued_fraction(x); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i, a.front()); + CHECK_EQUAL(int64_t(4), a.back()); + } + + for (int64_t i = -20; i < 20; ++i) { + Real x = i + Real(1)/Real(8); + auto cfrac = centered_continued_fraction(x); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i, a.front()); + CHECK_EQUAL(int64_t(8), a.back()); + } + + for (int64_t i = -20; i < 20; ++i) { + Real x = i + Real(3)/Real(4); + auto cfrac = centered_continued_fraction(x); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i+1, a.front()); + CHECK_EQUAL(int64_t(-4), a[1]); + } + + for (int64_t i = -20; i < 20; ++i) { + Real x = i + Real(7)/Real(8); + auto cfrac = centered_continued_fraction(x); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(2), a.size()); + CHECK_EQUAL(i + 1, a.front()); + CHECK_EQUAL(int64_t(-8), a[1]); + } +} + +template +void test_simple() +{ + std::cout << "Testing rational numbers on type " << boost::core::demangle(typeid(Real).name()) << "\n"; + { + Real x = Real(649)/200; + // ContinuedFraction[649/200] = [3; 4, 12, 4] + auto cfrac = centered_continued_fraction(x); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(4), a.size()); + CHECK_EQUAL(int64_t(3), a[0]); + CHECK_EQUAL(int64_t(4), a[1]); + CHECK_EQUAL(int64_t(12), a[2]); + CHECK_EQUAL(int64_t(4), a[3]); + } + + { + Real x = Real(415)/Real(93); + // [4; 2, 6, 7]: + auto cfrac = centered_continued_fraction(x); + auto const & a = cfrac.partial_denominators(); + CHECK_EQUAL(size_t(4), a.size()); + CHECK_EQUAL(int64_t(4), a[0]); + CHECK_EQUAL(int64_t(2), a[1]); + CHECK_EQUAL(int64_t(6), a[2]); + CHECK_EQUAL(int64_t(7), a[3]); + } + +} + +template +void test_khinchin() +{ + using std::sqrt; + auto rt_cfrac = centered_continued_fraction(sqrt(static_cast(2))); + Real K0 = rt_cfrac.khinchin_geometric_mean(); + CHECK_ULP_CLOSE(Real(2), K0, 10); +} + + +int main() +{ + test_integral(); + test_integral(); + test_integral(); + test_integral(); + test_halves(); + test_halves(); + test_halves(); + test_halves(); + test_simple(); + test_simple(); + test_simple(); + test_simple(); + test_khinchin(); + #ifdef BOOST_HAS_FLOAT128 + test_integral(); + test_halves(); + test_simple(); + test_khinchin(); + #endif + return boost::math::test::report_errors(); +}