Lecture 6 :  Computational Experiments with Mathematica

For the 2020 UNCG Summer School in Computational Number Theory: Ergodic Theory with Applications to Continued Fractions
http://go.uncg.edu/ergodictheory
Written by Clifford Smyth:
http://www.uncg.edu/~cdsmyth/
cdsmyth@uncg.edu

Tips for Using Mathematica

Reference“An Elementary Introduction to the Wolfram Language, Second Edition”

For a basic yet comprehensive introduction to much of the functionality of Mathematica,see "An Elementary Introduction to the Wolfram Language, Second Edition" by Steven Wolfram at
 https://www.wolfram.com/language/elementary-introduction/2nd-ed/

Reference the Help Documentation

To search for commands and interactive tutorials on how to use them,select the "Help" option in the menu bar of the application.

Interrupting Evaluation and Saving Work Often

Save your notebook often.  It is easy to write code that unexpectedly hangs.
The keyboard commands “Command .” (in Mac) or “Ctrl .” will terminate a run.
But then Mathematica will try to save its often enormously complicated internal “state”, so this quitting can take far longer than you may be willing to wait. 
To instantly regain control, you will then have to use the menu bar options “Evaluation” then “Quit Kernel”. 
But then you will lose all computational results that weren’t saved to the visible screen in the notebook or to data files.

Examples:
The command
x=(long computation)
sets the variable x be the result of the the long computation and prints a copy of x to the notebook.  That print-out will stay in the notebook after quitting the kernel and be used to instantly reinitialize x to its previous value without going through the computation again.  It’s even better to have saved it to a data file.

On the other hand the command
x=(long computation);
sets x as before but the ; symbol at the end suppresses printing out a copy of a to the notebook.  The result will be lost after quitting the kernel.

Writing Fast Code

Use the very fast built-in code already available for  thousands of problems.

If you want to do something, spend your first five minutes searching Mathematica’s help documentation with its built in text searching capability instead of just plunging in.

It might very well already have been done.  Also, if it has been done, the documentation is generally good and provides you with many worked examples.

For example, searching for “continued fractions” quickly leads you to a tutorial on the built-in ContinuedFraction command.

You probably won’t be able to write code from scratch that can get the first 1000 continued fraction digits of Euler’s E number faster than .002 seconds.

Timing[
ContinuedFraction[E,1000]
]

Learn how to exploit Mathematica’s very fast list/structure manipulation and function mapping commands.

If you are familiar with procedural programming with traditional program flow commands: For, If/Then, While, etc. that is all available in Mathematica.

However, if you want to write code that is many times faster you need to learn and take advantage of Mathematica’s built in list manipulation and function mapping capabilities.

Here are two examples: getting a list of primes up to 1 million and then getting the gaps between those consecutive primes.

If you use list manipulation and function mapping commands instead of traditional program flow commands you can do these things 10 and 15 times faster.

Here’s an example of how to get the primes numbers up to 1000 traditionally using Mathematica’s built in PrimeQ testing function and For and If/Then commands.

primes={};
For[i=2,i≤1000,i=i+1,If[PrimeQ[i],primes=Append[primes,i]]]

primes

Getting the primes up to 1 million this way takes about 12 seconds.

Timing[
primes={};
For[i=2,i≤1000000,i=i+1,If[PrimeQ[i],primes=Append[primes,i]]];
primes
]

Getting the same list by first creating a list of the first 1 million positive integers and then selecting out only the primes using Mathematica’s Select function is almost 10 times faster.

Timing[
primes=Select[Range[1000000],PrimeQ[#]&]
]

Once this list of primes is created creating a list of the gaps between the consecutive primes up to 1 million takes .2 seconds if you do that procedurally.

Timing[
gaps={};
For[i=1,i≤Length[primes]-1,i=i+1,Append[gaps,primes[[i+1]]-primes[[i]]]];
gaps
]

Here’s you can split a list into overlapping sublists of length 2.

Partition[{1,4,9,16,25},2,1]

Here’s how you can map a function onto each element x on the list.  The particular function here being mapped onto x takes the second element of x and subtracts the first element of x.

Map[(#[[2]]-#[[1]])&,Partition[{1,4,9,16,25},2,1]]

Note that using Mathematica’s built in Partition and Map commands creates the same set of gaps about 15 times faster.

Timing[
gaps=Map[(#[[2]]-#[[1]])&,Partition[primes,2,1]]
]

Sample Code

The Base 10 Shift Map

Consider the base 10 digit shift map T:[0,1) -> [0,1) given by T(x) = (10 x) mod 1.
Here is the Mathematica implementation.

T[x_]:=10 x-Floor[10 x]

Here is a quick way to create an orbit of a function.

NestList[f,1/Pi,4]

Here are the first 10 iterates of the base 10 shift map applied to 1/Pi.

NestList[T,1/Pi,10]

Here is the 15 digit accuracy decimal expansions of those iterates.

N[NestList[T,1/Pi,10],15]

Here is the 25 digit accuracy decimal expansion of 1/Pi in which you can see those iterates.

N[1/Pi,25]

The Base b Shift Map.

Now we consider the base b shift map T_b: [0,1) -> [0,1) given by T_b(x) = (b x) mod 1.

T[b_,x_]:=b x - Floor[b x]

N[NestList[T[3,#]&,1/Pi,10]]

It is not clear from the decimal output that T_3 is really a shift map though.

Here are the first 20 digits of the base 3 decimal expansion of 1/Pi.  The first element of the output is the list of digits and the second element, -1, indicates that 1/Pi  = 3^{-1} (2.212100....) in base 3

RealDigits[1/Pi,3,20]

We’ll write our own version of RealDigits, that, when applied to x in [0,1), outputs the digits {a,b,c,...} of the base b expansion of x, x=.abc... even if it has an initial sequence of 0s

realdigits[x_,b_,n_]:=
(* Assumes that x is in [0,1) *)
Block[{digits,exponent,output},
{digits,exponent}=RealDigits[x,b,n];
output=Take[PadLeft[digits,n-exponent],n];
Return[output]
]

realdigits[1/Pi,3,15]

Now we can see that T_3 is a shift map.

Map[realdigits[#,3,5]&,N[NestList[T[3,#]&,1/Pi,10]]]

The Gauss Map or Continued Fraction Digit Shift Map and Orbits

Here is a straightforward implementation of the Gauss map on [0,1), to 30 decimal accuracy.

GaussMap[x_]:=N[1/x-Floor[1/x],20]

Due to round-off errors you quickly run into trouble in getting orbits of the GaussMap when it is implemented this way.

NestList[GaussMap,1/Pi,22]

Since the Gauss map is the shift map on the continued fraction expansion, the Mathematica ContinuedFraction and FromContinuedFraction commands can be used to get orbits of the Gauss Map to any desired accuracy.

Here are the first 30 continued fraction digits of 1/Pi.  The leading 0 is the integer part.

ContinuedFraction[1/Pi,30]

This is the rational approximant to 1/Pi corresponding to the preceding partial expansion.

FromContinuedFraction[%]

It can be seen to be a good approximation by looking at the corresponding decimal expansions.

N[%,40]

N[1/Pi,40]

If you want the first 20 shifts of the Gauss map of x in [0,1) with 10 continued fraction digits of accuracy in each you can do the following.
Take the first 30 continued fraction digits.
Drop the leading 0, the integer part.
Partition that list of digits into 20 sublists of length 10 each, overlapping by 1.
Insert the leading 0, the integer part, back at the front of those 10 digits expansions.
Then convert those expansions back to the corresponding rational approximants of the iterates of the Gauss map.

ContinuedFraction[1/Pi,30]

Drop[ContinuedFraction[1/Pi,30],1]

Partition[Drop[ContinuedFraction[1/Pi,30],1],10,1]

Map[Prepend[#,0]&,Partition[Drop[ContinuedFraction[1/Pi,30],1],10,1]]

Map[FromContinuedFraction,Map[Prepend[#,0]&,Partition[Drop[ContinuedFraction[1/Pi,30],1],10,1]]]

Map[N,%]

Here is that all coded up in one function.

GaussMapOrbit[x_,n_,d_]:=
(* GaussMapOrbit[x,n,d] will return x and the first n iterates of the Gauss map on x,each represented as a d-digit continued fraction approximant.  The code does the following.
a) Gets n + d -1 digits of the coninued fraction digts of x after the initial integer part of 0.
b)Partitions this into n windows of d digits each, each shifted by 1. These are the d-digit expansions of iterates 0,...,n-1 of the Gauss map.
c) Adds back the initial integer parts of 0s to the d digit expansions.
d) Transforms these expansions back to d-digit rational approximants of the iterates.
*)
Map[FromContinuedFraction,Map[Prepend[#,0]&,Partition[Drop[ContinuedFraction[x,n+d],1],d,1]]]

GaussMapOrbit[1/Pi,20,10]

Coding via Transformations versus Coding via Procedures

If you are used to coding things procedurally it can be quite a game-changer to see something coded up with list manipulations and function mappings. 

The paradigm is different.  You are now viewing code as a series of transformations applied to the input until it is transformed into the output.

There are some benefits to coding in this fashion.

Debugging and writing happen together and bugs can be caught more quickly.  You apply the transformations you think you need to apply to move from some trial input that leads to some known output and can check that you are getting the right results after each step.  Once you are satisfied each step is correct you can build the final code. 

The final code can be much shorter than a procedural equivalent.

The final code can be much faster than a procedural equivalent.

Many times, auxiliary variables whose only purpose is to carry out the transformations can be eliminated.

However, clarity can suffer greatly.  The high-level structure of what is being done can actually be more apparent from code written this way.  But that will usually only be the case if you record what you are doing in comments that are a part of the function.  GaussMapOrbit without the included comments would be very hard to parse if you had to return to it a long time after you had written it. 

The major advantages to this transformational approach to coding are the speed in which you can build up code and the speed of the code you write.

The Orbit of the Base b Shift Map

This command gives orbits of the base b shift map applied to x in [0,1) as rational number approximants.  See description of how the GaussMapOrbit function was coded to appreciate how this code works and how it might have been designed from the ground up.

ShiftMapOrbit[x_,b_,n_,d_]:=
(* Assumes x is in [0,1) and returns iterates 0 through n-1 of the base b shift map on x, up to d digit accuracy *)
Map[b^(-d)FromDigits[#,b]&,Partition[realdigits[x,b,n+d-1],d,1]]

ShiftMapOrbit[(1/3)^5,3,7,5]

Sample Explorations

Distribution of the Orbits of Shift Maps

Here’s a way to use Mathematica to visualize an orbit in [0,1).

ListPlot[{.1,.4,.7,.5,.9,.05}]

These are the first 10 000 iterates of the base 10 shift Map applied to 1/Pi

ListPlot[ShiftMapOrbit[1/Pi,10,10000,5]]

These are the first 10 000 iterates of the Gauss map applied to 1/Pi.

ListPlot[GaussMapOrbit[1/Pi,10000,20]]

Digit Frequency and Digit Correlation

If you divide [0,1) into 10 intervals, the jth iterate of the base 10 shift map applied to x landing in the interval [i/10,(i+1)/10) indicates that the jth decimal digit of x is i.
Plotting the corresponding histogram shows that digits in the base 10 expansion of 1/Pi are fairly evenly distributed.

Histogram[ShiftMapOrbit[1/Pi,10,10000,5],10]

This piece of code uses Mathematica’s generalized inner product function to turn the orbits of two maps into a single orbit of the product of the two maps.

Inner[List,{1,2,3,4},{-1,-2,-3,-4},List]

We use that idea to examine the orbit (T^k(1/Sqrt[2]), T^k(1/Sqrt[3])) where T is the base 10 shift map.

We plot the histogram of these pairs as they fall into the regions [i,(i+1)/10) \times [j/10,(j+1)/10).

In effect, we are recording the number of times the kth decimal digits of 1/Sqrt[2] and 1/Sqrt[3] are (i,j).

From the histogram we see that there does not appear to be any clear correlation between the decimal digits of these two irrational real numbers.

Histogram3D[Inner[List,ShiftMapOrbit[1/Sqrt[2],10,100000,5],ShiftMapOrbit[1/Sqrt[3],10,100000,5],List],10]

Even if they are irrational, when the two numbers are commensurable though something interesting happens.  This “correlation” between the shift maps you see in the histogram below can be easily proved though.

Histogram3D[Inner[List,ShiftMapOrbit[1/(5Sqrt[2]),10,100000,5],ShiftMapOrbit[1/(3Sqrt[2]),10,100000,5],List],100]

Below is a histogram of the first 100 000 pairs of continued fraction digits 2^{-1/3} and 3^{-1/3}.   There appears to be no visually obvious correlation between the digit sequences.

Histogram3D[Inner[List,ContinuedFraction[1/2^(1/3),100000],ContinuedFraction[1/3^(1/3),100000],List],{1,20,1}]