A spigot stream for Sqrt[n]
A spigot stream for Sqrt[n]
Although the dragon scroller is great for right now, it’s not the only digital sequence we would like to display on our new peripheral. Irrational algebraic numbers are also a good target, for example:
In[]:=
N[Sqrt[2],50]
Out[]=
1.4142135623730950488016887242096980785696718753769
To get a procedure for spigot streaming these digits out, let’s consider the following equation:
n==(a+b+c...)(a+b+c...)n==+2ab++2ac+2bc+...
2
a
2
b
2
c
with a,b,c, etc. as positive digital rationals in a chosen base. The equation above implies a series of inequalities:
n>n->2ab+n---2ab>2(a+b)c+...
2
a
2
a
2
b
2
a
2
b
2
c
It’s already apparent from just looking at these inequalities that we can find each a, b, c, etc. by incrementing while the condition is satisfied. As we get more and more digits, we just need to keep the remainder on the left hand side as well as a running total on the right hand side. The code is not too bad:
In[]:=
SqrtInit[base_][n_]:=Association["Base"->base,"Total"->0,"Remainder"->n,"Output"->0,"Exponent"->Length[IntegerDigits[Floor[Sqrt[n]],base]]]
In[]:=
IterateSqrtDigit[state_]:=Module[{exp=state["Exponent"]-1,tot=state["Total"],max,next=0,den,rem},max=state["Remainder"];den=state["Base"]^exp;rem=next*den;While[max>rem^2+2tot*rem,rem=++next*den;];rem=rem-1*den;Association["Base"->state["Base"],"Remainder"->max-rem^2-2*tot*rem,"Output"->next-1,"Exponent"->exp,"Total"->tot+rem]]
To spot-check correctness, we can just look at the first few terms:
In[]:=
#Output&/@Rest[NestList[IterateSqrtDigit,SqrtInit[3][2],20]]First[RealDigits[Sqrt[2],3,20]]Union[%-%%]
Out[]=
{1,1,0,2,0,1,1,2,2,1,2,2,2,0,0,1,2,1,2,2}
Out[]=
{1,1,0,2,0,1,1,2,2,1,2,2,2,0,0,1,2,1,2,2}
Out[]=
{0}
Presently, we’ve left the neopixels running in the office, so instead of outputting to the peripheral, let’s use monitor to output directly to the notebook:
In[]:=
Module[{state=SqrtInit[3][2],digits=Association[],ind=0},CheckAbort[Monitor[While[True,state=IterateSqrtDigit[state];digits[++ind]=state["Output"];Pause[1/4];],Values[digits]],Values[digits]]]
Out[]=
{1,1,0,2,0,1,1,2,2,1,2,2,2,0,0,1,2,1,2,2}
Here we can sit and watch the approximation grow digits and real time, and when we abort the open loop, it finally outputs everything we’ve calculated so far, very nice. If we take out the pause, the sequence will grow much faster, and we can easily get 10,000 digits:
In[]:=
Module[{state=SqrtInit[3][2],digits=Association[],ind=0},CheckAbort[Monitor[While[True,state=IterateSqrtDigit[state];digits[++ind]=state["Output"];],Length@digits],Short@Values[digits]]]
Out[]//Short=
{1,1,0,2,0,1,1,2,2,1,2,2,2,0,0,1,2,1,2,2,10670,0,1,0,1,1,1,0,2,1,0,1,2,0,1,2,1,2,2,1}
It’s not the fastest calculation, but at least it’s not using too much memory. We can check the memory footprint fairly easily, and it looks to be growing only linearly
ListPlot[ByteCount/@NestList[IterateSqrtDigit,SqrtInit[3][2],10000]]
Out[]=
Since we’re not just copying someone else’s standard calculation of a transcendental number, it’s really not clear whether the algorithm we’ve just written actually classifies as a spigot. That uncertainty is okay, and honestly, the classification of an algorithm doesn’t matter as much as its behavior.
Dragons in Pi
Dragons in Pi
A few years ago, Danielle invited me to the “Find My PiPday” site. It’s a neat sideshow for the “me generation”, but instead of anthropocentrism, let’s try to find the unicode for dragon:
In[]:=
ToCharacterCode["龍"]
Out[]=
{40845}
By just giving the website a numeric input, it produces a date-formatted result that makes sense in one context but not the other:
Indeed, if we have a neopixel peripheral, we don’t really need a computer screen, much less a fancy graphic to attract the attention of prospective students. What we really need is a first-principles calculation that gets the same number 96,174 without needing to use a black-box web form.
So let’s go back to Beuker's Pi day lecture, and our followup demonstration. The basic idea is to take one of these approximation sequences, and iterate through increasingly precise rationals, while popping out digits whenever they sufficiently agree. The code is not too bad:
In[]:=
PiInit[base_]:=Association["Index"->2,"Annihilator"->{-4+4n,2-4n,-n},"Integers"->{1,-2,8},"Rationals"->{0,8,-24},"Approximation"->{0,4,3},"Output"->0,"Place"->Length[IntegerDigits[3,base]]];
In[]:=
IterateValues[ann_][state_]:=Module[{ind=state["Index"]+1,int=Rest[state["Integers"]],rat=Rest[state["Rationals"]],app=Rest[state["Approximation"]],den=Last[ann],next=state},AppendTo[int,ReplaceAll[-Most[ann].int/den,n->ind]];AppendTo[rat,ReplaceAll[-Most[ann].rat/den,n->ind]];AppendTo[app,Divide[-Last[rat],Last[int]]];next["Index"]=ind;next["Integers"]=int;next["Rationals"]=rat;next["Approximation"]=app;next]
It’s unproven whether the convergence criteria for re-iterating values will produce correct results for all bases, and we suspect it could error for large radix. However, for base = 10, the algorithm seems to work fine:
Now that we have a spigot-style streamer, let’s try to find the unicode for 龍:
Since we’re not using the AGM, it takes a noticeable amount of time to get even 10,000 digits. Even though this calculation seems to be lasting quite a long time, at 34,000 digits the rate of enumeration is still decently high, maybe slightly lower than 100 digits per second. From this observation, we could estimate the total search time will last for 30 minutes to an hour. Depending how rates slow down due to increasing memory footprint, the actual time-to-halt might end up closer to two or three hours. This is a good reason to take a meditation break, and return later to examine the result.
Indeed, if we just check the last values, we find 40845 starting at index 96,174:
The total time to calculate was about an hour, with an average iteration rate of about 1400 digits per minute:
We can also check that the convergence criteria never failed to produce a correct digit:
What’s import to realize is that our experiment is not about finding oddities in the digits of Pi. The purpose of this memo is just to make a neat display that scrolls through digits and find then find some digit-producing algorithms that pair nicely with the peripheral. Throughout the entire run to 96,000, iterations stayed well faster than one per quarter second, so I’d say our Pi spigot is a winner--exactly what we want to use for open-loop outputting.
Comparing with standards
Comparing with standards
In an article on “Unbounded spigot algorithms for the digits of Pi”, Jeremy Gibbons gives a concise but “deliberately obfuscated” spigot for calculating Pi digits. In Mathematica it can be written as:
The first few digits are produced as follows:
Removing bounds on recursion limit allows us to calculate up to 1000 digits and compare with our calculation from the previous section:
Despite being slower to 1,000, our new spigot algorithm has a smaller memory footprint. This gives us hope that refinement could improve our times significantly. Alas, Pi day has come and gone, so potential optimization will have to remain a task for some other day...