{"id":9470,"date":"2023-07-17T17:57:55","date_gmt":"2023-07-17T21:57:55","guid":{"rendered":"https:\/\/blogs.mathworks.com\/community\/?p=9470"},"modified":"2023-07-17T18:05:01","modified_gmt":"2023-07-17T22:05:01","slug":"continued-fractions-and-computations","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/community\/2023\/07\/17\/continued-fractions-and-computations\/","title":{"rendered":"Continued Fractions and Computations"},"content":{"rendered":"
\n
Today's guest article is by Wesley Hamilton, a STEM Outreach engineer here at MathWorks. See his earlier piece on solving Sudoku puzzles<\/a>.<\/em><\/div>\n
In this blog post we'll explore the fascinating world of continued fractions<\/span>, like $ \\sqrt{2} = 1 + \\frac{1}{2 + \\frac{1}{2 + \\frac{1}{2 + \\cdots}}} $. While these objects have been studied for several centuries, their presence in modern math classes is relatively rare (unless you're studying number theory, probably). One of the historical uses of these objects has been computing rational approximations of irrational numbers, such as the 5th century approximations for \u03c0<\/span><\/a> of $ \\frac{22}{7} $ and $ \\frac{355}{113} $ by Chinese astronomer Zu Chongzhi. Here, we'll focus on something even less prevalent: their computational aspects. The goal is to start exploring statistical properties of continued fractions, generating conjectures, and then pointing to several topics that can serve as computational playgrounds to get one's hands dirty and explore even more.<\/div>\n
Before defining what a continued fraction even is, we'll start by defining a few \"important\" numbers to use as examples throughout. Next we'll define continued fractions and compute some by hand, before introducing the computational algorithm to compute these things. Afterwards we'll explore their statistical properties, before discussing some conjectures and neat facts, and round out the discussion by including links to other interesting computational topics to explore.<\/div>\n
\n
Table of Contents<\/span><\/div>\n
Some Numbers
\n<\/a>Continued Fractions
\n<\/a> Infinite Continued Fractions
\n<\/a> Continued Fraction Notation
\n<\/a> Computing the Coefficients
\n<\/a> Computing 's Continued Fraction
\n<\/a>Statistics of the Coefficients
\n<\/a>Other Topics to Explore<\/a><\/div>\n<\/div>\n

Some Numbers<\/h2>\n
So that we have a rich class of examples, let's set up a few \"important\" numbers to use throughout the rest of this blog. Some of the numbers are important in some area(s) of math, though the whole collection should give us a nice, diverse set of examples to explore with. Also some of the names used are likely not standard, but describe how the numbers are constructed fairly well.<\/div>\n
First, let's set up some parameters that'll determine the precision we'll be using throughout, along with a random number seed to ensure reproducibility no matter which machine (and when) these examples are run.<\/div>\n
\n
\n
d = 10000; %set the computational precision \"very\" high<\/span><\/span><\/div>\n<\/div>\n
\n
rng(1); %set the seed for the random number generator<\/span><\/span><\/div>\n<\/div>\n<\/div>\n
And with that, let's set up our collection of numbers. For those interested in learning about each number, check out the code block at the end of this post. If you're content with knowing these are just being used as examples, feel free to skip seeing the definitions.<\/div>\n
\n
\n
[num1, num2, num3, num4, num5, num6] = generateExamples(d);<\/span><\/div>\n<\/div>\n<\/div>\n

Continued Fractions<\/h2>\n
So what exactly is a continued fraction? As the name implies, a continued fraction is a fraction that continues<\/span>. What are fractions that don't continue? Most fractions you might think of, such as $ \\frac{11}{13} $.<\/div>\n
This isn't the only way to write this fraction though. One way, using Egyptian fractions<\/a>, is $ \\frac{11}{13} = \\frac{1}{2} + \\frac{1}{3} + \\frac{1}{78} $, which we can verify using MATLAB's Symbolic Math Toolbox<\/a> (which will also be used significantly later on):<\/div>\n
\n
\n
eFrac = sym(1\/2) + sym(1\/3) + sym(1\/78)<\/span><\/div>\n
\n
\n
eFrac\u00a0=<\/div>\n

<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n

Another way is $ \\frac{11}{13} = \\frac{1}{1 + \\frac{1}{ 5 + \\frac{1}{2}}} $. We can check this just as we did above:<\/div>\n
\n
\n
cFrac = 1\/(1 + 1\/(5 + sym(1\/2)))<\/span><\/div>\n
\n
\n
cFrac\u00a0=<\/div>\n

<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n

This is what we call a finite (simple) continued fraction:<\/span> continued fraction<\/span> because it continues, finite<\/span> because there are finitely many terms, and simple<\/span> because the numerators that appear are all 1. We also note that all of the coefficients we want should be positive integers; variations with negative integers also exist and the resulting continued fractions will have different properties.<\/div>\n
This begs two natural questions:<\/div>\n
    \n
  1. Are there continued fractions that have infinitely many terms?<\/li>\n
  2. Are there continued fractions with numerators other than 1?<\/li>\n<\/ol>\n
    Yes and yes, which we'll demonstrate by example:<\/div>\n
      \n
    1. The golden ratio<\/a> $ \\phi = \\frac{1 + \\sqrt{5}}{2} = 1+ \\frac{1}{1 + \\frac{1}{1 + \\cdots}} \\approx 1.618... $,<\/li>\n
    2. $ \\frac{4}{\\pi} = 1 + \\frac{1^2}{2 + \\frac{3^2}{2 + \\frac{5^2}{2 + \\cdots}}} \\approx 1.2732... $.<\/li>\n<\/ol>\n
      In this post we'll focus on infinite continued fractions, which we'll refer to as continued fractions for simplicity. Likewise, we won't touch on non-simple continued fractions, but rest assured they exist and often express interesting patterns and are fascinating mathematical objects to explore; for simplicity, whenever we say continued fraction<\/span> we're talking about a simple<\/span> continued fraction.<\/div>\n
      So how does one compute the continued fraction for a number? The general algorithm is a sort of Euclid's algorithm, wherein we take off the integer portion of a number, invert the fraction, and repeat. Let's see how this works with $ \\frac{11}{13} $:<\/div>\n
        \n
      1. $ \\frac{11}{13} $ has no integer part, so we have $ \\frac{11}{13} = 0 + \\frac{11}{13} $.<\/li>\n
      2. Next we invert the fractional part of the last expression, giving $ \\frac{11}{13} = 0 + \\frac{1}{\\frac{13}{11}} $.<\/li>\n
      3. For the newest \"denominator\", we repeat the process of taking out an integer part: since $ \\frac{13}{11} = 1 + \\frac{2}{11} $, we get $ \\frac{11}{13} = 0 + \\frac{1}{1 + \\frac{2}{11}} $.<\/li>\n
      4. Having isolated the integer part, we now invert the new fractional part $ \\frac{2}{11} $, giving $ \\frac{11}{13} = 0 + \\frac{1}{1 + \\frac{1}{\\frac{11}{2}}} $.<\/li>\n
      5. Since $ \\frac{11}{2} = 5 + \\frac{1}{2} $, we substitute this in to get the continued fraction $ \\frac{11}{13} = 0 + \\frac{1}{1 + \\frac{1}{5 + \\frac{1}{2}}} $.<\/li>\n<\/ol>\n
        It shouldn't be too surprising that any finite continued fraction corresponds to a rational number (just apply the algorithm backwards from the end of the continued fraction), but it may not be obvious that every rational number has a finite continued fraction. As a comparison, the decimal representation of $ \\frac{1}{2} $ is $ 0.5 $, while the decimal representation of $ \\frac{1}{3} $ is $ 0.333... $. It turns out that for continued fractions, there will be infinitely many coefficients if the number is irrational, and likewise if the number is irrational its continued fraction will have infinitely coefficients.<\/div>\n

        Infinite Continued Fractions<\/h3>\n
        We've already seen one example of an infinite continued fraction, $ \\phi = 1+ \\frac{1}{1 + \\frac{1}{1 + \\cdots}} $. There are two ways to find the coefficients: one involves the algorithm from above, while the other leverages the fact that \u03d5<\/span> is a solution to $ 0 = \\phi^2 - \\phi - 1 $. Let's see both methods in action:<\/div>\n
          \n
        1. Using the approximation $ \\phi =1.618033988749895 $, we get $ \\phi = 1 +0.618033988749895 $. Inverting the fraction gives $ \\phi = 1 + \\frac{1}{1.618033988749894} $, since $ \\frac{1}{0.618033988749895} \\approx 1.618033988749894 $. Notice that this new denominator is incredibly close to our original approximation for \u03d5<\/span>, and only differs in the last digit we've recorded. Continuing on, we get $ \\phi = 1 + \\frac{1}{1 + \\frac{1}{\\frac{1}{0.618033988749894}}} $, etc.<\/li>\n
        2. The second method doesn't rely on approximations, but includes some of the observations in the last method. In particular, we saw that $ \\frac{1}{\\phi - 1} \\approx \\phi $, which upon rearranging gives $ 1 = \\phi^2 - \\phi $, which is precisely the equation $ 0 = \\phi^2 - \\phi - 1 $. What we'll do now is start with the quadratic equation for \u03d5<\/span> and rearrange it to get us the continued fraction fairly easily. So, starting with $ 0 = \\phi^2 - \\phi - 1 $, we rearrange to get $ \\phi^2 = \\phi + 1 $, and divide by \u03d5<\/span> to end up with $ \\phi = \\frac{\\phi + 1}{\\phi} $. Cleaning up the right-hand side of this last equality gives $ \\phi = 1 + \\frac{1}{\\phi} $. What this expression suggests is that, since we have $ \\phi = $ something<\/span>, let's plug that something in to the left-hand side: $ \\phi = 1+ \\frac{1}{1 + \\frac{1}{\\phi}} $. Repeating this again gives $ \\phi =1 + \\frac{1}{1 + \\frac{1}{1+\\frac{1}{\\phi}}} $, and so on.<\/li>\n<\/ol>\n
          This second method can be used for a lot of roots of quadratic equations, though for other irrational numbers there may or may not be some nice properties we can exploit. Before discussing how the general algorithm works in general, and implementing it in MATLAB, we'll see one more example of finding a continued fraction.<\/div>\n
          Let's find the continued fraction for $ \\sqrt{2} $. First, recall that $ \\sqrt{2} $ is a solution to the quadratic equation $ x^2 - 2 = 0 $. We can try to write this in the same form that worked for \u03d5<\/span>, namely $ x = \\frac{2}{x} $, but this won't actually get us any closer. Instead, we're going to add 0 to $ \\sqrt{2} $ and hope for the best: $ \\sqrt{2} = \\sqrt{2} + 1 - 1 = 1 + \\left( \\sqrt{2} - 1\\right) $. This step is the first step of the general algorithm, since $ \\sqrt{2} \\approx 1.414 $, so we've successfully isolated the integer part. Next we invert the fractional part to get $ \\frac{1}{\\sqrt{2}-1} = \\frac{\\sqrt{2}+1}{(\\sqrt{2}-1)(\\sqrt{2}+1)} = \\sqrt{2} + 1 $. Thus, our continued fraction should look like $ \\sqrt{2} = 1 + \\frac{1}{1 + \\sqrt{2} } $. At this point we have the recurrence we need, and we can keep substituting $ \\sqrt{2} $ back in on the left-hand side to get $ \\sqrt{2} = 1 + \\frac{1}{1 + 1 + \\frac{1}{1 + \\sqrt{2}} } = 1 + \\frac{1}{2 + \\frac{1}{1 + \\sqrt{2}}} = 1 + \\frac{1}{2 + \\frac{1}{2 + \\frac{1}{ 1 + \\sqrt{2}}}} $, and so on.<\/div>\n<\/div>\n
          \n

          Continued Fraction Notation<\/h3>\n
          As you may have noticed, writing out a continued fraction can get fairly unwieldly. Since we're working with simple continued fractions, the only coefficients we have to keep track of are the \"integer pieces\" in each level of the continued fraction. For $ \\sqrt{2} $ the integer piece is 1, but all of the subsequent coefficients are 2. So that we don't have to write out $ \\sqrt{2} = 1 + \\frac{1}{ 2 + \\frac{1}{ 2 + \\frac{1}{2 + \\cdots} } } $ each time, we'll use the condensed notation $ \\sqrt{2} = [1; 2,2,...] $. The semicolon just indicates that the first coefficient is \"special\" since it's the actual integer part. As another example, $ \\phi =1 + \\frac{1}{1 + \\frac{1}{1+\\frac{1}{1 + \\cdots}}} = [1; 1, 1, ...] $.<\/div>\n
          As a third example, $ \\pi = [3; 7, 15, 1, 292, 1, 1,...] = 3 + \\frac{1}{7 + \\frac{1}{15 + \\frac{1}{1 + \\frac{1}{292 + \\frac{1}{1 + \\cdots}}}}} $. Next, we'll implement the algorithm for computing continued fractions and verify this last representation.<\/div>\n
          More generally, if $ a_0, a_1, a_2, ... $ are positive integers, then the continued fraction $ a_0 + \\frac{1}{a_1 + \\frac{1}{a_2 + \\cdots}} $ is written $ [a_0; a_1, a_2, ...] $.<\/div>\n

          Computing the Coefficients<\/h3>\n
          As a quick reminder, the algorithm to find the coefficients for a continued fraction involves iteratively (1) taking out the integer part of a quantity, (2) inverting the fractional part, and repeating this indefinitely; in the case of using a computer, we repeat this process until our level of precision is met or exceeded. Let's see what this looks for \u03d5<\/span> and \u03c0<\/span>.<\/div>\n

          Computing \u03d5<\/span>'s Continued Fraction<\/h3>\n
          As a test case to verify the algorithm, and its accuracy, we'll start with \u03d5<\/span>'s continued fraction coefficients. Recall that $ \\phi = [1; 1, 1, ...] $, so as soon as the algorithm returns a coefficient that is not<\/span> 1, we'll know that's around the limitation of what we can accurately compute.<\/div>\n
          \n
          \n
          phi = vpa((1 + sqrt(sym(5)))\/2,d) % define phi as a symbolic variable<\/span><\/span><\/div>\n
          \n
          \n

          phi\u00a0=\u00a0<\/span><\/p>\n

          1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939113748475408807538689175212663386222353693179318006076672635443338908659593958290563832266131992829026788067520876689250171169620703222104321626954862629631361443814975870122034080588795445474924618569536486444924104432077134494704956584678850987433944221254487706647809158846074998871240076521705751797883416625624940758906970400028121042762177111777805315317141011704666599146697987317613560067087480710131795236894275219484353056783002287856997829778347845878228911097625003026961561700250464338243776486102838312683303724292675263116533924731671112115881863851331620384005222165791286675294654906811317159934323597349498509040947621322298101726107059611645629909816290555208524790352406020172799747175342777592778625619432082750513121815628551222480939471234145170223735805772786160086883829523045926478780178899219902707769038953219681986151437803149974110692608867429622675756052317277752035361393621076738937645560606059216589466759551900400555908950229530942312482355212212415444006470340565734797663972394949946584578873039623090375033993856210242369025138680414577995698122445747178034173126453220416397232134044449487302315417676893752103068737880344170093954409627955898678723209512426893557309704509595684401755519881921802064052905518934947592600734852282101088194644544222318891319294689622002301443770269923007803085261180754519288770502109684249362713592518760777884665836150238913493333122310533923213624319263728910670503399282265263556209029798642472759772565508615487543574826471814145127000602389016207773224499435308899909501680328112194320481964387675863314798571911397815397807476150772211750826945863932045652098969855567814106968372884058746103378105444390943683583581381131168993855576975484149144534150912954070050194775486163075422641729394680367319805861833918328599130396072014455950449779212076124785645916160837059498786006970189409886400764436170933417270919143365013715766011480381430626238051432117348151005590134561011800790506381421527093085880928757034505078081454588199063361298279814117453392731208092897279222132980642946878242748740174505540677875708323731097591511776297844328474790817651809778726841611763250386121129143683437670235037111633072586988325871033632223810980901211019899176841491751233134015273384383723450093478604979294599158220125810459823092552872124137043614910205471855496118087642657651106054588147560443178479858453973128630162544876114852021706440411166076695059775783257039511087823082710647893902111569103927683845386333321565829659773103436032322545743637204124406408882673758433953679593123221343732099574988946995656473600729599983912881031974263125179714143201231127955189477817269141589117799195648125580018455065632952859859100090862180297756378925999164994642819302229355234667475932695165421402109136301819472270789012208728736170734864999815625547281137347987165695274890081443840532748378137824669174442296349147081570073525457070897726754693438226195468615331209533579238014609273510210119190218360675097308957528957746814229543394385493155339630380729169175846101460995055064803679304147236572039860073550760902317312501613204843583648177048481810991602442523271672190189334596378608787528701739359303013359011237102391712659047026349402830766876743638651327106280323174069317334482343564531850581353108549733350759966778712449058363675413289086240632456395357212524261170278028656043234942837301725574405837278267996031739364013287627701243679831144643694767053127249241047167001382478312865650649343418039004101780533950587724586655755229391582397084177298337282311525692609299594224000056062667867435792397245408481765197343626526894488855272027477874733598353672776140759171205132693448375299164998093602461784426757277679001919190703805220461232482391326104327191684512306023627893545432461769975753689041763650254785138246314658336383376023577899267298863216185839590363998183845827644912459809370430555596137973432613483049494968681089535696348281781288625364608420339465381944194571426668237183949183237090857485026656803989744066210536030640026081711266599541993687316094572288810920778822772036366844815325617284117690979266665522384688311371852991921631905201568631222820715599876468423552059285371757807656050367731309751912239738872246825805715974457404842987807352215984266766257807706201943040054255015831250301753409411719101929890384472503329880245014367968441694795954530459103138116218704567997866366174605957000344597011352518134600656553520347888117414994127482641521355677639403907103870881823380680335003804680017480822059109684420264464021877053401003180288166441530913939481564031928227854824145105031888251899700748622879421558957428202166570621880905780880503246769912972872103870736974064356674589202586565739785608595665341070359978320446336346485489497663885351045527298242290699848853696828046459745762651434359050938321243743333870516657149005907105670248879858043718151261004403814880407252440616429022478227152724112085065788838712493635106806365166743222327767755797399270376231914704732395512060705503992088442603708790843334261838413597078164829553714321961189503797714630007555975379570355227144931913217255644012830918050450089921870512118606933573153895935079030073672702331416532042340155374144268715405511647961143323024854404094069114561398730260395182816803448252543267385759005604320245372719291248645813334416985299391357478698957986439498023047116967157362283912018127312916589952759919220318372356827279385637331265479985912463275030060592567454979435088119295056854932593553187291418011364121874707526281068698301357605247194455932195535961045283031488391176930119658583431442489489856558425083410942950277197583352244291257364938075417113739243760143506829878493271299751228688196049835775158771780410697131966753477194792263651901633977128473907933611119140899830560336106098717178305543540356089529290818464143713929437813560482038947912574507707557510300242072662900180904229342494259060666141332287226980690145994511995478016399151412612525728280664331261657469388195106442167387180001100421848302580916543383749236411838885646851431500637319042951481469424314608952547072037405566913069220990804819452975110650464281054177552590951871318883591476599604131796020941530858553323877253802327276329773721431279682167162344211832018028814127474431688472184593927814354740999990722332030592629766112383279833169882539312620065037028844782866694044730794710476125586583752986236250999823233597155072338383324408152577819336426263043302658958170800451278873115935587747217256494700051636672577153920984095032745112153687300912199629522765913163709396860727134269262315475330437993316581107369643142171979434056391551210810813626268885697480680601169189417502722987415869917914534994624441940121978586013736608286907223651477139126874209665137875620591854328888341742920901563133283193575622089713765630978501563154982456445865424792935722828750608481453351352181729587932991171003247622205219464510536245051298843087134443950724426735146286179918323364598369637632722575691597239543830520866474742381511079273494836952396479268993698324917999502789500060459661313463363024949951480805329017902975182515875049007435187983511836032722772601717404535571658855578297291061958193517105548257930709100576358699019297217995168731175563144485648100220014254540554292734588371160209947945720823780436871894480563689182580244499631878342027491015335791072733625328906933474123802222011626277119308544850295419132004009998655666517756640953656197897818380451030356510131589458902871861086905893947136801484570018366495647203294334374298946427412551435905843484091954870152361403173913903616440198455051049121169792001201999605069949664030350863692903941007019450532016234872763232732449439630480890554251379723314751852070910250636859816795304818100739424531700238804759834323450414258431406361272109602282423378228090279765960777108493915174887316877713522390091171173509186006546200990249758527792542781659703834950580106261553336910937846597710529750223173074121778344189411845965861029801877874274456386696612772450384586052641510304089825777754474115332076407588167751497553804711629667771005876646159549677692705496239398570925507027406997814084312496536307186653371806058742242598165307052573834541577054292162998114917508611311765773172095615656478695474489271320608063545779462414531066983742113798168963823533304477883169339728728918103664083269856988254438516675862289930696434684897514840879039647604203610206021717394470263487633654393195229077383616738981178124248365578105034169451563626043003665743108476654877780128577923645418522447236171374229255841593135612866371670328072171553392646325730673063910854108868085742838588280602303341408550390973538726134511962926415995212789311354431460152730902553827104325966226743903745563612286139078319433570590038148700898661315398195857442330441970856696722293142730741384882788975588860799738704470203166834856941990965480298249319817657926829855629723010682777235162740783807431877827318211919695280051608791572128826337968231272562870001500182929757729993579094919640763442861575713544427898383040454702710194580042582021202344580630345033658147218549203679989972935353919681213319516537974539911149424445183033858841290401817818821376006659284941367754317451605409387110368715211640405821934471204482775960541694864539878326269548013915019038995931306703186616706637196402569286713887146631189192685682691995276457997718278759460961617218868109454651578869122410609814197268619255478789926315359472922825080542516906814010781796021885330762305563816316401922454503257656739259976517530801427160714308718862859836037465057134204670083432754230277047793311183666903232885306873879907135900740304907459889513647687608678443238248218930617570319563803230819719363567274196438726258706154330729637038127515170406005057594882723856345156390526577104264594760405569509598408889037620799566388017861855915944111725092313279771138<\/span><\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          tempR = phi; % tempR is the number we have at the latest stage of the algorithm<\/span><\/span><\/div>\n<\/div>\n
          \n
          numTerms = 100; % let's start with 1000 coefficients computed<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %initialize an array of zeros so we have space allocated at the start for<\/span><\/span><\/div>\n<\/div>\n
          \n
          %everything<\/span><\/span><\/div>\n<\/div>\n
          \n
          A = zeros(1,numTerms);<\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %the main for loop where we compute the coefficients<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:numTerms<\/span><\/div>\n<\/div>\n
          \n
          %compute the integer piece of wherever we're at in the process<\/span><\/span><\/div>\n<\/div>\n
          \n
          A(ii) = floor(tempR); % floor finds the nearest integer below tempR<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %set up to get the next integer piece, etc.<\/span><\/span><\/div>\n<\/div>\n
          \n
          tempR = 1\/(tempR - floor(tempR)); % remove the integer piece and invert the fraction<\/span><\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          % let's see what was computed<\/span><\/span><\/div>\n<\/div>\n
          \n
          A <\/span><\/div>\n
          \n
          \n
          \n
          \n
          A = <\/span>1\u00d7100<\/span><\/div>\n<\/div>\n
          \n
          1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          The 96th coefficient is 2, not 1, so with this precision we're getting ~100 coefficients accurately. While this may not seem like too many coefficients, using the continued fraction with these 95 coefficients would let us compute \u03d5<\/span> to within about 40 digits. This follow's from Loch's theorem<\/a>, since for \u03d5<\/span> it takes 2.39 continued fraction coefficients to approximate 1 digit of accuracy for \u03d5<\/span>.<\/div>\n
          Next, let's do the same with $ \\pi = [3; 7, 15, 1, 292, 1, 1,...] $.<\/div>\n
          \n
          \n
          symPi = vpa(sym(pi),d)<\/span><\/div>\n
          \n
          \n

          symPi\u00a0=\u00a0<\/span><\/p>\n

          3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420198938095257201065485863278865936153381827968230301952035301852968995773622599413891249721775283479131515574857242454150695950829533116861727855889075098381754637464939319255060400927701671139009848824012858361603563707660104710181942955596198946767837449448255379774726847104047534646208046684259069491293313677028989152104752162056966024058038150193511253382430035587640247496473263914199272604269922796782354781636009341721641219924586315030286182974555706749838505494588586926995690927210797509302955321165344987202755960236480665499119881834797753566369807426542527862551818417574672890977772793800081647060016145249192173217214772350141441973568548161361157352552133475741849468438523323907394143334547762416862518983569485562099219222184272550254256887671790494601653466804988627232791786085784383827967976681454100953883786360950680064225125205117392984896084128488626945604241965285022210661186306744278622039194945047123713786960956364371917287467764657573962413890865832645995813390478027590099465764078951269468398352595709825822620522489407726719478268482601476990902640136394437455305068203496252451749399651431429809190659250937221696461515709858387410597885959772975498930161753928468138268683868942774155991855925245953959431049972524680845987273644695848653836736222626099124608051243884390451244136549762780797715691435997700129616089441694868555848406353422072225828488648158456028506016842739452267467678895252138522549954666727823986456596116354886230577456498035593634568174324112515076069479451096596094025228879710893145669136867228748940560101503308617928680920874760917824938589009714909675985261365549781893129784821682998948722658804857564014270477555132379641451523746234364542858444795265867821051141354735739523113427166102135969536231442952484937187110145765403590279934403742007310578539062198387447808478489683321445713868751943506430218453191048481005370614680674919278191197939952061419663428754440643745123718192179998391015919561814675142691239748940907186494231961567945208095146550225231603881930142093762137855956638937787083039069792077346722182562599661501421503068038447734549202605414665925201497442850732518666002132434088190710486331734649651453905796268561005508106658796998163574736384052571459102897064140110971206280439039759515677157700420337869936007230558763176359421873125147120532928191826186125867321579198414848829164470609575270695722091756711672291098169091528017350671274858322287183520935396572512108357915136988209144421006751033467110314126711136990865851639831501970165151168517143765761835155650884909989859982387345528331635507647918535893226185489632132933089857064204675259070915481416549859461637180270981994309924488957571282890592323326097299712084433573265489382391193259746366730583604142813883032038249037589852437441702913276561809377344403070746921120191302033038019762110110044929321516084244485963766983895228684783123552658213144957685726243344189303968642624341077322697802807318915441101044682325271620105265227211166039666557309254711055785376346682065310989652691862056476931257058635662018558100729360659876486117910453348850346113657686753249441668039626579787718556084552965412665408530614344431858676975145661406800700237877659134401712749470420562230538994561314071127000407854733269939081454664645880797270826683063432858785698305235808933065757406795457163775254202114955761581400250126228594130216471550979259230990796547376125517656751357517829666454779174501129961489030463994713296210734043751895735961458901938971311179042978285647503203198691514028708085990480109412147221317947647772622414254854540332157185306142288137585043063321751829798662237172159160771669254748738986654949450114654062843366393790039769265672146385306736096571209180763832716641627488880078692560290228472104031721186082041900042296617119637792133757511495950156604963186294726547364252308177036751590673502350728354056704038674351362222477158915049530984448933309634087807693259939780541934144737744184263129860809988868741326047215695162396586457302163159819319516735381297416772947867242292465436680098067692823828068996400482435403701416314965897940924323789690706977942236250822168895738379862300159377647165122893578601588161755782973523344604281512627203734314653197777416031990665541876397929334419521541341899485444734567383162499341913181480927777103863877343177207545654532207770921201905166096280490926360197598828161332316663652861932668633606273567630354477628035045077723554710585954870279081435624014517180624643626794561275318134078330336254232783944975382437205835311477119926063813346776879695970309833913077109870408591337464144282277263465947047458784778720192771528073176790770715721344473060570073349243693113835049316312840425121925651798069411352801314701304781643788518529092854520116583934196562134914341595625865865570552690496520985803385072242648293972858478316305777756068887644624824685792603953527734803048029005876075825104747091643961362676044925627420420832085661190625454337213153595845068772460290161876679524061634252257719542916299193064553779914037340432875262888963995879475729174642635745525407909145135711136941091193932519107602082520261879853188770584297259167781314969900901921169717372784768472686084900337702424291651300500516832336435038951702989392233451722013812806965011784408745196012122859937162313017114448464090389064495444006198690754851602632750529834918740786680881833851022833450850486082503930213321971551843063545500766828294930413776552793975175461395398468339363830474611996653858153842056853386218672523340283087112328278921250771262946322956398989893582116745627010218356462201349671518819097303811980049734072396103685406643193950979019069963955245300545058068550195673022921913933918568034490398205955100226353536192041994745538593810234395544959778377902374216172711172364343543947822181852862408514006660443325888569867054315470696574745855033232334210730154594051655379068662733379958511562578432298827372319898757141595781119635833005940873068121602876496286744604774649159950549737425626901049037781986835938146574126804925648798556145372347867330390468838343634655379498641927056387293174872332083760112302991136793862708943879936201629515413371424892830722012690147546684765357616477379467520049075715552781965362132392640616013635815590742202020318727760527721900556148425551879253034351398442532234157623361064250639049750086562710953591946589751413103482276930624743536325691607815478181152843667957061108615331504452127473924544945423682886061340841486377670096120715124914043027253860764823634143346235189757664521641376796903149501910857598442391986291642193994907236234646844117394032659184044378051333894525742399508296591228508555821572503107125701266830240292952522011872676756220415420516184163484756516999811614101002996078386909291603028840026910414079288621507842451670908700069928212066041837180653556725253256753286129104248776182582976515795984703562226293486003415872298053498965022629174878820273420922224533985626476691490556284250391275771028402799806636582548892648802545661017296702664076559042909945681506526530537182941270336931378517860904070866711496558343434769338578171138645587367812301458768712660348913909562009939361031029161615288138437909904231747336394804575931493140529763475748119356709110137751721008031559024853090669203767192203322909433467685142214477379393751703443661991040337511173547191855046449026365512816228824462575916333039107225383742182140883508657391771509682887478265699599574490661758344137522397096834080053559849175417381883999446974867626551658276584835884531427756879002909517028352971634456212964043523117600665101241200659755851276178583829204197484423608007193045761893234922927965019875187212726750798125547095890455635792122103334669749923563025494780249011419521238281530911407907386025152274299581807247162591668545133312394804947079119153267343028244186041426363954800044800267049624820179289647669758318327131425170296923488962766844032326092752496035799646925650493681836090032380929345958897069536534940603402166544375589004563288225054525564056448246515187547119621844396582533754388569094113031509526179378002974120766514793942590298969594699556576121865619673378623625612521632086286922210327488921865436480229678070576561514463204692790682120738837781423356282360896320806822246801224826117718589638140918390367367222088832151375560037279839400415297002878307667094447456013455641725437090697939612257142989467154357846878861444581231459357198492252847160504922124247014121478057345510500801908699603302763478708108175450119307141223390866393833952942578690507643100638351983438934159613185434754649556978103829309716465143840700707360411237359984345225161050702705623526601276484830840761183013052793205427462865403603674532865105706587488225698157936789766974220575059683440869735020141020672358502007245225632651341055924019027421624843914035998953539459094407046912091409387001264560016237428802109276457931065792295524988727584610126483699989225695968815920560010165525637568<\/span><\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          tempR = symPi; % tempR is the number we have at the latest stage of the algorithm<\/span><\/span><\/div>\n<\/div>\n
          \n
          numTerms = 100; % let's start with 100 coefficients computed<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %initialize an array of zeros so we have space allocated at the start for<\/span><\/span><\/div>\n<\/div>\n
          \n
          %everything<\/span><\/span><\/div>\n<\/div>\n
          \n
          A = zeros(1,numTerms);<\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %the main for loop where we compute the coefficients<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:numTerms<\/span><\/div>\n<\/div>\n
          \n
          %compute the integer piece of wherever we're at in the process<\/span><\/span><\/div>\n<\/div>\n
          \n
          A(ii) = floor(tempR); % floor finds the nearest integer below tempR<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %set up to get the next integer piece, etc.<\/span><\/span><\/div>\n<\/div>\n
          \n
          tempR = 1\/(tempR - floor(tempR)); % remove the integer piece and invert the fraction<\/span><\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          % let's see what was computed<\/span><\/span><\/div>\n<\/div>\n
          \n
          A <\/span><\/div>\n
          \n
          \n
          \n
          \n
          A = <\/span>1\u00d7100<\/span><\/div>\n<\/div>\n
          \n
          3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2 1 84 2 1 1 15 3 13 1 4 2 6 6 99 1 4 4 2 4 8 1 3 4 2 12 2 3 1 8 1<\/div>\n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          Just as expected! In case there's any doubt, the continued fraction coefficients for \u03c0<\/span> can be externally verified using the Online Encyclopedia of Integer Sequences as A001203<\/a>.<\/div>\n
          Since we'll use the algorithm over and over, it's been encoded in the function getCF<\/span> available at the end of this post.<\/div>\n
          We can also visualize these coefficients using bar graphs, where the bars are only used a quick graphical aid to get a sense of the distribution, and in particular how large, the coefficients can get. Here we're using a higher precision than above since we're going to start going big with the number of coefficients we compute.<\/div>\n
          \n
          \n
          piCF = getCF(vpa(sym(pi),d),100);<\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % first 10 coefficients<\/span><\/span><\/div>\n<\/div>\n
          \n
          bar(piCF(1:10))<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 10 continued fraction coefficients of \\pi\"<\/span>)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %first 100 coefficients<\/span><\/span><\/div>\n<\/div>\n
          \n
          bar(piCF)<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 100 continued fraction coefficients of \\pi\"<\/span>)<\/span><\/div>\n
          \n
          \n
          <\/div>\n
          \"\"<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          Before getting to the core of this blog post - statistical properties of continued fraction coefficients - let's see what some of the other numbers we defined at the start of this post look like.<\/div>\n
          \n
          \n
          tic<\/span><\/div>\n<\/div>\n
          \n
          num1CF = getCF(num1,100);<\/span><\/div>\n<\/div>\n
          \n
          toc<\/span><\/div>\n
          \n
          \n
          Elapsed time is 0.133267 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          bar(num1CF)<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 100 continued fraction coefficients of the Fibonacci constant\"<\/span>)<\/span><\/div>\n
          \n
          \n
          \"\"<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          tic<\/span><\/div>\n<\/div>\n
          \n
          num2CF = getCF(num2,100);<\/span><\/div>\n<\/div>\n
          \n
          toc<\/span><\/div>\n
          \n
          \n
          Elapsed time is 0.125300 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          bar(num2CF)<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 100 continued fraction coefficients of Champernowne's constant\"<\/span>)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          tic<\/span><\/div>\n<\/div>\n
          \n
          num3CF = getCF(num3,100);<\/span><\/div>\n<\/div>\n
          \n
          toc<\/span><\/div>\n
          \n
          \n
          Elapsed time is 0.128072 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          bar(num3CF)<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 100 continued fraction coefficients of the Prime constant\"<\/span>)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          tic<\/span><\/div>\n<\/div>\n
          \n
          num4CF = getCF(num4,100);<\/span><\/div>\n<\/div>\n
          \n
          toc<\/span><\/div>\n
          \n
          \n
          Elapsed time is 0.134937 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          bar(num4CF)<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 100 continued fraction coefficients of the Squares constant\"<\/span>)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          tic<\/span><\/div>\n<\/div>\n
          \n
          num5CF = getCF(num5,100);<\/span><\/div>\n<\/div>\n
          \n
          toc<\/span><\/div>\n
          \n
          \n
          Elapsed time is 0.125739 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          bar(num5CF)<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 100 continued fraction coefficients of Liouville's constant\"<\/span>)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          tic<\/span><\/div>\n<\/div>\n
          \n
          num6CF = getCF(num6,100);<\/span><\/div>\n<\/div>\n
          \n
          toc<\/span><\/div>\n
          \n
          \n
          Elapsed time is 0.123741 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          bar(num6CF)<\/span><\/div>\n<\/div>\n
          \n
          title(\"First 100 continued fraction coefficients of a random decimal\"<\/span>)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          Some of the coefficients get quite large, but most tend to stay closer (relatively speaking) to 1. While many of these coefficients may seem random, they have some pretty remarkable statistical properties as we'll explore next.<\/div>\n

          Statistics of the Coefficients<\/h2>\n
          Since we've mentioned statistics, let's start by taking means of the coefficients and seeing what pops out. In particular, since we've been constructing continued fractions one coefficient at a time, we're going to take running means <\/span>of the coefficients and see what their long-term behaviour is. Recall that the mean of numbers $ a_1, a_2, ..., a_n $ is $ \\frac{a_1 + \\cdots + a_n}{n} $.<\/div>\n
          We'll start by considering the running means of \u03c0<\/span>'s coefficients.<\/div>\n
          \n
          \n
          Means = zeros(1,100); % initialize an array to hold the means <\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(ii) = mean(piCF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          Means<\/span><\/div>\n
          \n
          \n
          \n
          \n
          Means = <\/span>1\u00d7100<\/span><\/div>\n<\/div>\n
          \n
          3.0000 5.0000 8.3333 6.5000 63.6000 53.1667 45.7143 40.1250 35.8889 32.4000 29.7273 27.3333 26.3077 24.5714 23.0000 21.6250 20.4706 19.4444 18.5263 17.7000 16.9048 19.9545 19.1739 18.4167 17.7200 17.6154 17.0741 16.9286 16.3793 15.9667 15.5161 15.2188 14.9394 17.4118 16.9429 16.5833 16.2432 15.8684 15.5641 15.3750 15.0244 14.7381 14.4884 14.2045 14.1556 13.8913 13.6596 13.3958 13.2857 13.0400<\/div>\n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          Let's visualize these results with a line plot:<\/div>\n
          \n
          \n
          plot(Means)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          The plot suggests that the means are tending towards the value 10, and inspecting the array of means tells us that the actual, apparent, limiting value is about 8. Let's do this same procedure for the other continued fractions we computed to use as examples and plot them all on the same graph.<\/div>\n
          \n
          \n
          Means = zeros(7,100); % initialize an array to hold the means <\/span><\/span><\/div>\n<\/div>\n
          \n
          % this time 7 rows since we have pi and 6 other examples <\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the means for pi <\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(1,ii) = mean(piCF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the means for the Fibonacci constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(2,ii) = mean(num1CF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the means for Champernowne's constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(3,ii) = mean(num2CF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the means for the Prime constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(4,ii) = mean(num3CF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the means for the Squares constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(5,ii) = mean(num4CF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the means for Liouville's constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(6,ii) = mean(num5CF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the means for random decimal<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(7,ii) = mean(num6CF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % plot all lines together<\/span><\/span><\/div>\n<\/div>\n
          \n
          plot(Means')<\/span><\/div>\n<\/div>\n
          \n
          ylim([0 100])<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          We can generate even more examples, and evidence for conjectures, by constructing more random decimals and plotting their running means.<\/div>\n
          \n
          \n
          numExperiments = 100; % the number of random numbers to generate<\/span><\/span><\/div>\n<\/div>\n
          \n
          Means = zeros(numExperiments,100); % initialize an array to hold the means <\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          tic % to time the computation<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>jj = 1:numExperiments<\/span><\/div>\n<\/div>\n
          \n
          %copying the code from above that generates a random decimal<\/span><\/span><\/div>\n<\/div>\n
          \n
          tempDigits = randi(10,[1 d]) - 1; %generate an array of randomly selected digits between 1 to 9<\/span><\/span><\/div>\n<\/div>\n
          \n
          temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
          \n
          temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
          \n
          randDec = vpa(\"0.\"<\/span>+temp,d);<\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % compute its CF<\/span><\/span><\/div>\n<\/div>\n
          \n
          randCF = getCF(randDec,100);<\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % compute its running means<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Means(jj,ii) = mean(randCF(1:ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          end <\/span><\/span><\/div>\n<\/div>\n
          \n
          toc % display how much time elapsed for the data generation<\/span><\/span><\/div>\n
          \n
          \n
          Elapsed time is 13.492355 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          plot(Means')<\/span><\/div>\n<\/div>\n
          \n
          ylim([0 100])<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          For most of the continued fractions, the means seem to tend to a value between 5 and 20, and all seem to be greater than 2.5. This is interesting, and with more coefficients we may generate more evidence for this observation - we may even generate a<\/div>\n
          Conjecture: <\/span>the running means of (most) continued fractions coefficients is at least 2.5<\/div>\n
          (and prove it, later on).<\/div>\n
          Here though we'll pivot, and instead of exploring means, we'll explore geometric <\/span>means of the coefficients. For numbers $ a_1, a_2, ..., a_n $, their geometric mean<\/span> is the quantity $ (a_1a_2 \\cdots a_n)^{1\/n} $.<\/div>\n
          Let's start by computing the running geometric means of the continued fraction coefficients of \u03c0<\/span>.<\/div>\n
          \n
          \n
          Geo = ones(1,99); % initialize the array to store geometric means<\/span><\/span><\/div>\n<\/div>\n
          \n
          % we're taking 1 less than the size of the CF since we're ignoring the<\/span><\/span><\/div>\n<\/div>\n
          \n
          % initial integer part<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          %the main for loop, where we loop through and compute the geometric means<\/span><\/span><\/div>\n<\/div>\n
          \n
          %of the sequences<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 2:length(piCF)<\/span><\/div>\n<\/div>\n
          \n
          Geo(ii-1) = prod(nthroot(piCF(2:ii),ii)); %note the indexing<\/span><\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          Geo<\/span><\/div>\n
          \n
          \n
          \n
          \n
          Geo = <\/span>1\u00d799<\/span><\/div>\n<\/div>\n
          \n
          2.6458 4.7177 3.2011 7.8943 5.5945 4.3746 3.6377 3.4037 3.0113 3.0103 2.7462 3.1127 3.0159 2.8019 2.6272 2.5854 2.5488 2.5164 2.4877 2.3821 2.8008 2.7601 2.6458 2.5448 2.7245 2.7342 2.8908 2.7869 2.8207 2.7895 2.8571 2.9221 3.2411 3.1340 3.1553 3.1756 3.1372 3.1568 3.2311 3.1399 3.1365 3.1543 3.1218 3.2167 3.1836 3.1796 3.1039 3.1644 3.0923 3.0246<\/div>\n<\/div>\n<\/div>\n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          plot(Geo)<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          The coefficients seem to converge to a value just below 3. Let's see what the other example numbers indicate.<\/div>\n
          \n
          \n
          Geo = ones(7,99); % initialize an array to hold the geometric means <\/span><\/span><\/div>\n<\/div>\n
          \n
          % this time 7 rows since we have pi and 6 other examples <\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the geometric means for pi <\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(1,ii) = prod(nthroot(piCF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the geometric means for the Fibonacci constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(2,ii) = prod(nthroot(num1CF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the geometric means for Champernowne's constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(3,ii) = prod(nthroot(num2CF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the geometric means for the Prime constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(4,ii) = prod(nthroot(num3CF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the geometric means for the Squares constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(5,ii) = prod(nthroot(num4CF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the geometric means for Liouville's constant<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(6,ii) = prod(nthroot(num5CF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % add the geometric means for random decimal<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 1:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(7,ii) = prod(nthroot(num6CF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % plot all lines together<\/span><\/span><\/div>\n<\/div>\n
          \n
          plot(Geo')<\/span><\/div>\n<\/div>\n
          \n
          ylim([0 10])<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          We can also run the same experiment as above, generating a bunch of random decimals and seeing what their running geometric means are.<\/div>\n
          \n
          \n
          numExperiments = 100; % the number of random numbers to generate<\/span><\/span><\/div>\n<\/div>\n
          \n
          Geo = zeros(numExperiments,99); % initialize an array to hold the geometric means <\/span><\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          tic % to time the computation<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>jj = 1:numExperiments<\/span><\/div>\n<\/div>\n
          \n
          %copying the code from above that generates a random decimal<\/span><\/span><\/div>\n<\/div>\n
          \n
          tempDigits = randi(10,[1 d]) - 1; %generate an array of randomly selected digits between 1 to 9<\/span><\/span><\/div>\n<\/div>\n
          \n
          temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
          \n
          temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
          \n
          randDec = vpa(\"0.\"<\/span>+temp,d);<\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % compute its CF<\/span><\/span><\/div>\n<\/div>\n
          \n
          randCF = getCF(randDec,100);<\/span><\/div>\n<\/div>\n
          \n
          <\/div>\n<\/div>\n
          \n
          % compute its running geometric means<\/span><\/span><\/div>\n<\/div>\n
          \n
          for <\/span>ii = 2:100<\/span><\/div>\n<\/div>\n
          \n
          Geo(jj,ii-1) = prod(nthroot(randCF(2:ii),ii));<\/span><\/div>\n<\/div>\n
          \n
          end<\/span><\/span><\/div>\n<\/div>\n
          \n
          end <\/span><\/span><\/div>\n<\/div>\n
          \n
          toc % display how much time elapsed for the data generation<\/span><\/span><\/div>\n
          \n
          \n
          Elapsed time is 14.002318 seconds.<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          \n
          plot(Geo')<\/span><\/div>\n<\/div>\n
          \n
          ylim([0 10])<\/span><\/div>\n
          \n
          \n
          <\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
          Remarkably, the running geometric means tend towards a value between 2 and 3 much more noticeably than the running means did. Of course with more continued fraction coefficients we'd get more data, and hopefully stronger evidence, for a conjecture, like<\/div>\n
          Conjecture:<\/span> the running geometric means of (most) continued fraction coefficients converge to 2.5.<\/div>\n
          The value 2.5 was chosen as the visual mean of the lines in our experiments. So the question remains: is our conjecture true?<\/div>\n
          Unfortunately it's not actually true. The actual value is $ \\approx 2.68545... $<\/div>\n
          This result is known as Khinchin's constant<\/a>, and was actually conjectured and proven in the 1930s. What's even wilder is the exact value is known, and is the infinite product$ \\prod_{r=1}^\\infty \\left( 1 + \\frac{1}{r(r+2)}\\right)^{\\log_2 r} $. The proof of this goes well beyond this post and involves more advanced mathematical techniques developed to study dynamical systems, but the linked Wikipedia article is a great starting point for intrepid explorers curious to learn more.<\/div>\n
          How does this connect to the exploration of means? Well, the AM-GM inequality<\/a> tells us that the arithmetic mean (the usual mean) is always greater than or equal to the geometric mean of the same set of numbers. Since Khinchin's constant is what most running geometric means converge to, the running (arithmetic) means of continued fraction coefficients should converge to values greater than Khinchin's constant.<\/div>\n

          Other Topics to Explore<\/h2>\n
          This discussion ignored many (arguably) more standard topics when it comes to continued fractions, like:<\/div>\n
            \n
          1. convergents, or what you get when you truncate a continued fraction after finitely many terms,<\/li>\n
          2. rates of convergence for convergents, or why $ \\frac{22}{7} $ and $ \\frac{355}{113} $ are excellent approximations of \u03c0<\/span>,<\/li>\n
          3. the fact that continued fractions are the \"best\" rational approximants for numbers,<\/li>\n
          4. MATLAB's built-in rational approximation tools<\/a>, which can find continued fractions with possibly negative exponents,<\/li>\n
          5. Loch's theorem<\/a>, which informs how many coefficients (or levels) of a continued fraction you need for another digit of accuracy (~1 coefficient per digit),<\/li>\n
          6. how roots of quadratic polynomials have periodic continued fractions<\/a> (the coefficients repeat), and are the only kinds of irrational numbers with periodic coefficients,<\/li>\n<\/ol>\n
            and so on. Indeed, when looking at convergents a result similar to Khinchin's constant exists called Levy's constant<\/a>.<\/div>\n
            There even exist (non-simple) continued fraction expansions for some well-known functions, like<\/div>\n
            $ \\sin(x) = \\frac{x}{1 + \\frac{x^2}{2\\cdot 3 - x^2 + \\frac{2\\cdot 3 x^2}{4\\cdot 5 - x^2 + \\frac{4\\cdot 5 x^2}{6\\cdot 7 - x^2 + ...}} } } $. How this expression is obtained goes far beyond the discussion here, but hopefully serves as motivation to keep exploring these fascinating objects (the Wikipedia page on Euler's continued fraction<\/a> might be a good starting point).<\/div>\n
            Happy exploring!<\/div>\n
            \n
            \n
            function <\/span>[num1, num2, num3, num4, num5, num6] = generateExamples(d)<\/span><\/div>\n<\/div>\n
            \n
            %%%%%%%%%%<\/span><\/span><\/div>\n<\/div>\n
            \n
            %Number - Fibonacci's Constant<\/span><\/span><\/div>\n<\/div>\n
            \n
            %0.1123581321..., constructed by concatenating the Fibonacci numbers 1, 1, 2, 3, 5...<\/span><\/span><\/div>\n<\/div>\n
            \n
            numTerms = 70; %set the number of terms in the Fibonacci sequence to use<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits = zeros(1,numTerms); %generate an empty array to store the first d terms of the Fibonacci sequence<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits(1) = 1;%initialize the initial condition<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits(2) = 1;%initialize the initial condition<\/span><\/span><\/div>\n<\/div>\n
            \n
            %the main for loop where we fill in the Fib. sequence<\/span><\/span><\/div>\n<\/div>\n
            \n
            for <\/span>ii = 3:numTerms<\/span><\/div>\n<\/div>\n
            \n
            tempDigits(ii) = tempDigits(ii-1) + tempDigits(ii-2);<\/span><\/div>\n<\/div>\n
            \n
            end<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
            \n
            num1 = vpa(\"0.\"<\/span>+temp,d); % tell MATLAB to use Variable Precision Arithmetic (vpa) for high precision<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %%%%%%%%%%<\/span><\/span><\/div>\n<\/div>\n
            \n
            %Number 2 - Champernowne's Constant<\/span><\/span><\/div>\n<\/div>\n
            \n
            %0.1234567891011..., constructed by concatenating the digits 1, 2, 3, ...<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits = 1:d; %generate an array of the first d digits<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
            \n
            num2 = vpa(\"0.\"<\/span>+temp,d); % tell MATLAB to use Variable Precision Arithmetic (vpa) for high precision<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %%%%%%%%%%<\/span><\/span><\/div>\n<\/div>\n
            \n
            %Number 3 - Prime Constant<\/span><\/span><\/div>\n<\/div>\n
            \n
            %0.235711..., formed by concatenating the prime numbers 2, 3, 5, 7, 11...<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits = primes(d); %generate an array of the primes below d<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
            \n
            num3 = vpa(\"0.\"<\/span>+temp,d); % tell MATLAB to use Variable Precision Arithmetic (vpa) for high precision<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %%%%%%%%%%<\/span><\/span><\/div>\n<\/div>\n
            \n
            %Number 4 - Squares Concatenated Constant<\/span><\/span><\/div>\n<\/div>\n
            \n
            %0.1491625..., formed by concatenating the integers squared 1, 4, 9, 16...<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits = 1:d; %generate an array of the first d digits<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits = tempDigits.^2; %square each of the terms in the array (hence the .^)<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
            \n
            num4 = vpa(\"0.\"<\/span>+temp,d); % tell MATLAB to use Variable Precision Arithmetic (vpa) for high precision<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %%%%%%%%%%<\/span><\/span><\/div>\n<\/div>\n
            \n
            %Number 5 - Liouville's Constant<\/span><\/span><\/div>\n<\/div>\n
            \n
            %0.1100010000000000000000010... formed by adding a 1 in the decimal place and zeros elsewhere.<\/span><\/span><\/div>\n<\/div>\n
            \n
            numTerms = 8; %number of factorial terms to use<\/span><\/span><\/div>\n<\/div>\n
            \n
            LiouvilleIndeces = factorial(1:numTerms); %first identify where the ones should go<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits = zeros(1,max(LiouvilleIndeces)+1); %generate an array of the first d digits<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits(LiouvilleIndeces) = 1; %square each of the terms in the array (hence the .^)<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
            \n
            num5 = vpa(\"0.\"<\/span>+temp,d); % tell MATLAB to use Variable Precision Arithmetic (vpa) for high precision<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %%%%%%%%%%<\/span><\/span><\/div>\n<\/div>\n
            \n
            % Number 6 - Random Decimal<\/span><\/span><\/div>\n<\/div>\n
            \n
            % A number whose decimal digits are all uniformly chosen between <\/span><\/span><\/div>\n<\/div>\n
            \n
            tempDigits = randi(10,[1 d]) - 1; %generate an array of randomly selected digits between 1 to 9<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = num2str(tempDigits); %convert the array to a string<\/span><\/span><\/div>\n<\/div>\n
            \n
            temp = strrep(temp,\" \"<\/span>,\"\"<\/span>); %remove empty space in the string<\/span><\/span><\/div>\n<\/div>\n
            \n
            num6 = vpa(\"0.\"<\/span>+temp,d); % tell MATLAB to use Variable Precision Arithmetic (vpa) for high precision<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            end<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            function <\/span>A = getCF(r,numTerms)<\/span><\/div>\n<\/div>\n
            \n
            % get the continued fraction coefficients<\/span><\/span><\/div>\n<\/div>\n
            \n
            %<\/span><\/span><\/div>\n<\/div>\n
            \n
            % the input r is the number we want the CF of<\/span><\/span><\/div>\n<\/div>\n
            \n
            % numTerms is however many coefficients we want to compute<\/span><\/span><\/div>\n<\/div>\n
            \n
            %<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempR = r; % tempR is the number we have at the latest stage of the algorithm<\/span><\/span><\/div>\n<\/div>\n
            \n
            % if numTerms isn't specified, default to computing 10 coefficients<\/span><\/span><\/div>\n<\/div>\n
            \n
            if <\/span>~exist('numTerms'<\/span>,'var'<\/span>)<\/span><\/div>\n<\/div>\n
            \n
            numTerms = 10;<\/span><\/div>\n<\/div>\n
            \n
            end<\/span><\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %initialize an array of zeros so we have space allocated at the start for<\/span><\/span><\/div>\n<\/div>\n
            \n
            %everything<\/span><\/span><\/div>\n<\/div>\n
            \n
            A = zeros(1,numTerms);<\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %the main for loop where we compute the coefficients<\/span><\/span><\/div>\n<\/div>\n
            \n
            for <\/span>ii = 1:numTerms<\/span><\/div>\n<\/div>\n
            \n
            %compute the integer piece of wherever we're at in the process<\/span><\/span><\/div>\n<\/div>\n
            \n
            A(ii) = floor(tempR);<\/span><\/div>\n<\/div>\n
            \n
            <\/div>\n<\/div>\n
            \n
            %set up to get the next integer piece, etc.<\/span><\/span><\/div>\n<\/div>\n
            \n
            tempR = 1\/(tempR - floor(tempR)); % remove the integer piece and invert the fraction<\/span><\/span><\/div>\n<\/div>\n
            \n
            end<\/span><\/span><\/div>\n<\/div>\n
            \n
            end<\/span><\/span><\/div>\n<\/div>\n<\/div>\n<\/div>\n","protected":false},"excerpt":{"rendered":"
            <\/div>\n

            \nToday's guest article is by Wesley Hamilton, a STEM Outreach engineer here at MathWorks. See his earlier piece on solving Sudoku puzzles.
            \nIn this blog post we'll explore the fascinating world of... read more >><\/a><\/p>\n","protected":false},"author":69,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/9470"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/users\/69"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/comments?post=9470"}],"version-history":[{"count":14,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/9470\/revisions"}],"predecessor-version":[{"id":9569,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/9470\/revisions\/9569"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/media?parent=9470"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/categories?post=9470"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/tags?post=9470"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}