Python 2, 110 bajtów
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
Maksymalna liczba cyfr do sprawdzenia jest pobierana od standardowego wejścia. 10 000 cyfr kończy się w ciągu około 2 sekund przy PyPy 5.3.
Przykładowe użycie
$ echo 10000 | pypy pi-runs.py
3
33
111
9999
99999
999999
Coś użytecznego
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
W tym celu przeszedłem z Chudnovsky'ego na Ramanujan 39. Chudnovsky'emu zabrakło pamięci w moim systemie krótko po 100 milionach cyfr, ale Ramanujan dotarł do 400 milionów w zaledwie około 38 minut. Myślę, że to kolejny przypadek, w którym wygrywa wolniejszy wzrost terminów, przynajmniej w systemie z ograniczonymi zasobami.
Przykładowe użycie
$ 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
Szybsze generatory bez ograniczeń
Implementacja referencji podana w opisie problemu jest interesująca. Wykorzystuje nieograniczony generator, pobrany bezpośrednio z papierowego nieograniczonego algorytmu czopów dla cyfr Pi . Według autora dostarczone implementacje są „celowo niejasne”, więc postanowiłem wprowadzić nowe implementacje wszystkich trzech algorytmów wymienionych przez autora, bez celowego zaciemniania. Dodałem także czwarty, oparty na 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)
Notatki
Powyżej jest 6 implementacji: dwie implementacje referencyjne dostarczone przez autora (oznaczone _ref) i cztery, które obliczają partie, generując wiele cyfr jednocześnie ( _md). Wszystkie wdrożenia zostały potwierdzone do 100 000 cyfr. Wybierając wielkości partii, wybrałem wartości, które z czasem powoli tracą precyzję. Na przykład g1_mdgeneruje 10 cyfr na partię, z 33 iteracjami. Jednak da to tylko ~ 9,93 poprawnych cyfr. Gdy skończy się precyzja, warunek sprawdzania zakończy się niepowodzeniem, co spowoduje uruchomienie dodatkowej partii do uruchomienia. Wydaje się to być bardziej wydajne niż powoli zyskiwać dodatkową, niepotrzebną precyzję w czasie.
- g1 (Leibniz / Euler) Przechowywana jest
dodatkowa zmienna jreprezentująca 2*i+1. Autor robi to samo w implementacji referencyjnej. Obliczanie nosobno jest znacznie prostsze (i mniej niejasne), ponieważ wykorzystuje bieżące wartości q, ra tzamiast następnego.
- g2 (Lambert)
Czek n == q/sjest wprawdzie dość luźny. To powinno przeczytać n == (q*(k+2*j+4)+r)/(s*(k+2*j+4)+t), gdzie jjest 2*i-1i kjest i*i. Przy wyższych iteracjach ri tterminy stają się coraz mniej znaczące. Jak na razie jest dobry na pierwsze 100 000 cyfr, więc prawdopodobnie jest dobry dla wszystkich. Autor nie podaje żadnej referencyjnej implementacji.
- g3 (Gosper)
Autor przypuszcza, że nie trzeba sprawdzać, nczy nie zmieni się w kolejnych iteracjach, a jedynie spowalnia algorytm. Chociaż prawdopodobnie jest to prawda, generator przechowuje o około 13% więcej poprawnych cyfr niż obecnie generuje, co wydaje się nieco marnotrawstwem. Ponownie dodałem czek i czekam, aż 50 cyfr będzie poprawnych, generując je wszystkie naraz, z zauważalnym wzrostem wydajności.
- g4 (Ramanujan 39)
Obliczony jako
Niestety, snie zeruje się ze względu na początkowy skład (3528 ÷), ale wciąż jest znacznie szybszy niż g3. Konwergencja wynosi ~ 5,89 cyfr na termin, naraz generowanych jest 3511 cyfr. Jeśli to trochę za dużo, generowanie 271 cyfr na 46 iteracji jest również dobrym wyborem.
Czasy
Podjęte w moim systemie, wyłącznie w celach porównawczych. Czasy podano w sekundach. Jeśli pomiar czasu trwał dłużej niż 10 minut, nie przeprowadzałem żadnych dalszych testów.
| 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
Interesujące jest to, że g2ostatecznie wyprzedza g3, pomimo wolniejszej konwergencji. Podejrzewam, że dzieje się tak, ponieważ operandy rosną znacznie wolniej, wygrywając na dłuższą metę. Najszybsza implementacja g4_mdjest około 235 razy szybsza niż g3_refimplantacja na 500 000 cyfr. To powiedziawszy, nadal istnieje znaczne obciążenie związane z przesyłaniem strumieniowym cyfr w ten sposób. Obliczanie wszystkich cyfr bezpośrednio przy użyciu Ramanujan 39 ( źródło python ) jest około 10 razy szybsze.
Dlaczego nie Chudnovsky?
Algorytm Chudnovsky'ego wymaga pierwiastka kwadratowego o pełnej precyzji, w którym szczerze mówiąc nie jestem pewien, jak pracować - zakładając, że w ogóle może być. Ramanujan 39 jest pod tym względem wyjątkowy. Jednak metoda wydaje się sprzyjać formułom podobnym do Machina, takim jak te stosowane przez y-crunchera, więc może to być droga warta poznania.