Angelic and Evil Numbers

Copyright Chris Lomont, 2004

This notebook explains how to compute the "probability that an irrational number has initial digits summing to 666".

This was motivated by the fact mentioned on the website www . mathpuzzle . com on Sept 29,2004, that adding the initial (after the decimal) 144 digits of π  gives a sum of 666, and it was also noted this works for summing the first 146 digits after the decimal point of the Golden Ratio ϕ.  Here we note the built-in Mathematica Glaisher constant also has a sum of 666 for the first 153 digits after the decimal point. The question was posed: how likely is this? To analyze this problem, first we make a function to find the sum of m digits, taking every nth digit after the decimal point:

digitSum[r_, n_, m_]  :=  
Apply[
        Plus, 
Take[
RealDigits[
FractionalPart[
N[r, 
                  100 + n m  +  
                    Length[IntegerDigits[
                        IntegerPart[r]] ] ]]][[1]],   (* 
            note we take many more digits *) 
{1, 
            n m,  n}]
]

digitSum[π, 1, 144]

digitSum[GoldenRatio, 1, 
      146]

digitSum[Glaisher, 1, 
      153]

We call such numbers Evil. Note that this does not work for the natural base e, so not all numbers are Evil:

{digitSum[E, 1, 141], digitSum[E, 1, 142]}

However, if we add every 15th digit of e, we find by adding 137 digits, that e is somewhat evil....

digitSum[E, 15, 137]

Like a good numerologist, we find similarly the smallest skip making each common Mathematica constant evil, ranking them by total digits needed to detect their evilness. Note it took a lot of digits to find 666 in the Catalan number.

digitSum[Degree, 4, 153]

digitSum[EulerGamma, 6, 139]

digitSum[(Sqrt[2]), 6, 152]

digitSum[Khinchin, 7, 155]

digitSum[E, 15, 137]

digitSum[Catalan, 28, 144]

Mysterious! Are all numbers inherently evil? All this really shows is not to trust numerology books like the Bible Code. It's easy to find number patterns if you look and fiddle with data.

To show that it is easy to find these, we will show the probability of a sequence of random digits adding to an arbitrary positive integer k is roughly 1/5,and give a way to compute the exact probability p[k]. Thus we are quite likely to find evil numbers (for some digit skip) almost everywhere. In fact, I conjecture they are dense in the set of real numbers.

First some reductions, all intuitively stated, although they can be rigorously proven. Consider a random variable X that returns a digit from 0 to 9 inclusive with uniform probability. Then consider all sequences ((Subscript[X, 1]), (Subscript[X, 2]),  ... , (Subscript[X, n]),  ... ). Note this returns an irrational number with probability 1 (rational numbers with period m must have repeating sequences after the first m digits, which has 0 chance of happening). We want the probability p[k], given a positive integer k, that there is an integer m such that the sequence ((Subscript[X, 1]), (Subscript[X, 2,]) ... , (Subscript[X, m]))sums to k. For an example let's compute p[2] the hard way. We can get 2 for the sum from two types of sequences: all zeros ended by a 2, and a sequence of all zeros and a single 1 in any order, followed by a 1. Computing these probabilities takes a bit of work, since we need two infinite sums.

p[2] = (1/10)Sum[((((1/10)))^m), {m, 0, Infinity}]+Sum[(m - 1) ((((1/10)))^m), {m, 2, 
            Infinity}]

Note these sums can be done by hand. The proof to do the first one with a geometric series can easily be adapted and used twice to compute the second sum. But there is a better way.

Note the 0's in the sequence do not affect the sum, so we want to remove 0s from the problem. First, note the probability is 0 of the variable X returning all 0s  (note that it could happen, it just has arbitrarily small probability, thus 0). Secondly, note that for a fixed positive integer n, the probability is 0 that the sequence sum never surpasses n, since that would mean that after at most n digits X returns all 0s, again which has 0 probability. Thus we are safe in removing the 0s, since they will not affect the final probabilities p[k]. So redefine X such that X only returns digits 1 to 9, with uniform probability.

To show how useful this is, let's recompute p[2] the new way. There are now only two ways to do it-  a sequence of {2} and a sequence {1,1} with probability:

p[2]  =  (1/9) + (1/9) (1/9)

Note we get the same answer as before (which is nice). To compute the general formula, we make some observations. For k<10, the probability of getting a digit sum of k is the probability of getting k on the first digit, plus the probabilities of getting a smaller value times the probability of jumping to k from the smaller value. This defines the probabilities recursively. We note p[1] = (1/9). Thus we have for k<10

Clear[p]

p[1]  =  (1/9);

p[k_ /; k < 10]  :=  p[k]  =  
          (1/9) + (1/9) Sum[p[j], {j, 1, k - 1}];

Table[p[k], {k, 1, 9}]

We note that this is a geometric series (which unfortunately is not repeated for k>9). There should be a nice way to see these probabilities directly, but I have not found it.

Table[((10^(n - 1))/(9^n)), {n, 1, 9}]

Next we use a similar observation to compute p[k] for k>9. In this case, to get a sum of k means before the last digit d is added obtaining k, we were at k-d. So for each last possible digit 1-9, we compute the probability of reaching k-d by any means, times the 1/9 probability of getting last digit d. This gives

p[
          k_ /; k > 9]  :=  p[k]  =  (1/9) 
            Sum[p[k - d], {d, 1, 9}];

We used the Mathematica trick of reassigning each value with p[k_] := p[k]=... which speeds up computation greatly by storing precomputed values, at the cost of storage space. Unfortunately, if we try to evaluate p[666], we get an error since we exceed the recursion limit. So we "precompute" values, say the first 1000, using a trick to do 100 at a time. Plus this is much faster than trying to compute p[666] directly.

Timing[Table[p[100 n], {n, 1, 10}]][[1]]

We just created the first 1000 values in memory. Now we get the exact answer! What is the probability of a uniformly chosen infinite string of digits containing a prefix sum of 666? It is.....

p[666]

Hmmm, not very nice. However we show how close it is to 20%. We compute to 75 digits, otherwise we only get noise since the machine precision is usually around (10^(-17)).

N[p[666]]

1/5  -  N[p[666], 75]

So the probability of a given number being Evil right off the bat is 20%. Thus is it very likely that by fiddling with the spacing between digits summed that any irrational number will turn out to be evil. What is the probability of summing to 7 and being Angelic?

N[p[7]]

At least this gives some hope - a random number is more likely to be Angelic (20.9%) than Evil (20%). That was close! All the constants tested above for (somewhat) Evilness become (somewhat) Angelic with less digits than it takes to find Evil, except the Catalan constant, since the first digit after the decimal is 9, and so it has no Angelic possibilities. Also the Angelic numbers cannot be dense, since some irrationals begin 0.8... or 0.9... So we find out (again) that we can find patterns wherever we look. Here are the constants again, showing Angelic properties with minimal digits.

digitSum[E, 1, 1]

digitSum[GoldenRatio, 1, 2]

digitSum[EulerGamma, 3, 2]

digitSum[Glaisher, 2, 4]

digitSum[π, 6, 2]

digitSum[(Sqrt[2]), 4, 3]

digitSum[Khinchin, 8, 2]

digitSum[Degree, 18, 2]

How close is p[k] to 20% for the largest k we computed? Note we have to fiddle with the precision to get good answers. We find the maximum distance the p[k] varies from 20% for k from 900 to 1000.

Max[Abs[Table[N[p[k], 105] - 1/5, {k, 900, 1000}]]]

Very close.... It seems ALL the probabilities are converging to 20%. Proving this convergence would be a good challenge.

Some pictures to show the self similar nature of the graph.

To show some of how this probability behaves, let's graph it for various intervals of k.

sh[k1_, k2_] := ListPlot[
        Table[{k, p[k]}, {k, k1, k2}], PlotJoined -> True]

sh[1, 100]

So it looks wild at first, but perhaps is smooths out? Let's look closer.

sh[1, 20]

Is it  jumpy early on. Does it stabilize? Let's look at some intervals.

sh[20, 40]

sh[40, 80]

sh[80, 120]

sh[120, 160]

sh[160, 200]

sh[200, 600]

The first several graphs are similar at different scales - fractals?  The graph tapers off to near a probability of 1/5, but still oscillates. It is likely we have run up against the numerical limits of Mathematica, and have not actually stabilized. To test this hypothesis, we compute a intervals as above, and measure the difference from 1/5 to see that the fractal oscillations are still occurring:

ListPlot[Table[p[k] - 1/5, {k, 550, 600}], 
      PlotJoined -> True]

ListPlot[Table[p[k] - 1/5, {k, 950, 1000}], 
      PlotJoined -> True]

So it seems this probability function is oscillates and converges to 20% as we thought, even at large k values.

A pressing question: can we find a closed form? Mathematica has a built in function to solve for this, but it may be a very hard recurrence to solve since each p[k] relies on 9 previous terms for k>9. Using the RSolve command ran for several hours, but had no results by then....

eqs  =  
      Join[Table[
          q[i] == p[i], {i, 1, 
            9}], {q[k]  ==  (1/9) 
              Sum[q[k - d], {d, 1, 9}]}]

We can show two nice forms for the generating function (found by hand) of p[k]. Finding a closed form from this is a lot of work (I tried).

First note we can replace the recursion by p[k+1] = (10 p[k]-p[k-9])/9 for k ≥ 10 by computing p[k+1]-p[k], and simplifying. This eases the work of finding the generating function. Even this simpler form does not work well in RSolve.

f[t_]  :=  ((1 - (t^9))/((t^10) - 10  t + 9))

CoefficientList[Series[f[t], {t, 0, 999}], 
        t] == Table[p[k], {k, 1, 1000}]

f2[
        t_] := ((((t^6) + (t^3) + 1) ((t^2) + t + 
                1))/((1 - t) (1 (t^8) + 2 (t^7) + 3 (t^6) + 
                4 (t^5) + 5 (t^4) + 6  (t^3) + 7  (t^2) + 8  (t^1) + 
                9)))

Reduce[f[t]  ==  f2[t], t]

So the probability p[k] is now found as the kth coefficient in the Taylor Series of f[t] .

Here are some misc computations used while trying to find a closed form, which show a nice closed form for k < 20. These formulae are able to be found for higher k, but they get very messy.

q[k_ /;  k < 10]  :=  ((10^(k - 1))/(9^k))

q[k_ /;  
          10  ≤  k  <  
            20]  :=  (((10^(k - 1)) - k (10^(k - 11)) (9^9))/(9^k))

Union[Table[q[k] == p[k], {k, 1, 19}]]

And a final sanity check - test some random digit sequences against the probability formula.

digSeq[n_] := Table[
        Random[Integer, {1, 9}], {n}]

tests  =  Table[digSeq[200], {10000}] ;

totals  =  Table[0, {200*9}] ;

countSeq[
        s_] := Module[{curSum = 0, 
          pos}, 
For[pos  =  1,  
          pos  <=  Length[s],  pos++, 
          curSum  +=  
            s[[pos]]; totals[[curSum]]++]; 
Return[{}];
]  (* 
        count the sums found in a given sequence,  
        and add the counts to totals *)

Timing[countSeq[#] &  /@  tests][[1]]

probs  =  totals/Length[tests] ;

err  =  
        Take[probs, {1, 200}] - 
          Table[p[i], {i, 1, 200}] ;

N[Max[Abs[err]]]

th[a_, b_]  :=  
      ListPlot[Take[probs, {a, b}], PlotJoined -> True]

th[1, 20]

sh[1, 20]

th[1, 100]

sh[1, 100]

Although the experimental data has a lot of noise, it is comforting that it also seems to center around a probability of 20% for large k.

THE END. I hope you have enjoyed this. Visit www . math . purdue . edu/~clomont or www . lomont . org.