Dodekaeder    MB Matheblog # 5 Inhalt Blog
voriger Eintrag
zur Leitseite
Index der gesamten Website


2020-11-14


Partitionen natürlicher Zahlen
Formel von Hardy-Ramanujan-Rademacher
                    Kommentare sind willkommen.


Srinivasa Ramanujan und Godfrey Harold Hardy publizierten 1918 eine Formel für die Partitionsfunktion \(~P(n)\), die die Anzahl der Partitionen einer natürlichen Zahl \(~n~\) angibt [1]. Am Anfang dieser Arbeit stand die bewundernswerte zahlentheoretische Intuition Ramanujans; Hans Rademacher veröffentlichte 1937 eine Weiterentwicklung der Formel in eine besser handhabbare Gestalt [3]; beide Varianten findet man unter den Stichworten Formel von Hardy-Ramanujan-Rademacher und Approximation der Partitionsfunktion.

Fredrik Johansson gab in [2] eine Wertung der Formel: "One of the most astonishing number-theoretical discoveries of the 20th century is the Hardy-Ramanujan-Rademacher formula first given as an asymptotic expansion by Hardy and Ramanujan in 1917 and subsequently refined to an exact representation by Rademacher in 1936, which provides a direct and computationally efficient expression for the single value P(n)."

Diese Formel für \(~P(n)~\) ist eine asymptotische Summenformel. Man kann sie mit etwas Wohlwollen als "geschlossene Formel" akzeptieren, weil sie die Berechnung von \(~P(n)~\) "direkt" und exakt und mit endlich vielen Summanden erlaubt. Damit ist sie von ganz anderer Natur als konkurrierende Verfahren mittels Rekursion oder erzeugender Funktion, bei denen der Algorithmus die Gesamtheit der Werte \(~P(1)~...~P(n)~\) ausgibt. Die Anzahl der Summanden in der Formel ist nicht konstant, daher mag man die Bezeichnung der Formel als "geschlossen" in Frage stellen.  -  Ganz unabhängig von diesen Bezeichnungsfragen stellt die Formel von Hardy-Ramanujan-Rademacher einen Meilenstein der Kombinatorik, der Zahlentheorie und der Numerik dar.

Die ersten Glieder der Folge \(~P(n)~\) findet man in OEIS A000041 .

Nun zur Formel selbst: Wir geben direkt die "rechnergeeignete" Umformung an (Johansson [2] (1.4), (1.5), (2.1)) :

\[\text{(1)}~~~P(n)=~\left\langle\frac{4}{24n-1} \sum_{k=1}^{k_n} A_{nk}\cdot U_{nk}\right\rangle\] \(\langle..\rangle\) steht für ganzzahlige Rundung. \[A_{nk}=\sum_{L=0 \atop {n+L(3L+1)/2~=~0~\text{mod}~k}}^{2k-1} (-1)^L~\text{cos}\left(\frac{6L+1}{6k} \pi \right)\] \[U_{nk}=\text{cosh}~C_{nk} - \frac{\text{sinh}~C_{nk}}{C_{nk}}~~~~~\text{mit}~~~~~C_{nk}=\frac{\pi~\sqrt{24n-1}}{6k}\] \(k_n\) kann als \(~\text{max}~\{~9,~\lceil\sqrt{n}~\rceil~\}~\) genommen werden; mit weniger Summanden kommt man aus, wenn man Johansson [2] (1.8) anwendet: \[\text{(2)}~~~\text{Wähle kleinstes}~~k_n~~\text{mit}~~\frac{44~\pi^2}{225 \sqrt{3~k_n}}+\frac{\pi}{75}~\sqrt{\frac{2~k_n}{n-1}}~~\text{sinh}~\left(\frac{\pi}{k_n}~\sqrt{\frac{2n}{3}}\right)\lt\frac{1}{2}\] Wegen \(~\lim_{x \to \infty} (a/x +b~x~\text{sinh}(c/x^2))=0~\) findet man mit \(x=\sqrt{k_n}~\) schnell ein passendes \(~k_n\).

In der Summe für \(~A_{nk}~\) kommen wegen der Einschränkung \(~n+L(3L+1)/2 = 0 ~\text{mod}~k~\) nur  ~\(\sqrt{k}~\) Summanden vor (Johansson [2] 2.1); die zugehörigen \(~L~\) werden weiter unten im Programm (Liste  t2 ) als  "erlaubte L"  bezeichnet.

Im Programm wird die Summation nur über Summanden \(~k~\) durchgeführt, bei denen \(~A_{nk}~\) nicht die leere Summe  (\(=0~\))  ist, denn es kommt häufig vor, dass es keine erlaubten \(~L~\) gibt. So ist im Programm  "benötigte k" (Liste  t3 ) zu verstehen.


Die Hardy-Ramanujan-Rademacher-Formel ist in verschiedenen Programmiersprachen implementiert worden. Johansson macht in [2] detaillierte Angaben zu Optimierung, Komplexität und numerischer Evaluation des Algorithmus. Er zeigt, dass eine Berechnung von \(~P(10^{15})~\) schon 2012 in etwa einer halben Stunde Rechenzeit möglich war. So ist es denn auch nicht verwunderlich, dass der Mathematica-Befehl  PartitionsP  erstaunlich schnell ist (wir wissen nichts über den verwendeten Algorithmus).  -  Diese Angaben lassen vermuten, dass auch eine einfache, nicht-optimierte Programmierung von (1) brauchbare Ergebnisse liefert, d.h. dass damit \(~P(n)~\) für "sehr große" (also weit über "Alltagsbedürfnisse" hinausgehende) \(~n~\) schnell berechnet werden kann.

Um einen Vergleich mit den optimierten Programmen in [2] ch. 4 zu ermöglichen, wird hier Mathematica 7 eingesetzt. Es soll gezeigt werden, wie (1) direkt in ein kurzes, leicht durchschaubares Programm umgesetzt werden kann  -  fernab aller Effizienzerwägungen der theoretischen Informatik:

(* Hardy-Ramanujan-Rademacher: Berechnung von P(n) *)
n = 53; (* n > 1 *)   kn = 9; (* Anzahl Summanden *)

(* Liste t2 der erlaubten L : *)
t1 = Table[Mod[n+L(3
 L+1)/2,k],{k,kn},{L,0,2k-1}];
t2 = Table[Flatten[Position[t1[[k]],0]]-1,{k,kn}];

(* Liste t3 der benötigten k : *)
anzt2 = Table[Length[t2[[k]]],{k,kn}];
t3 = Complement[Range[kn],Flatten[Position[anzt2,0]]];

A = Table[Sum[(-1)^t2[[k]][[j]]
 Cos[(6 t2[[k]][[j]]+1) Pi/(6 k)],{j,anzt2[[k]]}]/.k->t3[[i]],{i,Length[t3]}];
c = Pi
 Sqrt[24 n-1]/(6 k);
U = Table[Cosh[c]-Sinh[c]/c/.k->t3[[i]],{i,Length[t3]}];
pn = Round[(4/(24
 n-1)) Plus@@(A U)]

329931


Das Ergebnis \(~P(n)=329.931~\) wird ohne bemerkbare Wartezeit ausgegeben.

Falls man für große \(~n~\) mit weniger als \(~\text{max}~\{~9,~\lceil\sqrt{n}~\rceil~\}~\) Summanden auskommen möchte, setzt man (2) in Mathematica um:

(* Bestimmung kn *)
n
 = 53;
Do[If[44
 Pi^2/(225 Sqrt[3 N1])+Pi Sqrt[2 N1] Sinh[Pi Sqrt[2 n]/(N1 Sqrt[3])]/(75 Sqrt[n-1])<.5, kn = N1; Break[]], {N1,4 n}]
kn

9


Das Ergebnis \(~k_n=9~\) wird natürlich ohne bemerkbare Wartezeit ausgegeben; aber auch noch für \(~n=10^{12}~\) dauert die Berechnung von \(~k_n~\) nur wenige Sekunden.


Das Programm wird leichter durchschaubar, wenn man sich die Tabellen  t1, t2, t3  ausgeben lässt.

t1  enthält für jeden Summanden  k  die Werte von  Mod[n+L(3L+1)/2,k] :

t1

In  t2  stehen für jedes  k  diejenigen  L ,  die in  t1  den modulo-Wert  0  haben; darunter steht  t3 ,  also die Menge der  k  aus  t2 ,  bei denen die  L  nicht die leere Menge bilden:

t2, t3


Dass dieses Programm für praktische Zwecke hinreichend schnell ist, soll im Folgenden gezeigt werden. Da die Programmläufe auf einem handelsüblichen Notebook Baujahr 2017 durchgeführt wurden, ist mit modernen Rechnern und Mathematica-Versionen sicherlich eine noch höhere Geschwindigkeit zu erwarten.

Die Laufzeit steigt keineswegs monoton mit \(~n~\) an. Diese Beobachtung korrespondiert mit den sehr schwankenden Anzahlen der im Programm verwendeten \(~k~\) und \(~L~\) (in den Tabellen  t3  und  t2 ).  Dies wird anhand eines Beispiels deutlich:

    n    Stellen P(n) sec  kn   Anz. k  Anz. L
7.000.000   2.940     23   942   463    2226
7.000.001   2.940     13   942   289    1098
7.000.002   2.940     20   942   430    1900
7.000.003   2.940     19   942   406    1776
7.000.004   2.940     19   942   436    1780
7.000.005   2.940     41   942   639    2870
7.000.006   2.940     15   942   329    1370
7.000.007   2.940     16   942   363    1578
7.000.008   2.940     18   942   404    1748
7.000.009   2.940     17   942   413    1604


Die drittletzte Tabellenspalte zeigt die mit dem Programm  "Bestimmung kn" (s.o.) berechnete Anzahl der Summanden in (1). Das sind schon deutlich weniger als bei der gröberen Obergrenze \(~\text{max}~\{~9,~\lceil\sqrt{n}~\rceil~\}=2646~\).  Aber auch diese \(~942~\) Summanden werden nicht alle benötigt, wie die vorletzte Spalte zeigt, die die Werte aus  t3  enthält. In der inneren Summe \(~A_{nk}~\) in (1) entfallen die meisten \(~L~\);  die Anzahl der benötigten \(~L~\),  aufsummiert über alle \(~k~\),  findet sich in der letzten Spalte.  -  Man sieht in der dritten Spalte die stark schwankende Rechenzeit bei wachsendem \(~n~\);  das ist auch nicht verwunderlich, da die Anzahl der zu berechnenden Summanden in den beiden letzten Spalten sehr volatil ist und mit der Rechenzeit hoch korreliert.

Hier ist eine Tabelle für größere \(~n~\):

      n    Stellen P(n) sec    kn     Anz. k  Anz. L
 10.000.000   3.515      26   1.110     466   2.020
 20.000.000   4.974     129   1.529     844   4.388
 30.000.000   6.094     218   1.844   1.062   5.744
 40.000.000   7.038     273   2.107   1.177   5.850
 50.000.000   7.869     336   2.336   1.226   6.080
 60.000.000   8.621     367   2.543   1.162   5.952
 70.000.000   9.312     469   2.731   1.283   7.102
 80.000.000   9.956     403   2.906   1.075   5.132
 90.000.000  10.560     616   3.070   1.386   7.652
100.000.000  11.132     600   3.224   1.276   6.652


Trotz der wesentlich größeren Sprünge in \(~n~\) wächst die Rechenzeit in der dritten Spalte nicht monoton.
In diesem Bereich \(~n=10^7 ... 10^8~\) wächst die Stellenzahl von \(~P(n)~\) wie  ~\(1,11~\sqrt{n}~\),  das war zu erwarten wegen der asymptotischen Hardy-Ramanujan-Formel\[~P(n)\approx\frac{e^{\pi~\sqrt{\frac{2~n}{3}}}}{4~\sqrt{3}~n}~~~\Rightarrow~~~\text{Stellenzahl}~P(n)\approx 1,114~\sqrt{n}-log_{10}~n-0,84~\] \(~k_n~\) wächst in diesem Bereich wie  ~\(0,3~\sqrt{n}~\).



Literatur

[1]  Hardy, G. H., Ramanujan, S., Asymptotic formulae in combinatory analysis, Proc. London Math. Soc., 2nd Series, 17 (1918), p. 75 - 115.  Reprinted in Collected Papers of Srinivasa Ramanujan, Amer. Math. Soc. (2000), p. 276 - 309

[2]  Johansson, F., Efficient implementation of the Hardy-Ramanujan-Rademacher formula, LMS J. Comp. Math., 15 (2012), p. 341 - 359

[3]  Rademacher, H., On the partition function p(n), Proc. London Math. Soc., 2nd Series, 43/4 (1937), p. 241 - 254


Kommentare sind willkommen.



Stand 2020-08-13


Inhalt Blog   |    voriger Eintrag


Manfred Börgens   |   Zur Leitseite                                     Datenschutz







ab 2020-06-25
www.boergens.de/manfred/blog/blog005.htm