Python 2, 110 bytes
n=input()
x=p=7*n|1
while~-p:x=p/2*x/p+2*10**n;p-=2
l=m=0
for c in`x`:
l=l*(p==c)+1;p=c
if l>m:m=l;print p*l
The maximum number of digits to check is taken as from stdin. 10,000 digits finishes in about 2s with PyPy 5.3.
Sample Usage
$ echo 10000 | pypy pi-runs.py
3
33
111
9999
99999
999999
Something Useful
from sys import argv
from gmpy2 import mpz
def pibs(a, b):
if a == b:
if a == 0:
return (1, 1, 1123)
p = a*(a*(32*a-48)+22)-3
q = a*a*a*24893568
t = 21460*a+1123
return (p, -q, p*t)
m = (a+b) >> 1
p1, q1, t1 = pibs(a, m)
p2, q2, t2 = pibs(m+1, b)
return (p1*p2, q1*q2, q2*t1 + p1*t2)
if __name__ == '__main__':
from sys import argv
digits = int(argv[1])
pi_terms = mpz(digits*0.16975227728583067)
p, q, t = pibs(0, pi_terms)
z = mpz(10)**digits
pi = 3528*q*z/t
l=m=0
x=0
for c in str(pi):
l=l*(p==c)+1;p=c
if l>m:m=l;print x,p*l
x+=1
I've switched from Chudnovsky to Ramanujan 39 for this. Chudnovsky ran out of memory on my system shortly after 100 million digits, but Ramanujan made it all the way to 400 million, in only about 38 minutes. I think this is another case were slower growth rate of terms wins out in the end, at least on a system with limited resources.
Sample Usage
$ python pi-ramanujan39-runs.py 400000000
0 3
25 33
155 111
765 9999
766 99999
767 999999
710106 3333333
22931752 44444444
24658609 777777777
386980421 6666666666
Faster Unbounded Generators
The reference implementation given in the problem description is interesting. It uses an unbounded generator, taken directly from the paper Unbounded Spigot Algorithms for the Digits of Pi. According to the author, the implementations provided are "deliberately obscure", so I decided to make fresh implementations of all three algorithms listed by the author, without deliberate obfuscation. I've also added a fourth, based on Ramanujan #39.
try:
from gmpy2 import mpz
except:
mpz = long
def g1_ref():
# Leibniz/Euler, reference
q, r, t = mpz(1), mpz(0), mpz(1)
i, j = 1, 3
while True:
n = (q+r)/t
if n*t > 4*q+r-t:
yield n
q, r = 10*q, 10*(r-n*t)
q, r, t = q*i, (2*q+r)*j, t*j
i += 1; j += 2
def g1_md():
# Leibniz/Euler, multi-digit
q, r, t = mpz(1), mpz(0), mpz(1)
i, j = 1, 3
z = mpz(10)**10
while True:
n = (q+r)/t
if n*t > 4*q+r-t:
for d in digits(n, i>34 and 10 or 1): yield d
q, r = z*q, z*(r-n*t)
u, v, x = 1, 0, 1
for k in range(33):
u, v, x = u*i, (2*u+v)*j, x*j
i += 1; j += 2
q, r, t = q*u, q*v+r*x, t*x
def g2_md():
# Lambert, multi-digit
q, r, s, t = mpz(0), mpz(4), mpz(1), mpz(0)
i, j, k = 1, 1, 1
z = mpz(10)**49
while True:
n = (q+r)/(s+t)
if n == q/s:
for d in digits(n, i>65 and 49 or 1): yield d
q, r = z*(q-n*s), z*(r-n*t)
u, v, w, x = 1, 0, 0, 1
for l in range(64):
u, v, w, x = u*j+v, u*k, w*j+x, w*k
i += 1; j += 2; k += j
q, r, s, t = q*u+r*w, q*v+r*x, s*u+t*w, s*v+t*x
def g3_ref():
# Gosper, reference
q, r, t = mpz(1), mpz(180), mpz(60)
i = 2
while True:
u, y = i*(i*27+27)+6, (q+r)/t
yield y
q, r, t, i = 10*q*i*(2*i-1), 10*u*(q*(5*i-2)+r-y*t), t*u, i+1
def g3_md():
# Gosper, multi-digit
q, r, t = mpz(1), mpz(0), mpz(1)
i, j = 1, 60
z = mpz(10)**50
while True:
n = (q+r)/t
if n*t > 6*i*q+r-t:
for d in digits(n, i>38 and 50 or 1): yield d
q, r = z*q, z*(r-n*t)
u, v, x = 1, 0, 1
for k in range(37):
u, v, x = u*i*(2*i-1), j*(u*(5*i-2)+v), x*j
i += 1; j += 54*i
q, r, t = q*u, q*v+r*x, t*x
def g4_md():
# Ramanujan 39, multi-digit
q, r, s ,t = mpz(0), mpz(3528), mpz(1), mpz(0)
i = 1
z = mpz(10)**3511
while True:
n = (q+r)/(s+t)
if n == (22583*i*q+r)/(22583*i*s+t):
for d in digits(n, i>597 and 3511 or 1): yield d
q, r = z*(q-n*s), z*(r-n*t)
u, v, x = mpz(1), mpz(0), mpz(1)
for k in range(596):
c, d, f = i*(i*(i*32-48)+22)-3, 21460*i-20337, -i*i*i*24893568
u, v, x = u*c, (u*d+v)*f, x*f
i += 1
q, r, s, t = q*u, q*v+r*x, s*u, s*v+t*x
def digits(x, n):
o = []
for k in range(n):
x, r = divmod(x, 10)
o.append(r)
return reversed(o)
Notes
Above are 6 implementations: the two reference implementations provided by the author (denoted _ref
), and four that compute terms in batches, generating multiple digits at once (_md
). All implementations have been confirmed to 100,000 digits. When choosing batch sizes, I chose values that slowly lose precision over time. For example, g1_md
generates 10 digits per batch, with 33 iterations. However, this will only produce ~9.93 correct digits. When precision runs out the check condition will fail, triggering an extra batch to be run. This seems to be more performant than slowly gainly extra, unneeded precision over time.
- g1 (Leibniz/Euler)
An extra variable j
is kept, representing 2*i+1
. The author does the same in the reference implementation. Calculating n
separately is far simpler (and less obscure), because it uses the current values of q
, r
and t
, rather than the next.
- g2 (Lambert)
The check n == q/s
is admittedly quite lax. That should read n == (q*(k+2*j+4)+r)/(s*(k+2*j+4)+t)
, where j
is 2*i-1
and k
is i*i
. At higher iterations, the r
and t
terms become increasingly less significant. As is, it's good for the first 100,000 digits, so it's probably good for all. The author provides no reference implementation.
- g3 (Gosper)
The author conjectures that it is unnecessary to check that n
will not change in subsequent iterations, and that it only serves to slow the algorithm. While probably true, the generator is holding onto ~13% more correct digits than it has currently generated, which seems somewhat wasteful. I've added the check back in, and wait until 50 digits are correct, generating them all at once, with a noticable gain in performance.
- g4 (Ramanujan 39)
Calculated as
![](../../I/static/images/306493845cb2e55da6f85e1f81dbe93c0609d4b35e597d32502d6f50a6d75d0f)
Unfortunately, s
doesn't zero out, due to the initial (3528 ÷) composition, but it's still significantly faster than g3. Convergence is ~5.89 digits per term, 3511 digits are generated at a time. If that's a bit much, generating 271 digits per 46 iterations is also a decent choice.
Timings
Taken on my system, for purposes of comparison only. Times are listed in seconds. If a timing took longer than 10 minutes, I didn't run any further tests.
| g1_ref | g1_md | g2_md | g3_ref | g3_md | g4_md
------------+---------+---------+---------+---------+---------+--------
10,000 | 1.645 | 0.229 | 0.093 | 0.312 | 0.062 | 0.062
20,000 | 6.859 | 0.937 | 0.234 | 1.140 | 0.250 | 0.109
50,000 | 55.62 | 5.546 | 1.437 | 9.703 | 1.468 | 0.234
100,000 | 247.9 | 24.42 | 5.812 | 39.32 | 5.765 | 0.593
200,000 | 2,158 | 158.7 | 25.73 | 174.5 | 33.62 | 2.156
500,000 | - | 1,270 | 215.5 | 3,173 | 874.8 | 13.51
1,000,000 | - | - | 1,019 | - | - | 58.02
It's interesting that g2
eventually overtakes g3
, despite a slower rate of convergence. I suspect this is because the operands grow at a significantly slower rate, winning out in the long run. The fastest implmentation g4_md
is approximately 235x faster than the g3_ref
implmentation on 500,000 digits. That said, there's still a significant overhead to streaming digits in this way. Calculating all digits directly using Ramanujan 39 (python source) is approximately 10x as fast.
Why not Chudnovsky?
The Chudnovsky algorithm requires a full precision square root, which I'm honestly not sure how to work in - assuming it could be at all. Ramanujan 39 is somewhat special in this regard. However, the method does seem like it might be conducive to Machin-like formulas, such as those used by y-cruncher, so that might be an avenue worth exploring.
1Related - that runs of 9s caused headaches for many answers :P – Mego – 2016-07-13T21:26:46.490
Are you allowed to start output with the empty sequence? – LegionMammal978 – 2016-07-14T00:03:59.830
2Also, the next term of the sequence appears to be 3333333 in digits 10^-710100 through 10^-710106. The value for n = 8 does not appear in the first 5 000 000 digits. – LegionMammal978 – 2016-07-14T00:08:45.880
Do we have to do it for base 10, or can we choose our base (integer>1)? (I'm thinking this may help....)
– msh210 – 2016-07-14T05:11:33.930@msh210 I want to do that exact thing, except I can't figure out how to implement the formula to find the nth digit of
π
... – R. Kap – 2016-07-14T05:22:32.510@msh210 It must be base 10. – mbomb007 – 2016-07-14T13:28:24.650
@R.Kap Just Google "nth digit of pi". Here's the Spigot algorithm that I used.
– mbomb007 – 2016-07-14T13:29:26.0104Two more terms: 44444444 at digits 10^-22931745 through 10^-22931752 and 777777777 at digits 10^-24658601 through 10^-24658609. The value for n = 10 does not appear in the first 100 000 000 digits. – LegionMammal978 – 2016-07-14T18:02:06.347
1One more term: 6666666666 at 10^-386980412. The 11th term does not appear in the first 2 000 000 000 digits. – primo – 2016-07-19T04:25:24.903
@primo How long did that take to find? – mbomb007 – 2016-07-19T13:50:23.493
@mbomb007 400m took only 38 minutes, surprisingly. I've updated the code in my answer. – primo – 2016-07-19T13:53:25.267