Dates everywhere in pi(e)!

[mtrott (Michael Trott)]

Some statistical and numerological musings about the occurrences of dates in the digits of pi

---

Notebook is attached at the end of this post.

The Pi Day celebration continues here at Wolfram. Yesterday, Stephen Wolfram wrote a blog post showing us how to get our own piece of pi at mypiday.com. And on Saturday 3.14.15 at 9:26:53.589, the celebration continues at SXSW with a countdown event to the pi day of the century moment.

For those interested in more fun with the dates in digits of Pi, I’ve used the Wolfram Language to do some analyses here. Please feel free to use this as a launchpad for your own explorations, or suggest other aspects for investigation.

---

In the last blog, Stephen discussed the unique position of this year’s pi-day within all pi-days of this century and gave various examples of the occurrence of dates in the (decimal) digits sequence of pi. In this blog, we will look at the statistics of the distribution of all possible dates/birthdays from the last 100 years within the (first ten million decimal) digits of pi. We will find that 99.998% of all digits occur in a date and that we find millions of dates within the first ten million digits of pi.

Here, we will concentrate on dates than can be described with maximally 6 digits. This means, we will be able to uniquely encode all dates between Saturday March 14 2015 and Monday March 15 1915 - a time range of 36,525 days.

We start with a graphical visualization of the topic at hand to set the mood.

piMask = ImagePad[ColorConvert[ImageCrop[
     Image[
      Rasterize[
       Style[TraditionalForm[Pi], 30, Bold, FontFamily -> "Times", 
        Background -> White], ImageSize -> 850]]], "Grayscale"], 80, 
   White];

dateBackground = ImageCrop[Image[Graphics[Table[
      Inset[
       Rotate[DatePlus[DateObject[{1915, 3, 14}], 
         RandomInteger[{1, 36525}]],
        2 Pi RandomReal[]], RandomReal[{-1, 1} 20, 2]], {2000}]], 
    ImageSize -> 900], {850, 850}];

Show[ImageCrop[ImageAdd[{dateBackground, piMask}]], ImageSize -> 400]

Get all dates for the last 100 years

So, this year’s pi day is, like every year, on March 14.

piDay = DateObject[{2015, 3, 14}]

Since the centurial pi day of the 20th century, 36,525 days have gone by.

piDay100YearsAgo = DateObject[{1915, 3, 14}]

DateDifference[piDay, piDay100YearsAgo]

Quantity[-36525, "Days"]

We generate a list of all the 36,525 dates on interest.

datesLast100Years = Table[DatePlus[piDay, d], {d, 0, -36524, -1}];
Length[datesLast100Years]

36525

{First[datesLast100Years], Last[datesLast100Years]}

For later use, we define a function dateNumber that for a given date returns the number of the date; we start counting on Mar 15 1915.

MapIndexed[(dateNumber[#[[1]]] = #2[[1]]) &, datesLast100Years];

For the months January to September, we do not require a date to be written as 01, 02, …, 12 for the months and similar for days 01 to 09 we allow a date specified using just 1 … 9. So, for some dates, there are multiple digits sequences that represent it. The function makeDateTuples generates all tuples of single-digit integers that represent a date. One could use slightly different conventions and minimal chnages of the code and always enforce short dates or always enforce zeros. With the optional inclusion of zeros for days and month, we get more possible matches and a richer result, so we will use these in the following. (And, if you prefer a European date format in the form day-month-year, then some larger adjustments have to me made to the definition of makeDateTuples.)

makeDateTuples[DateObject[{year_, month_, day_}, ___]] := 
 Module[{monthDay, y = Take[IntegerDigits[year], -2] },
  monthDay =
   Which[month < 10 && day < 10, {{month, day}, {0, month, day},{month, 0, day}, {0, month, 0, day}},                
    month >= 10 && day < 10, {{month, day}, {month, 0, day}},
    month < 10 && day >= 10, {{month, day}, {0, month, day}},
    month >= 10 && day >= 10, {{month, day}}];
  Flatten[Append[#, y]] & /@ Map[IntegerDigits , monthDay, {-1}]]

Some examples with four, two, and one representation:

makeDateTuples[DateObject[{1955, 5, 5}]]

{{5, 5, 5, 5}, {0, 5, 5, 5, 5}, {5, 0, 5, 5, 5}, {0, 5, 0, 5, 5, 5}}

makeDateTuples[DateObject[{1964, 6, 10}, CalendarType -> "Gregorian"]]

{{6, 1, 0, 6, 4}, {0, 6, 1, 0, 6, 4}}

makeDateTuples[DateObject[{1954, 12, 1}]]

`{{1, 2, 1, 5, 4}, {1, 2, 0, 1, 5, 4}}’

makeDateTuples[DateObject[{1981, 11, 10}]]

{{1, 1, 1, 0, 8, 1}}

The next plot shows which days from the last year are representable with 4, 5, and 6 digits.The first nine days of the months January to September just need four or five digits to be represented and the last days of October, November and December need six.

DateListPlot[
 Flatten[Function[d, {d, Length[#]} & /@ makeDateTuples[d]] /@ 
   Table[ DatePlus[DateObject[{2014, 1, 1}], d], {d, 0, 365, 1}], 1], 
 Joined -> False, PlotStyle -> PointSize[0.002], PlotRange -> {0, 7}]

For a fast (constant time), repeated recognition of a tuple as a date, we define two functions dateQ and dateOf. dateOf gives a normalized form of a date digits sequence. We start with generating pairs of tuples and their date interpretations.

tuplesAndDates = {#[[1, 1]], Last /@ #} & /@ Split[
    Sort[Flatten[
      Function[d, {#, d} & /@ makeDateTuples[d]] /@ 
       Take[datesLast100Years, All], 1]], #1[[1]] === #2[[1]] &];

Here are some examples.

RandomSample[tuplesAndDates, 4]

Most (77,350) tuples can be uniquely interpreted as dates, some (2,700) have two possible date interpretations.

Tally[Length[Last[#]] & /@ tuplesAndDates]

{{1, 77350}, {2, 2700}}

Here are some of the digits sequences with two date interpretations.

RandomSample[#, 4] &@Cases[tuplesAndDates, {_, {_, _}}]

Function[{tuple, dates}, dateQ[tuple] = True; 
   datesOf[tuple] = (First /@ dates)] @@@ tuplesAndDates;

dateQ[_] := False

Here are the two date interpretations of the sequence {1,2,1,5,4} as Jan 21 1954 or as Dec 1 1954 recovered by using the function datesOf.

DateObject /@ datesOf[{1, 2, 1, 5, 4}]

These are the counts for the 4, 5, and 6–digits representations of dates.

{#[[1, 1]], Length[Union[Last /@ #]]} & /@ 
 Split[Sort[{ Length[#1[[1, 1]]], #2} & @@@ 
    Take[DownValues[datesOf], All]], #1[[1]] === #2[[1]] &]

{{4, 8100}, {5, 28225}, {6, 36525}}

And the counts for the number of definitions set up for the function datesOf.

{#[[1, 1]], Length[#]} & /@ 
 Split[Sort[{ Length[#1[[1, 1]]], #2} & @@@ 
    Take[DownValues[datesOf], All]], #1[[1]] === #2[[1]] &]

{{4, 8100}, {5, 35425}, {6, 36525}}

Find all dates in the digits of pi

For all further calculations, we will use the first ten million decimal digits of pi. We allow for an easy substitution of pi by another constant.

constant =  Pi  (* \[ExponentialE] *) (* ChampernowneNumber[10] *);
base = 10;
constantDigitSequence = RealDigits[constant, base, 10^7][[1]];

Instead of using the full digit sequence as a string, we will use the digit sequence split into (overlapping) tuples. Then we can independently and quickly operate onto each tuple. And we index each tuple with the index representing the digit number. E.g.

Transpose[{#, Range[Length[#]]}] &[
 Partition[ Take[constantDigitSequence, 20], 4, 1]]

{{{3, 1, 4, 1}, 1}, {{1, 4, 1, 5}, 2}, {{4, 1, 5, 9}, 3}, {{1, 5, 9, 2}, 4}, {{5, 9, 2, 6}, 5}, {{9, 2, 6, 5}, 6}, {{2, 6, 5, 3}, 7}, {{6, 5, 3, 5}, 8}, {{5, 3, 5, 8}, 9}, {{3, 5, 8, 9}, 10}, {{5, 8, 9, 7}, 11}, {{8, 9, 7, 9}, 12}, {{9, 7, 9, 3}, 13}, {{7, 9, 3, 2}, 14}, {{9, 3, 2, 3}, 15}, {{3, 2, 3, 8}, 16}, {{2, 3, 8, 4}, 17}}

Using the above-defined dateQ and dateOf functions we can now quickly find all digits sequences that have a date interpretation.

getDates[{tupleList_, posList_}] :=  
 Module[{cL = dateQ /@ tupleList},
   Transpose[{Pick[tupleList, cL], Pick[posList, cL]}]]

getDateLists[k_] := 
 With[{l = Partition[ Take[constantDigitSequence, All], k, 1], 
   p = Partition[#, 10^6, 10^6, 1, {}] &}, 
  Flatten[getDates /@  Transpose[{p[l], p[Range[Length[l]]]}], 1]]

makeDate[{tuple_, pos_}] := {#, pos, tuple} & /@ datesOf[tuple]

fourDigitDates = Flatten[makeDate /@ getDateLists[4], 1];

fiveDigitDates = Flatten[makeDate /@ getDateLists[5], 1];

sixDigitDates = Flatten[makeDate /@ getDateLists[6], 1];

Here are some of the date interpretations found. Each sublist is of the form {date, startingDigit, digitSequenceRepresentingTheDate}.

RandomSample[fourDigitDates, 6]

RandomSample[fiveDigitDates, 6]

RandomSample[sixDigitDates, 6]

We have found about 8.1 million dates represented as 4-digits, about
3.8 million dates as 5 digits and about 365 thousand dates represented as six digits totaling more than 12 million dates.

Length /@ {fourDigitDates, fiveDigitDates, sixDigitDates}

{8100779, 3811415, 365491}

Total[%]

12277685

Note that we could have used string-processing functions (especially StringPosition) to find the positions of the date sequences. And, of course, we would have obtained the same result.

constantDigitString = StringJoin[ToString /@ constantDigitSequence];

(strCaces = StringPosition[constantDigitString, "82959"]) // Short[#, 4] &

(#2 & @@@ Cases[fiveDigitDates, {{1959, 8, 29}, __}]) === (#1 & @@@ strCaces)

True

While the use of StringPosition is a good approach to deal with single date, dealing with all 35 thousand sequences would have taken much longer.

Quantity[First[Timing[StringPosition[constantDigitString, "82959"]]], 
   "Seconds"] Length[tuplesAndDates] // UnitConvert[#, "Minutes"] &

Quantity[First[Timing[StringPosition[constantDigitString, "82959"]]], 
   "Seconds"] Length[tuplesAndDates] // UnitConvert[#, "Minutes"] &

Quantity[190.48, "Minutes"]

We pause a moment and have a look at the counts found for the four tuples. Out of the 10,000 possible four tuples the 8,100 used appear each in average (1/10)^4*10^7=10^4 times based on the randomness of the digits of pi. And approximately, we expect a standard deviation of approximately 1000^(1/2) ~ 31.6. Some quick calculations and a plot confirm these numbers.

fourTupleCounts = Tally[fourDigitDates[[All, 3]]];
{Mean[Last /@ fourTupleCounts], StandardDeviation[Last /@ fourTupleCounts]} // N

{1000.1, 31.535}

The histogram of the counts shows the expected Bell curve.

Histogram[#2 & @@@ fourTupleCounts, 60]

And the following graphic shows how often each of the 4-tuples that represent dates were found in the 10 million decimal digits. We enumerate the 4-tuples by concatenating the digits; as a result we see ‘empty’ vertical stripes in the region where no 4-tuples are represented by dates.

ListPlot[Tooltip[{FromDigits[#1], #2}, {#1, #2}] & @@@ 
  fourTupleCounts, PlotRange -> All, AxesLabel -> {"4-tuple", "counts"}]

Now we continue to process the found date positions. We group the results into sublists of identical dates.

dateGroups = 
  Split[Sort[Join[fourDigitDates, fiveDigitDates, sixDigitDates]], #1[[1]] === #2[[1]] &];
Clear[fourDigitDates, fiveDigitDates, sixDigitDates]

Every date does indeed occur in the first 10 million digits, meaning we have 36,525 different dates found. (We will see later that we did not calculate much more digits than needed.)

Length[dateGroups]

36525

Here is how a typical member of dateGroups looks like.

dateGroups[[27835]] // Short[#, 4] &

Statistics of all dates

Now let’s do some statistics on the found dates. Here are the number of occurrences of each starting date positions in the first 10 million digits of pi.Interestingly, many dates appear hundreds of times. The periodically occurring vertical stripes result from the October-November-December months.

DateListPlot[{#[[1, 1]], Length[#]} & /@ Take[dateGroups, All], 
 PlotRange -> All, Filling -> Axis, Joined -> False, PlotRange -> All,
  PlotStyle -> PointSize[0.001], FillingStyle -> Directive[Thickness[0.001], Opacity[0.2]]]

For easier readability, we format the triples {date, startingPosition, dateDigitSequence} in a customized manner.

formatOccurrence[{date_, startingPosition_, dateDigitSequence_}] := 
 With[{s = Text[Style[#, Gray]] &}, 
  Framed[Row[{DateObject[date], s[" at "], startingPosition , 
     s[" in "], constant, s[" as "], 
     StringJoin[
      Riffle[ToString /@ dateDigitSequence, 
       "\[NegativeThinSpace]\[CenterDot]\[NegativeThinSpace]"]]}], 
   RoundingRadius -> 5, Background -> GrayLevel[0.96], 
   FrameMargins -> 1]]

(* list of data *)

formatOccurrence[l : {{_List, _, _List} ..}] := formatOccurrence /@ l

dateGroups[[16238, 1]]

formatOccurrence[%]

Most occurring dates in the first ten million digits of pi is Aug 6 1939 - it occurs 1,362 times.

SortBy[dateGroups, Length][[-1]] // formatOccurrence // Short[#, 0.4] &

Least occurring dates in the first ten million digits of pi. These three birthdays occur only once in the first 10 million digits:

Select[dateGroups, Length[#] === 1 &] // (formatOccurrence /@ #) &

And all of these birthdays occur only twice in the first 10 million digits:

Select[dateGroups, Length[#] === 2 &] // (formatOccurrence /@ #) & // Short[#, 0.3] &

Here is the distribution of the number of the date occurrences. The three peaks corresponding to the 6, 5, and 4 digit date representations (from left to right) are clearly distinct. The dates that are represented by 6-tuples occur each only a very few times and, as we have already seen above, appear in average about 1,200 times.

Histogram[Length /@ dateGroups, {10}, AxesLabel -> {"number of dates", "counts"}]

We can also accumulate by year and display the date interpretations per year (the smaller values at the beginning and end come from the truncation of the dates to ensure uniqueness.) The distribution is nearly uniform.

DateListPlot[{{#1, 1, 1}, #2} & @@@ (Tally[
    Flatten[{#1[[1]]} & @@@ Flatten[Take[dateGroups, All], 1]]]), 
 Filling -> Axis, Joined -> False, 
 PlotRange -> {{{1914, 1, 1}, {2016, 1, 1}}, {120, 125} 1000}]

Let’s have a look at the dates with some ‘neat’ date sequences and how often they occur. As the results in dateGroups are sorted by date, we can easily access a given date. When does the date 11-11-11 occur?

dateGroups[[36525 - dateNumber[{2011, 11, 11}] + 1]] // formatOccurrence

And where does the date 1-23-45 occur?

dateGroups[[36525 - dateNumber[{1945, 1, 23}] + 1]] // formatOccurrence // Short[#, 0.3] &

No date starts on its ‘own’ position (meaning there is no example such as Jan 1, 1945 (1-1-4-5) in position 1145).

Select[Flatten[Take[dateGroups, All], 1], FromDigits[#[[3]]] === #[[2]] &] `{}`

But one palindromic case exists: March 3, 1985 (3.3.8.5), which occurs at palindromic position 5833.

Select[Flatten[Take[dateGroups, All], 1], 
  FromDigits[Reverse@#[[3]]] === #[[2]] &] // formatOccurrence 

A very special date is Jan 9 1936: 1.9.3.6 appears at the position of the 1936th prime, 16747.

Select[Flatten[Take[dateGroups, All], 1], Prime[FromDigits[#[[3]]]] === #[[2]] &] // formatOccurrence 

Prime[1936]

16747

Let’s see what anniversaries happened on this day in history.

WolframAlpha["what happened Jan 9 1936", 
 IncludePods -> "NotableEventForDate", AppearanceElements -> {"Pods"},
  PodStates -> {"NotableEventForDate__More"}]

While no date appeared at its ‘own’ position, if we slightly relax this condition and search for all dates that overlap with its digits positions, we do find some dates.

Select[Flatten[Take[dateGroups, All], 1], 
  MemberQ[ #[[2]] + Range[0, Length[#[[3]]] - 1], 
    FromDigits[#[[3]]]] &] // formatOccurrence 

And at more than 100 positions within the first 10 million digits of pi we find the famous pi starting sequence 3,1,4,5,9 again.

({#2, Row[{StringJoin[ToString /@ Take[constantDigitSequence, {#2, #2 + Length[#3] - 1}]], 
       Style[StringJoin[ToString /@ Take[constantDigitSequence, {#2 + Length[#3], #2 + 
             Length[#3] + 2}]], Gray]}]} & @@@ 
   Cases[dateGroups[[-1]], {__, {3, 1, 4, 1, 5}}]) // Short[#, 4] &

Within the digits of pi we do not just find birthday, dates, but also the physical constants days, such as the [HBar]-day last time celebrated as the centurial instance on October 5 1945.

UnitConvert[Quantity[1, "ReducedPlanckConstant"]]

Position[dateGroups, {{{1945, 10, 5}, __} ..}, {1}]

{{11163}}

Here are the positions of the matching date sequences shown.

ListLogLinearPlot[MapIndexed[{#1[[2]], #2[[1]]} &, dateGroups[[11163]]]]

And here is an attempt to visualize the appearance of all dates. In the date-digit-plane, we place a point at the beginning of each date interpretation. We use a logarithmic scale for the digit position and as a result the number of points is much larger in the upper part of the graphic.

Graphics[{PointSize[0.001], Opacity[0.4], 
  Point[Developer`ToPackedArray@
    Flatten[MapIndexed[Function[{d, p}, 
       MapIndexed[({1. p[[1]], Log10[1. #[[2]]]}) &, d]], 
      Take[dateGroups, All]], 1]]},
 AspectRatio -> 1, Frame -> True, 
 FrameLabel -> {"dates since Mar 15 1915", "log10 digit position"}]

For the dates that appear early on in the digits sequence, the finite extension of the date over the digits can be visualized too. A date extends over four to six digits in the digit sequence. The next graphic shows all digits of all dates that start within the first 10,000 digits.

Graphics[{Thickness[0.001], Line[Flatten[MapIndexed[
     Function[{d, p}, 1. {{p[[1]], Log10[#2]}, {p[[1]], 
           Log10[#2 + Length[#3]]}} & @@@  
       Select[d, #[[2]] < 50000 &]], Take[dateGroups, All ]], 1]]},
 AspectRatio -> 1, Frame -> True, 
 FrameLabel -> {"dates since Mar 15 1915", "log10 digit position"}]

After coarse-graining, the distribution is quite uniform.

roundedDateOccurrences = Tally[ Flatten[ MapIndexed[Function[{d, p},
      MapIndexed[({Round[1. p[[1]], 100], 
          Ceiling[Log10[1. #[[2]]], 0.02]}) &, d]], Take[dateGroups, All]], 1]];

Histogram3D[
 WeightedData[1. First /@ roundedDateOccurrences, 
  1. Log10[Last /@ roundedDateOccurrences]],
 AxesLabel -> {"dates since Mar 15 1915", "log10 digit position", 
   "log10 birthday\nfrequency"},
 Ticks -> {Automatic, Automatic, False}]

So far we have taken a date an looked where this date starts in the digit sequence of pi. Now let’s look from the reverse direction: At the nth digit of [Pi], which dates are used, as a starting or following digits, to represent this date. To find the total counts of dates for each digits, we loop over the dates and accumulate the counts for each digit.

addDay[{date_, p_, l_}] := Do[lT[[p + j]]++, {j, 0, Length[l] - 1}]

lT = Table[0, {Length[constantDigitSequence]}];
(addDay /@ #) & /@ dateGroups;

Maximal 20 dates occur at a given digit.

{Min[lT], Max[lT], Mean[1. lT], StandardDeviation[1. lT]}

{0, 20, 5.36531, 1.45252}

Here are two intervals of two hundred digits each. We see that most digits are used in a date interpretation.

{ListPlot[Take[lT, {1, 200}], Filling -> Axis, ImageSize -> 240],
 ListPlot[Take[lT, 7777777 + {1, 200}], Filling -> Axis, ImageSize -> 240]}

Above, we noted that we have about 12 million dates in the digit sequence under consideration. The digit sequence that we used was only 10 million digits long and each date needs about 5 digits. This means the dates need about 60 million digits. This means that many of the 10 million digits must be shared and used in average about five times. Only 2,005 out of the first ten million digits are not used in any of the date interpretations, meaning that 99.98% of all digits are used for date interpretations (not all as starting positions).

{Count[lT, 0], Count[lT, 1], Count[lT, 20]}

{2005, 15424, 1}

And here is the histogram of the distribution of the number of dates present at a certain digit. The back-of-the-envelope number of an average of 6 dates per digits is clearly visible.

Histogram[lT]

The 2,005 positions that are not used are approximately uniformly distributed among the first ten million digits.

Histogram[Flatten[Position[lT, 0]], 200]

At a given digit, dates from various years overlap. The next graphic shows the range from the earliest to the latest year as a function of the digit position.

minMaxYears = {{#[[1, 1]] , Min[#[[2, 1]] & /@ #]} ,
     {#[[1, 1]] , Max[#[[2, 1]] & /@ #]} } & /@ (Union /@ Split[
      Sort[
       Flatten[Function[{d, p, l}, 
          Table[{p + j, d}, {j, 0, Length[l] - 1}]] @@@ 
         Flatten[Select[#, #[[2]] < 1000 &] & /@ dateGroups, 1], 
        1]], #1[[1]] === #2[[1]] &]);

ListLinePlot[Transpose[minMaxYears], Filling -> {1 -> {2}}, 
 PlotStyle -> Directive[Thickness[Small], Blue], AxesLabel -> {"digit", "year"}]

As the high coverage seems pretty unexpected, we select a random digits position and select all dates that use this digit.

selectOverlappingDates[pos0_] := 
 Flatten[Select[#, #[[2]] <= pos0 <= #[[2]] + Length[#[[3]]] - 1 &] & /@ dateGroups, 1]

datesWithPos0 = selectOverlappingDates[923456]

{{{1967, 3, 8}, 923454, {3, 8, 6, 7}}, {{1979, 8, 6}, 923455, {8, 6, 7, 9}}, {{1986, 1, 3}, 923453, {1, 3, 8, 6}}, {{1986, 2, 13}, 923452, {2, 1, 3, 8, 6}}, {{1986, 12, 13}, 923451, {1, 2, 1, 3, 8, 6}}, {{1991, 6, 7}, 923456, {6, 7, 9, 1}}}

And here is a visualization of the overlap of the dates.

overlappingDatePlot[pos0_, data_] :=
  With[{\[Lambda] = Length[data]},
  NumberLinePlot [Table[#2 + j, {j, 0, Length[#3]}] & @@@ data ,
   PlotLegends -> (formatOccurrence /@ data), Frame -> {True, True},
   Epilog -> {Table[
      Text[constantDigitSequence[[
        pos0 + j]], {pos0 + j, \[Lambda] + 1}], {j, -4, 4}],
     Lighter[Gray], Line[{{pos0, 0}, {pos0, \[Lambda] + 2}}]}, 
   ImagePadding -> 25,
   PlotRange -> {{pos0 - 6, pos0 + 6}, {-1, 8}}]];

overlappingDatePlot[923456, datesWithPos0]

The most-used digit is the 1 at position 2,645,274: 20 possible date interpretations meet at it.

Position[lT, 20]

{{2645274}}

Here are the digits in its neighborhood and the possible date interpretations.

Grid[Transpose[
  Prepend[Table[{2645274 + j, 
     constantDigitSequence[[2645274 + j]]}, {j, -3, 3}],
                  Text[Style[#, Gray]] & /@ {"position", "digit"}]], Dividers -> All]

datesWithPos0 = selectOverlappingDates[2645274];
formatOccurrence /@ datesWithPos0 

If we plot the years starting at a given digits for a larger amount of digits (say the first ten thousand), then we see the relatively dense coverage of date interpretations in the digits-date-plane.

digitsAndYears = {#[[1, 1]], Rest /@ #} & /@ 
   Split[Sort[{#2, #1[[1]], 
        Length[#3]} & @@@ Flatten[Select[#, #[[2]] <= 10000 &] & /@ 
dateGroups, 1]], #1[[1]] === #2[[1]] &];

Graphics[{PointSize[0.003], 
  Function[{d, da}, Point[ {d, #1} & @@@ da] ] @@@ digitsAndYears},
 Frame -> True, AspectRatio -> 1, FrameLabel -> {"digit", "year"}]

Let us now build a graph of dates that are ‘connected’. We will consider two dates connected if the two dates share a certain digit of the digit sequence (not necessarily the starting digit of a date).

digitsAndDates = Union /@ Split[
    Sort[Flatten[
      Function[{d, p, l}, Table[{p + j, d}, {j, 0, Length[l] - 1}]] @@@
        Flatten[Select[#, #[[2]] < 1000 &] & /@ dateGroups, 1], 
      1]], #1[[1]] === #2[[1]] &];

Graph[UndirectedEdge[DateObject[#1], DateObject[#2]] & @@@ 
  Union[Flatten[Subsets[Union[Last /@ #], {2}] & /@ digitsAndDates, 
    1]], VertexLabels -> Placed["Name", Tooltip]]

Here is the same as graph for the first 600 digits with communities emphasized.

CommunityGraphPlot@
 Graph[UndirectedEdge[DateObject[#1], DateObject[#2]] & @@@ 
   Union[Flatten[
     Subsets[Union[Last /@ #], {2}] & /@ Take[ digitsAndDates, 600], 
     1]], VertexLabels -> Placed["Name", Tooltip]]

We continue with calculating the mean distance between two occurrences of the same date.

meanBirthdaySpacings = {#[[1, 1]], 
     Mean[-1. Subtract @@@ Partition[#2 & @@@ #, 2, 1]]} & /@ 
   Take[Select[dateGroups, Length[#] > 1 &], All];

 DateListLogPlot[meanBirthdaySpacings, PlotRange -> All, 
 Filling -> Axis, Joined -> False, PlotRange -> All]

Grid[Prepend[{DateObject[#1], Quantity[#2, "Days"]} & @@@ 
   Take[SortBy[meanBirthdaySpacings, Last], 5],
  Text[Style[#, Gray]] & /@ {"birthday", "mean spacing"}],  Dividers -> Center]

The first occurrences of dates

Most interesting are the first occurrences of dates. So, lets extract these. We will work with two versions, one sorted by the date (the list firstOccurrences) they represent and one sorted by the starting position (the list firstOccurrencesSortedByOccurrence) in the digits of [Pi].

firstOccurrences = First /@ dateGroups;
firstOccurrencesSortedByOccurrence = SortBy[firstOccurrences, #[[2]] &];

Here are the possible date interpretations that start within the first 10 digits of [Pi].

Select[firstOccurrencesSortedByOccurrence, #[[2]] < 10 &] // formatOccurrence 

And here are the other extremes. The dates that appear deepest into the digits expansion.

Take[firstOccurrencesSortedByOccurrence, -12] // formatOccurrence 

We see that Wed Nov 23 1960 starts only at position 9,982,546(=2 7 713039) - so by starting with the first 10 million digits we were a bit lucky to catch it. A quick direct check of this record-setting date.

Take[constantDigitSequence, 9982546 + {0, 5}]

{1, 1, 2, 3, 6, 0}

So, who are (well-known) the lucky people associated with this number through their birthday?

WolframAlpha["famous people born Nov 23 1960", 
 IncludePods -> "BasicInformation:PeopleData", 
 AppearanceElements -> {"Pods"}, 
 TimeConstraint -> {30, Automatic, Automatic, Automatic}, 
 PodStates -> {"BasicInformation:PeopleData__Show full dates"}]

And, what were the moon phases on the top dozen out-pi-landish dates.

Grid[{formatOccurrence[#], 
    Show[ MoonPhase[#[[1]], "Icon"], ImageSize -> 32]} & /@ 
  Take[firstOccurrencesSortedByOccurrence, -12]]

And while Wed Nov 23 1960 is furthest out in the decimal digit sequence, the last prime date in the list is Oct 22 1995.

Select[Reverse[firstOccurrencesSortedByOccurrence], PrimeQ[#[[2]]] &, 1]

{{{1995, 10, 22}, 6901859, {1, 0, 2, 2, 9, 5}}}

In general, less than 10% of all first date appearances are prime.

Quantity[100. Count[#2 & @@@ firstOccurrencesSortedByOccurrence, 
    PrimeQ]/Length[firstOccurrencesSortedByOccurrence], "Percent"]

Quantity[0., "Percent"]

Here are the first-occurring date positions for the last few years. The burst in October, November, December of each year are caused by the need of five or six consecutive digits while January to September can be encoded with fewer digits of we skip optional zeros.

DateListPlot[Most /@ Take[firstOccurrences, -2000], PlotRange -> All, 
 Filling -> Axis, Joined -> False, AxesOrigin -> {DateObject[{1915, 3, 14}], 0}]

If we include all dates, we get, of course, a much denser filled graphic.

DateListPlot[Most /@ firstOccurrences, PlotRange -> All, 
 Filling -> Axis, Joined -> False, PlotRange -> All]

A logarithmic vertical axis shows that most dates occur between the thousands and millionths digit.

DateListLogPlot[Most /@ firstOccurrences, PlotRange -> All,  Joined -> False, PlotRange -> All]

To get a more intuitive understanding how overall uniformity and local randomness in the digits sequence (and as a result in the dates) we make a Voronoi tessellation of the day-digit plane based on points at the first occurrence of a date. The decreasing density for increasing digits results from the fact that we only take first date occurrences into account.

points = {QuantityMagnitude[
      DateDifference[DateObject[{1915, 3, 15}], DateObject[#]]], #2} & @@@  
   Take[firstOccurrencesSortedByOccurrence, All];

ListContourPlot[
  Append[#, RandomReal[]] & /@ Select[points, #[[2]] < 50000 &],  InterpolationOrder -> 0,
  FrameLabel -> {"days since Mar 15 1915", "digit"}] /. 
 p_Polygon :> {EdgeForm[Directive[Black, Thickness[0.001]]],   FaceForm[], p}

Easter Sunday positions are a good date to visualize as the vary in date over the years.

EasterDates = Flatten[{#1, StringSplit[#2]}] & @@@ Select[Partition[ 
       Flatten[ Drop[#, 2] & /@ Cases[
          
          Import["https://www.census.gov/srd/www/genhol/easter500.html", "Data"],
          {_String?(StringMatchQ[#, 
               " Date of Easter *"] &), __}, \[Infinity]]], 2], 
      1915 <= #[[1]] <= 2014 &] /.
    {"March" -> 3,   "April" -> 4} /. s_String :> ToExpression[s];

DateListLogPlot[
 Most /@ Sort[Cases[firstOccurrencesSortedByOccurrence, 
    Alternatives @@ ({#, __} & /@ EasterDates)]],
 Joined -> False, Filling -> Axis, PlotRange -> All]

The mean first occurrence as function of number of digits needed to specify a date depends, of course, on the number of digits needed to encode a date.

Grid[Prepend[{#[[1, 1]], 1. Mean[Last /@ #]} & /@

   Split[Sort[{ Length[#3], #2} & @@@ 
      firstOccurrencesSortedByOccurrence], #1[[1]] === #2[[1]] &],
  Text[Style[#, Gray]] & /@ {"n-digit date", "mean position"}],  Dividers -> Center]

The mean occurrence is at 239,083, but the due to the outliers at a few million of digits, the standard deviation is much larger.

{Mean[#2 & @@@ firstOccurrencesSortedByOccurrence],
  StandardDeviation[#2 & @@@ firstOccurrencesSortedByOccurrence]} // N

{239083., 557576.}

Here are the first occurrences of the ‘nice’ dates that are formed by repetition of a single digit.

jjDates = RotateLeft[
  Table[Cases[
     firstOccurrencesSortedByOccurrence, {{ 
       If[j == 1, 2000, 1900] + 10 j + j, j, j}, __}][[1]], {j, 9}]]

{{{1922, 2, 2}, 4903, {2, 2, 2, 2}}, {{1933, 3, 3}, 28468, {3, 3, 3, 3}}, {{1944, 4, 4}, 54526, {4, 4, 4, 4}}, {{1955, 5, 5}, 20981, {5, 0, 5, 5, 5}}, {{1966, 6, 6}, 21881, {6, 6, 6, 6}}, {{1977, 7, 7}, 1590, {7, 7, 7, 7}}, {{1988, 8, 8}, 4752, {8, 8, 8, 8}}, {{1999, 9, 9}, 763, {9, 9, 9, 9}}, {{2011, 1, 1}, 12701, {1, 1, 1, 1}}}

DateListPlot[{{#1, #2}} & @@@ jjDates, Joined -> False,
 PlotRange -> All, PlotLegends -> (DateObject[#1] & @@@ jjDates)]

The detailed distribution of the number of occurrences of first dates has a the largest density within the first few ten thousand digits.

{Histogram[#2 & @@@ firstOccurrencesSortedByOccurrence, 100,  ImageSize -> 280],
 Histogram[ Select[#2 & @@@ firstOccurrencesSortedByOccurrence, # < 250000 &],   100, ImageSize -> 280]}

A logarithmic axis shows the distribution much better, but because of the increasing bin sizes, the maximum has to be interpreted with care.

Histogram[#2 & @@@ firstOccurrencesSortedByOccurrence, {"Log",  80}, "Probability"]

And here is the cumulative distribution of the dates as a function of the digits positions. We see that the first 1% of the ten million digits covers already 60% of all dates.

ListLogLinearPlot[
 MapIndexed[{#1[[2]], 100.  #2[[1]]/36525} &,   firstOccurrencesSortedByOccurrence],
 AxesLabel -> { "digit", "date coverage (%)"}, Joined -> True,  GridLines -> Automatic]

More (slightly more) birthdays start at even positions than at odd positions.

Tally[Mod[#2, 2] & @@@ firstOccurrencesSortedByOccurrence]

{{1, 18136}, {0, 18389}}

We could do the same with mod 3, mod 4, [AliasDelimiter] … .The higher congruences, all considered again mod 2 are shown in the following image.

ArrayPlot[Table[ Mod[#2, 2] & @@@ Tally[Mod[#2, m] & @@@ 
     firstOccurrencesSortedByOccurrence], {m, 60}] ]

The actual number of first occurrences per year fluctuates around the mean value.

DateListPlot[{DateObject[{#[[1, 1]]}], Mean[Last /@ #]} & /@ 
  Split[Sort[{#1[[1]], #2 } & @@@ 
     Take[firstOccurrencesSortedByOccurrence, All]], #1[[1]] === #2[[ 1]] &],
 AxesOrigin -> {DateObject[{1915, 0, 0}], 0}]

The mean number of first date interpretations sorted by month clearly show the difference between the one-digit months and the two-digit months.

DateListPlot[{ {2015, #[[1, 1]]}, Mean[Last /@ #]} & /@ 
  Split[Sort[{#1[[2]], #2 } & @@@ 
     Take[firstOccurrencesSortedByOccurrence, All]], #1[[1]] === #2[[1]] &],
 Joined -> False, Filling -> Axis]

The mean number by day of the month (ranging from 1 to 31) is, in average, a slowly increasing function.

ListPlot[{ #[[1, 1]], Mean[Last /@ #]} & /@ 
  Split[Sort[{#1[[3]], #2 } & @@@ 
     Take[firstOccurrencesSortedByOccurrence, All]], #1[[1]] === #2[[1]] &],
 Joined -> False, Filling -> Axis, AxesOrigin -> {0, 0}]

Finally, here are the mean occurrence by weekday. Most first date occurrences happen for dates that are Wednesdays.

Grid[{ #[[1, 1]], Round@ Mean[Last /@ #]} & /@ 
  Split[Sort[{DayName[#1], #2} & @@@ 
     Take[firstOccurrencesSortedByOccurrence, All]], #1[[1]] === #2[[1]] &],
 Dividers -> Center, Alignment -> Left]

Above we observed that most numbers participate in a possible date interpretation. Only relatively few numbers participate in a first-occurring date interpretation — 121,470.

Length[Union[
  Flatten[Table[#2 + j, {j, 0, Length[#3] - 1}] & @@@  firstOccurrences, 1]]]

121470

Some of the position sequences overlap anyway and we can form network chains of the dates with overlapping digits sequences.

overlappingBirthdays = 
  Select[Partition[SortBy[firstOccurrencesSortedByOccurrence, #[[2]] &], 2, 1],
   #[[1, 2]] + Length[#[[1, 3]]] >= #[[2, 2]] &];

Graph[DirectedEdge[#1, #2] & @@@ Take[overlappingBirthdays, All], 
 GraphLayout -> "SpiralEmbedding"]

The next graphic shows the increasing gap sizes between consecutive dates.

NumberLinePlot[#2 & @@@ Take[firstOccurrencesSortedByOccurrence, All],
  PlotStyle -> Directive[Blue, PointSize[0.002]]]

Distribution of the gap sizes

Histogram[(#2[[2]] - #1[[2]]) & @@@ 
  Partition[firstOccurrencesSortedByOccurrence, 2, 1], {"Log", 60}, "Probability"]

Here are pairs of consecutively occurring date-interpretations that have the largest gap between them. The larger gaps were clearly visible in the penultimate graphic.

Grid[{Column[{formatOccurrence[#1], formatOccurrence[#2]}], #2[[
      2]] - #1[[2]]} & @@@ (Take[#, - 6] &@ 
    Reverse[SortBy[
      Partition[firstOccurrencesSortedByOccurrence, 2,  1], #[[1, 2]] - #[[2, 2]] &]]), Dividers -> Center, 
 Alignment -> {{Left, Right}, Center}]

Dates in other expansions and in other constants

Now, the real special dates are the ones where the concatenated continued fraction expansion position agrees with the decimal expansion position. And, interestingly, there is just one such date.

cfDigits = 
  Flatten[MapIndexed[
    Function[{d, p}, {#, p[[1]]} & /@ IntegerDigits[d]], 
    ContinuedFraction[Pi, 10000000]], 1];

Select[firstOccurrencesSortedByOccurrence, (First /@ 
       Take[cfDigits, #2 + {0, Length[#3] - 1}]) === #3 & @@ # &]

{{{2011, 1, 4}, 20494, {1, 4, 1, 1}}}

cfDigits[[20494]]

{1, 17793}

Take[ContinuedFraction[Pi, 17793 + 2], -5]

{1, 1, 14, 1, 1}

None of the calculations carried out so far were special to the digits in [Pi]. The digits of any other irrational numbers (or even sufficiently long rational numbers) contain date interpretations. Running some overnight searches, it is straightforward to find many numeric expressions that contain the dates of this year (2015). Here they are put together in some interactive demonstration.

We will now come to an end now our musings. As a last example, let us interpret digit positions as seconds since this year’s pi-time at March 14 9:26:53. How long would we have to wait until seeing the digit sequence 3\[CenterDot]1\[CenterDot]4\[CenterDot]1\[CenterDot]5 in the decimal expansion of other constants?

getPiStartPositions[\[Xi]_, seqString_] := 
  Module[{d = 5 10^StringLength[seqString], p = {}, c = 0},
   While[p === {} && c < 6, 
               
    p = StringPosition[
      StringJoin[ToString /@ RealDigits[\[Xi], 10, d][[1]]], 
      seqString , 1]; d = 2 d; c++];
               p];

\[Xi]s = {Pi, Sqrt[2], Sqrt[3], Sqrt[5], CubeRoot[3], Log[2], Log[3], 
   Cos[1], Tan[1], Sec[238], Cosh[1], Sinh[2], ArcSin[1/Pi], 
   ArcTan[1/E], EulerGamma, Pi^Pi, E^E, Pi^E, E^Pi, Zeta[2], Zeta[3], 
   ChampernowneNumber[10], ChampernowneNumber[2], 
   ChampernowneNumber[3], Gamma[1/2], Fibonacci[Pi]};

\[Pi]Posis = {#, getPiStartPositions[#, "31415"]} & /@ \[Xi]s;

Using the neat new function TimelinePlot that Brett described in his last blog post, we can easily show how long we would have to wait.

TimelinePlot[
 Association@(#1 -> 
      DatePlus[DateObject[{2015, 3, 14, 9, 26}], 
       Quantity[#2[[1, 1]], "Seconds"]] & @@@ \[Pi]Posis), PlotLayout -> "ReverseVertical"]

We let the reader explore the dates in the digits of pi more as well replace pi by another constant (for instance Euler’s constant E to justify the title of the blog), and maybe even 10 by another base. The overall, qualitative structures will be the same for almost all irrational numbers. (For a change, try ChampernowneNumber[10].) Will 10 million digits be enough to find every date in say E (where is October 21 2014)? Which special dates are hidden in other constants? This and many more things are left to explore.

[bernatep (Bernat Espigulé)]

Michael, thanks for sharing these numerological musings! Have a very happy Pi day!

[rowland (Todd Rowland)]

Great post, quite a lot of stuff here.

One of the big questions about Pi is the normality of its digit sequence (essentially how random are the digits) As far as I know, the digit sequence to all irrationals including Sqrt[2] are essentially random (see the discussion in chapter 4 of Stephen Wolfram’s A New Kind of Science )

Do Michael’s results on the dates show any sign of regularity? If we replace the pi digits by a random sequenceall of the pictures should look similar.

Is there a Random Day?

Another thing to wonder about is the distribution of birth days, which is complicated because it is not uniform for the general population, but the question is whether birthdate correlates with the likelihood of taking part in PiDay.