Re: Re: A problem in Pi digits as Lattice space filling
- To: mathgroup at smc.vnet.net
- Subject: [mg93976] Re: [mg93843] Re: A problem in Pi digits as Lattice space filling
- From: danl at wolfram.com
- Date: Sun, 30 Nov 2008 07:00:02 -0500 (EST)
- References: <ggj7co$j2i$1@smc.vnet.net>
> A friend and long time Mathematica programmer came up with this:
> Clear[f, n]
> f[0] = 1;
> f[n_Integer /; n > 0] := Module[{cnt, m, p},
> cnt = 0; m = Table[0, Evaluate[Apply[Sequence, Table[{10}, {n}]]]];
> p = Partition[RealDigits[Pi, 10, 20*10^n][[1]], n, 1];
> For[i = 1, i <= 20*10^n - n + 1, i++,
> If[m[[Apply[Sequence, p[[i]] + 1]]]++ == 0,
> If[++cnt == 10^n, Print[i]; Break[]]];
> ]
> ];
> Table[f[n], {n, 0, 4}]
> From In[9]:=
> 33
>
> From In[9]:=
> 606
>
> From In[9]:=
> 8554
>
> From In[9]:=
> 99847
> There my machine gives up in real time...
> I pretty much trust this and my older clunky programs.
> When run on a faster machine the answers are:
> starting at f[0]=1;
> 1,33,606,8554,99847,1369561, 14118308
> the hyperbolic probability limit:
>
> Limit[10^n/f[n],n->Infinity]
>
> appears to be zero from these numbers.
> That result would suggest that the Pi digits are really not normal at all.
I do not see that implication.
> But with just these few numbers it is hard to say.
> Since the BBP approach is actually based on the normal character of the
> Pi digits
The existance of the BBP formula probably has no bearing on normality of
Pi. That at least is what I heard Helaman Fergusen say about it, 11 years
ago.
> ( one of the assumptions necessary for a PSLQ digit solve to actually
> work! I'm not certain
> of the theory, but that is what the papers say.).
If there are integer relations to be found amongst the digits, PSLQ can
find them. I am not aware of a connection here to normality, or lack
thereof.
> I'm not willing to argue the theory of the point, as I really don't much
> understand it.
> The idea of the space filling character
> of the Pi digits came from them being normal.
> If the lattice diverges as it seems to toward Infinity,
> that would imply
> that the digits of Pi aren't really normal,
> I think.
No. That is to say, if the lattice behaves as this one is behaving, there
is nothing to indicate non-normality. I'll explain further below.
> As a result the BBP Pi digits base 16 forms are not derived
> correctly?
> Or at least there is a question,
> now, involved.
The derivation was proven, many years ago. It's not conjectural.
> Me, I'm an experimenter and not an axiomatic theoretical type.
> But in any case the question of the space filling limit seems important
> to my mind.
>
>
> Respectfully, Roger L. Bagula
> 11759Waterhill Road, Lakeside,Ca 92040-2905,tel: 619-5610814
> :http://www.geocities.com/rlbagulatftn/Index.html
> alternative email: rlbagula at sbcglobal.net
There are ways to investigate this. A relatively simple approach is to use
random large numbers, and check the behavior of the same lattice filling
as you are doing with digits of Pi. I will use my version of the filling
time code, from an earlier post (your friend's is similar but slightly
slower due to more work in the lookup phase).
latticeFillMark[n_?NumericQ, d_Integer] /; d > 0 :=
Catch[Module[{digits, vals, len = 10^d, hitlist, count = 0, elem},
digits = First[RealDigits[n, 10, 10*Ceiling[d/2]*len]];
vals = Map[FromDigits, Partition[digits, d, 1]];
hitlist = ConstantArray[0, len];
Do[elem = vals[[j]] + 1; If[hitlist[[elem]] == 1, Continue[]];
hitlist[[elem]] = 1; count++;
If[count == len, Throw[j]];, {j, Length[vals]}]; $Failed]]
Here we generate random reals, of sufficient length to enable the
experimental analysis.
randomFillMarks[digits_, len_] := Table[
latticeFillMark[
RandomReal[WorkingPrecision -> 10*Ceiling[digits/2]*10^digits],
digits]
, {len}]
We create a sample of 100 cases wherein we check how long to fill the
10x10 lattice of 100 values.
r2sample = N[randomFillMarks[2, 100]];
Now I show that the mean and standard deviation make a value of 606 quite
plausible.
Mean[r2sample]
Out[77]= 517.09
In[86]:= StandardDeviation[r2sample]
Out[86]= 130.531
Here I perform a similar experiment on 20 random values, this time going
after the 10x10x10 lattice of 1000 values.
r3sample = N[randomFillMarks[3, 20]]
Out[80]= {7145., 6530., 7058., 5903., 6476., 7655., 7303., 6636., \
11536., 7789., 6639., 8176., 6445., 6553., 6606., 5620., 6604., \
10084., 6408., 10454.}
In[82]:= Mean[r3sample]
Out[82]= 7381.
In[84]:= StandardDeviation[r3sample]
Out[84]= 1567.77
Again, the value for Pi (8554) is within one (estimated) standard
deviation of the experimental mean. And we see similarly for the
10x10x10x10 case as well.
r4sample = N[randomFillMarks[4, 20]];
In[182]:= Mean[r4sample]
Out[182]= 92447.8
In[183]:= StandardDeviation[r4sample]
Out[183]= 10790.9
There is an important caveat here. It is imperative that the random number
generator be extremely good, that is, not suffer from defects of
correlation between generated digits. The default RNG is based on a
cellular automaton and has, I believe, performed well on the l'Ecuyer test
suite. So I trust it. But if you do not, you could first do
SeedRandom[Method -> "Rule30CA"]
This uses the Rule 30 cellular automaton RNG. It's slower than the default
one, but has always seemed quite robust in terms of behavior.
Finally I will indicate an exact approach to these computations of
expected value and standard deviation. We assume the digits of our input
(Pi, in this case) are uniformly independently distributed (uid), and try
to analyze the expected value of our lattice filling thresholds. The code
below does this. In brief, it works as follows. We let p[j,n,d] denote the
quantity of cases of n uid values, each in range 1-d, filling exactly j
positions. This is seen to be the number filling at most that many, minus
the number that exactly fit 1, 2, ..., j-1 positions. Then we let q[j,n,d]
denote the probability of this happening (the number of cases, divided by
total number of such ensembles, which is d^n). Why p for quantity and q
for probability? I think my typing fingers got away from me (both of
them).
Clear[p, q]
p[0, n_, d_] = 0;
q[0, n_, d_] = 0;
p[j_, n_, d_] :=
p[j, n, d] =
Together[Binomial[d,
j]*(j^n - Sum[Binomial[j, k] p[k, n, k], {k, j - 1}])]
q[j_, n_, d_] := q[j, n, d] = p[j, n, d]/d^n
I have not been able to come up with any closed form. But this suffices
for computing at least the 10 and 10x10 cases exactly. The values of
interest are q[d,n,d] with d=10 or 100, and n>=d.
In[13]:= q10[n_] = q[10, n, 10]
Out[13]= 10^-n (-10 + 45 2^n + 45 2^(3 n) + 105 2^(1 + 2 n) -
10 3^(2 n) - 40 3^(1 + n) - 252 5^n + 35 6^(1 + n) - 120 7^n + 10^
n)
In[29]:= d10[n_] = Together[q10[n] - q10[n - 1]]
Out[29]= -(1/63) 2^(-2 -
n) 5^-n (-22680 + 2835 2^(3 n) + 2835 2^(4 + n) +
19845 2^(2 + 2 n) - 280 3^(2 n) - 7840 3^(2 + n) +
245 2^(4 + n) 3^(2 + n) - 63504 5^n - 12960 7^n)
q100[n_] = q[100, n, 100];
d100[n_] = Together[q100[n] - q100[n - 1]];
We can now compute expected values.
In[192]:= mean10 = Sum[n*d10[n], {n, 10, Infinity}]
Out[192]= 7381/252
In[193]:= N[mean10]
Out[193]= 29.2897
Recall that this is where Pi gave 33. Now we'll compute the standard
deviation.
In[194]:=
var10 = Sum[n^2*d10[n], {n, 10, Infinity}] - mean10^2;
In[195]:= std10 = Sqrt[var10];
In[196]:= N[std10]
Out[196]= 11.211
So 33 is well within a standard deviation of the mean. We'll repeat this
exact computation for the 10x10 case.
mean100 = Sum[n*d100[n], {n, 100, Infinity}];
In[204]:= N[mean100]
Out[204]= 518.738
In[205]:= var100 = Sum[n^2*d100[n], {n, 100, Infinity}] - mean100^2;
In[206]:= std100 = Sqrt[var100];
N[std100]
Out[207]= 125.822
Again we see that the computed value for Pi, 8554, is within a standard
deviation of the expected value.
This method runs out of steam around here due to computational complexity.
I have not as yet found a more tractable formulation. In particular I have
not found an exact form for the general case (that is, one for all d,n
pairs). I would guess this has been studied in the literature, and that if
one exists it is well known, just not to me.
Daniel Lichtblau
Wolfram Research