Im Dieser Artikel, wird eine schnelle rekursive Formulierung der Chudnovsky-Pi-Formel unter Verwendung von binärer Aufspaltung angegeben. In Python:
C = 640320
C3_OVER_24 = C**3 // 24
def bs(a, b):
if b - a == 1:
if a == 0:
Pab = Qab = 1
else:
Pab = (6*a-5)*(2*a-1)*(6*a-1)
Qab = a*a*a*C3_OVER_24
Tab = Pab * (13591409 + 545140134*a) # a(a) * p(a)
if a & 1:
Tab = -Tab
else:
m = (a + b) // 2
Pam, Qam, Tam = bs(a, m)
Pmb, Qmb, Tmb = bs(m, b)
Pab = Pam * Pmb
Qab = Qam * Qmb
Tab = Qmb * Tam + Pam * Tmb
return Pab, Qab, Tab
N = int(digits/DIGITS_PER_TERM + 1)
# Calclate P(0,N) and Q(0,N)
P, Q, T = bs(0, N)
one = 10**digits
sqrtC = sqrt(10005*one, one)
return (Q*426880*sqrtC) // T
Diese Methode ist bereits sehr schnell, aber es wird erwähnt, dass eine Implementierung auf der Website der GMP-Bibliothek, gmp-chudnovsky.c, faktorisiert auch Zähler und Nenner bei der binären Aufteilung. Da der Code optimiert und für mich schwer verständlich ist, was ist die allgemeine Idee dahinter, wie dies gemacht wird? Ich kann nicht sagen, ob Brüche vereinfacht werden, Zahlen in faktorisierter Form gehalten werden, anstatt vollständig multipliziert zu werden, oder beides.
Hier ist ein Codebeispiel der binären Aufteilung:
/* binary splitting */
void
bs(unsigned long a, unsigned long b, unsigned gflag, long int level)
{
unsigned long i, mid;
int ccc;
if (b-a==1) {
/*
g(b-1,b) = (6b-5)(2b-1)(6b-1)
p(b-1,b) = b^3 * C^3 / 24
q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
*/
mpz_set_ui(p1, b);
mpz_mul_ui(p1, p1, b);
mpz_mul_ui(p1, p1, b);
mpz_mul_ui(p1, p1, (C/24)*(C/24));
mpz_mul_ui(p1, p1, C*24);
mpz_set_ui(g1, 2*b-1);
mpz_mul_ui(g1, g1, 6*b-1);
mpz_mul_ui(g1, g1, 6*b-5);
mpz_set_ui(q1, b);
mpz_mul_ui(q1, q1, B);
mpz_add_ui(q1, q1, A);
mpz_mul (q1, q1, g1);
if (b%2)
mpz_neg(q1, q1);
i=b;
while ((i&1)==0) i>>=1;
fac_set_bp(fp1, i, 3); /* b^3 */
fac_mul_bp(fp1, 3*5*23*29, 3);
fp1[0].pow[0]--;
fac_set_bp(fg1, 2*b-1, 1); /* 2b-1 */
fac_mul_bp(fg1, 6*b-1, 1); /* 6b-1 */
fac_mul_bp(fg1, 6*b-5, 1); /* 6b-5 */
if (b>(int)(progress)) {
printf("."); fflush(stdout);
progress += percent*2;
}
} else {
/*
p(a,b) = p(a,m) * p(m,b)
g(a,b) = g(a,m) * g(m,b)
q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
*/
mid = a+(b-a)*0.5224; /* tuning parameter */
bs(a, mid, 1, level+1);
top++;
bs(mid, b, gflag, level+1);
top--;
if (level == 0)
puts ("");
ccc = level == 0;
if (ccc) CHECK_MEMUSAGE;
if (level>=4) { /* tuning parameter */
#if 0
long t = cputime();
#endif
fac_remove_gcd(p2, fp2, g1, fg1);
#if 0
gcd_time += cputime()-t;
#endif
}
if (ccc) CHECK_MEMUSAGE;
mpz_mul(p1, p1, p2);
if (ccc) CHECK_MEMUSAGE;
mpz_mul(q1, q1, p2);
if (ccc) CHECK_MEMUSAGE;
mpz_mul(q2, q2, g1);
if (ccc) CHECK_MEMUSAGE;
mpz_add(q1, q1, q2);
if (ccc) CHECK_MEMUSAGE;
fac_mul(fp1, fp2);
if (gflag) {
mpz_mul(g1, g1, g2);
fac_mul(fg1, fg2);
}
}
if (out&2) {
printf("p(%ld,%ld)=",a,b); fac_show(fp1);
if (gflag)
printf("g(%ld,%ld)=",a,b); fac_show(fg1);
}
}
Folgende Kommentare sind ausschlaggebend:
/*
g(b-1,b) = (6b-5)(2b-1)(6b-1)
p(b-1,b) = b^3 * C^3 / 24
q(b-1,b) = (-1)^b*g(b-1,b)*(A+Bb).
*/
/*
p(a,b) = p(a,m) * p(m,b)
g(a,b) = g(a,m) * g(m,b)
q(a,b) = q(a,m) * p(m,b) + q(m,b) * g(a,m)
*/
/*
p*(C/D)*sqrt(C)
pi = -----------------
(q+A*p)
*/
Beobachten:
- Wenn wir beide skalieren
p und q um einen Faktor von k in der Endrunde hat es keinen Einfluss auf das Ergebnis.
- Beim Berechnen der Zusammenführungsoperation, die dem zweiten Satz von Kommentaren entspricht, wenn wir jeden skalieren würden
p(a,m), g(a,m)und q(a,m) durch k dann würde sich dieser Faktor einfach durchsetzen p(a,b), g(a,b), q(a,b); ähnlich, wenn wir jeden skalieren würden p(m,b), g(m,b)und q(m,b) dann würde sich dieser Faktor durchsetzen.
-
Die Linie
fac_remove_gcd(p2, fp2, g1, fg1);
ist effektiv
k = gcd(p2, g1);
p2 /= k; // p(m,b) /= k
g1 /= k; // g(a,m) /= k
Dies hat den Nettoeffekt einer Herunterskalierung p(a,b), g(a,b)und q(a,b) damit gcd. Nach den beiden vorangegangenen Beobachtungen wird dieses Herunterskalieren bis zum Endergebnis sauber durchgezogen.
Nachschrift
Ich habe drei Möglichkeiten ausprobiert, die Faktorisierung in Python zu implementieren.
- Das Verfolgen der Faktorisierungen mit unveränderlichen Listen verlangsamt die Dinge schrecklich, weil es viel zu viel Arbeit gibt, die Listen zu pflegen.
- Verwenden von gmpys
gcd gibt keine Beschleunigung
- Vorberechnen einer Liste von Primzahlen bis zu
6 * N (d. h. die sich teilen könnten g) und das Testen jeder dieser Primzahlen verlangsamt die Dinge um den Faktor 2 bis 3.
Ich komme zu dem Schluss, dass eine Beschleunigung mit diesem Ansatz die Verwendung eines veränderlichen Zustands zum Verfolgen der Faktorisierungen erfordert, daher ist dies ein großer Erfolg für die Wartbarkeit.
Ich habe mir nicht den vollständigen Code angesehen, aber ich habe einen kurzen Blick darauf geworfen, um den Auszug, den Sie in Ihrer Frage angeben, besser zu verstehen.
Um einige Punkte aus Ihrer Frage zu beantworten, werfen Sie zuerst einen Blick auf dieses Stück Code:
typedef struct {
unsigned long max_facs;
unsigned long num_facs;
unsigned long *fac;
unsigned long *pow;
} fac_t[1];
Es definiert einen neuen Typ als Struktur (wenn Sie C überhaupt nicht kennen, nehmen wir an, es ist wie ein rudimentäres Pyhton-Objekt, das Variablen, aber keine Methoden einbettet). Dieser Typ ermöglicht es, ganze Zahlen als zwei ganze Zahlen und zwei Arrays (z. B. zwei Listen) zu behandeln:
- größte Faktor
- Reihe von Faktoren
- Liste der Faktoren
- Liste der (entsprechenden) Befugnisse für diese Faktoren
Gleichzeitig behält der Code die gleich Zahlen als große Ganzzahlen vom Typ libgmp (das ist gemeint mit mpz_t p und mpz_t g in den Argumenten der Funktion).
Nun, was ist mit der Funktion, die Sie zeigen. Es wird genannt fac_remove_gcd; die Initiale fac hat mit dem Namen des zuvor beschriebenen Typs zu tun; Die beiden folgenden Wörter sind leicht zu verstehen: Dividieren Sie zwei ganze Zahlen des Typs fac durch ihre gcd.
Der C-Code iteriert über die zwei Listen von Faktoren in beiden Listen; Es ist einfach, beide Listen zu synchronisieren, da die Faktoren geordnet sind (Abschnitt des Codes um die else if und else Erklärungen); Immer wenn zwei gemeinsame Faktoren erkannt werden (initial if Anweisung), ist das Dividieren eine Frage der einfachen Subtraktion: In beiden Potenzlisten wird für diesen Faktor die kleinste Potenz subtrahiert (z. B. mit a=2*5^3*17 und b=3*5^5*19, dem Wert 3 wird in den Listen der Befugnisse für subtrahiert a und b an der Position, die dem Faktor 5 entspricht, was zu a=2*5^0*17 und b=3*5^2*19 führt).
Während dieses Vorgangs wird eine Anzahl (derselben fac type) erstellt und aufgerufen fmul; es ist offensichtlich der ggT beider Zahlen.
Nach diesem Schritt rief der gcd auf fmul und von der sein fac type wird mit der aufgerufenen Funktion (ebenfalls im Code des Programms) in eine GMP Big Integer umgewandelt bs_mul. Dies ermöglicht es, den ggT als große ganze Zahl zu berechnen, um die neuen Werte der dividierten Zahlen in beiden Formen zu synchronisieren: große ganze Zahlen und Sonderformen fac Typ. Sobald der ggT als große ganze Zahl berechnet ist, ist es einfach, beide Anfangszahlen durch den ggT zu dividieren.
Somit wirken die Funktionen auf zwei Versionen jeder Anfangszahl.
Hoffe es kann helfen.
