| |
| Autor |
 Algorithmen für Bellsche Zahlen (schnell und reell) |
|
hyperG
Senior 
Dabei seit: 03.02.2017
Mitteilungen: 1677
 |  Themenstart: 2022-01-22
|
Hallo zusammen,
es geht um die Bellschen Zahlen: 1, 2, 5, 15, 52, 203,...
Quellen für Algorithmen gibt es ja viele wie:
https://de.wikipedia.org/wiki/Bellsche_Zahl
https://oeis.org/A000110
https://functions.wolfram.com/IntegerFunctions/BellB/
https://mathworld.wolfram.com/BellNumber.html
zunächst schnellste exakte Integer Algorithmus
Argument 500 ergibt eine 844stellige Zahl 1606...94229457772,
dessen Berechnungsdauer man schön vergleichen kann:
\sourceon mathematica
BellB[500] -> 0.0156 s
Ceiling[Sum[k^n/k!, {k, 1, 2*n}]/E] -> 0.078 s
(* https://functions.wolfram.com/IntegerFunctions/BellB/06/01/0002/ *)
Sum[(k^n/k!) Sum[(-1)^j/j!, {j, 0, n - k}], {k, 1, n}] -> 1.39 s
(* Summe der Stirling2 Zahlen, also n=500 *)
n!*SeriesCoefficient[Series[E^(E^t - 1), {t, 0, n}], n] -> 2.328 s
HypergeometricPFQ[ConstantArray[2,n-1],ConstantArray[1,n-1],1]/E -> >1 min
(* wegen ConstantArray nur ganzzahlige Argumente machbar! *)
\sourceoff
Die meisten werden nur Algorithmus 3 "Summe der Stirling2 Zahlen" kennen, aber der landet nur auf Platz 3.
Frage 1:
Welchen Algorithmus könnte Mathematica für BellB verwendet haben? Oder wurde einfach nur Algo 2 optimiert (z.B. bereits berechnete Stirling2 Zahlen in Array zwischenspeichern)?
Frage 2:
Kennt jemand weitere Algorithmen, die auch für Argumente größer 1000 brauchbar sind?
Algorithmen für reelle Argumente
Zwar wird es hier viele geben, die analog der Fakultät wieder schreien "gibt's nicht", -> aber genau da wird es für mich interessant (analog der Gammafunktion).
Alle Näherungsfunktionen die ich gefunden habe, sind sehr ungenau und erst für sehr große Argumente mit wenigen richtigen Dezimalstellen...
(orange im Diagramm)
a) exakte Integrale:
https://functions.wolfram.com/IntegerFunctions/BellB/07/01/0001/MainEq1.gif
Leider sehr aufwendig zu berechnen, wenn man mehr als 3 exakte Nachkommastellen haben möchte (rot im Diagramm)
b) unendliche Summe:
https://functions.wolfram.com/IntegerFunctions/BellB/06/03/0002/MainEq1.gif
Leider extrem große Zwischenergebnisse -> schlechte Konvergenz (grün)
Hier der Vergleich der 4 Algorithmen (blau ist nur eine Interpolation der Integer-Schnittstellen):
https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_BellB_reell_Diagramm.PNG
Frage 3:
kennt jemand brauchbare Algorithmen, die schnell & genau sind?
- hypergeometrische Funktionen {jedoch reelle Argumente}
- genauere Näherungsfunktionen {analog der Gammafunktion, wo die Stirlingsche Näherungsfunktion per weiterer Bernoulli-Terme beliebig genau gemacht werden kann}
- lässt sich aus der "generating Function" {https://mathworld.wolfram.com/BellNumber.html (11) } was machen?
Grüße Gerd
|
 Profil
|
hyperG
Senior 
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.1, vom Themenstarter, eingetragen 2022-01-22
|
Unter
https://arxiv.org/pdf/2110.01156.pdf (24)
fand ich eine sehr viel bessere & schnelle Näherungsformel:
\sourceon nameDerSprache
n ergebnis richtige Digits
3 5.017 2
6 203.2 3
500 1606073.. 6 ( < 10 ms !)
\sourceoff
|
Profil
|
cramilu
Aktiv
Dabei seit: 09.06.2019
Mitteilungen: 1485
Wohnort: Schwäbischer Wald, seit 1989 freiwilliges Exil in Bierfranken
 | Beitrag No.2, eingetragen 2022-01-22
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.3, vom Themenstarter, eingetragen 2022-01-22
|
\quoteon(2022-01-22 20:09 - cramilu in Beitrag No. 2)
... folgenden historischen Beitrag toll:
An Asymptotic Formula for the Bell Numbers
...
\quoteoff
Hi cramilu,
absolut identisch, nur dass in 3.8 ein R und in (24) ein w verwendet wurde:
https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_Bell_N_herung_identisch.png
für die LambertW Funktion. Leider nicht so gut wie die Stirlingsche Näherungsformel für Fakultät.
|
Profil
|
endy
Senior
Dabei seit: 10.01.2011
Mitteilungen: 3267
 | Beitrag No.4, eingetragen 2022-01-23
|
Hallo.
Zu Frage 1 :
\sourceon mathematica
Clear @ "Global`*"
a[n_] := BellB[n]
b[n_] := Sum[StirlingS2[n, k], {k, 0, n}]
\sourceoff
Ich vermute es ist der Algorithmus b.
Gruss endy
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.5, vom Themenstarter, eingetragen 2022-01-23
|
\quoteon(2022-01-23 00:02 - endy in Beitrag No. 4)
...
\sourceon mathematica
b[n_] := Sum[StirlingS2[n, k], {k, 0, n}]
\sourceoff
...
\quoteoff
Hi endy,
ja, das hatte ich ja auch vermutet. Und natürlich hat mathematica ja
die innere Summe bereits durch eine eigene Funktion "Stirling2" optimiert.
Bei Messungen (mit AbsolutTime[]-ta statt //Timing, da mathematica da gern tricks) gibt es aber immer noch gewaltige Unterschiede:
n=2000;
BellB -> 10^-8 s
Sum[StirlingS2 -> 0.008 s
Gab es nicht mal ein Befehl in Mathematica, wo man den inneren Algorithmus sichtbar machen konnte?
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.6, vom Themenstarter, eingetragen 2022-01-23
|
Ha, Mathematica arbeitet mit Cache!
Erster Lauf mit
n=4000
BellB -> 2.5 s
Sum[StirlingS2 -> 0.002 s
2. Lauf
BellB -> 10^-8 s
Sum[StirlingS2 -> 0.0015 s
Für ein Ergebnis mit 9707 Stellen sehr gut!
|
Profil
|
endy
Senior
Dabei seit: 10.01.2011
Mitteilungen: 3267
| Beitrag No.7, eingetragen 2022-01-23
|
Hallo.
@Nr. 5:
Der Befehl ist Trace und damit ist Frage Nr. 1 gelöst
\sourceon mathemtica
c[n_] := Total[Table[StirlingS2[n, k], {k, 0, n}]]
\sourceoff
Gruss endy
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.8, vom Themenstarter, eingetragen 2022-01-23
|
Danke endy!
Es kann durchaus sein, dass die Tabelle zunächst in einem Array abgelegt wird, welches sich dann minimal schneller aufaddieren lässt.
Ab Argument 20000 werden die Unterschiede besonders deutlich.
Wie immer übergebe ich die Ergebnisse dezimal in eine schnelle RAM-Disk,
um mit "csv statt w..." die Wandlung von Hex to Dezimal zu erzwingen aber die Nachteile von langsamer Festplatte & langsamer Ausgabe zu minimieren:
\sourceon mathematica
max1 = 20000; ta = AbsoluteTime[];
Export["r:\\Bell" <> ToString[max1] <> "Sum.csv", Sum[StirlingS2[max1, k], {k, 0, max1}], "CSV"] ; (*oder BellB[max1] *)
AbsoluteTime[] - ta
(*Ergebnisse für die 60551stellige Zahl: *)
a) nach Erststart ohne Cache:
- BellB -> 398.78 s
- Sum S2 -> 399.52 s
b) nochmalige Berechnung mit Zwischenergebnisse im Cache:
- BellB -> 0.078 s
- Sum S2 -> 87.198 s
c) beste Näherungsformel:
0.6 s, aber nur die ersten 9 Ziffern stimmen (schlechte asympt. Annäherung)
\sourceoff
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.9, vom Themenstarter, eingetragen 2022-01-23
|
Neuer Rekord mit dem Ceiling-Algorithmus durch kleine Optimierung:
\sourceon mathematica
Bsum3[n_] :=
Block[{$MaxExtraPrecision = Ceiling[n^1.14] + 1000},
Module[{p = 1, k = 1, end = n*7/6, sum = 0},
While[k <= end, sum += k^n/p; k++; p *= k]; Ceiling[sum/E]]];
bei Argument 20000 -> 285.427 s -> also etwa 114 s schneller!!
\sourceoff
Mir ist aufgefallen, dass die Summanden bei
N[(n*7/6)^n/(n*7/6)!, 22] kleiner als 10^-4400 sind und durch das Aufrunden nichts mehr zur Ganzzahl beitragen!
Die Aufsummierungsmenge wird kleiner & die Terme sind keine S2 Zahlen, sondern Brüche aus Potenz & Produkt.
Optimierung:
- MaxExtraPrecision könnte noch enger an der wirklich nötigen Genauigkeit angepasst werden
- Abbruch nur bei kleinen Argumenten bei 2*n -> je größer, um so dichter kann man an n; *7/6 oder 8/7 , 9/8,...?
- oder statt festen Endwert, besser als Abbruch die aktuelle Termgröße z.B. 10^(-n/5)
- statt E auf fast 61000 Stellen zu berechnen (oder geschummelt vorab zu speichern), müsste man diese irrationale Zahl in der Schleife parallel mit berechnen...
Grüße Gerd
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.10, vom Themenstarter, eingetragen 2022-01-25
|
Und nun mit ymp mit 20 Threads:
Argument 20000:
9,9 s
Also 398.78 s/ 9,9s = 40,28 mal schneller als Mathematica!!!
Wow!
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.11, vom Themenstarter, eingetragen 2022-01-25
|
Nochmals optimiert und größere Argumente: n=50000
\sourceon nameDerSprache
Mathematica/ YMP =
12847.7 s / 99.65 s = 129 mal schneller!
\sourceoff
Gewaltiger Unterschied, ob man 1,5 min oder über 3,5 Stunden auf das Ergebnis wartet!!!
Alle 168862 Stellen stimmen überein!
Die E-Berechnungszeit: 0.007 s spielt dabei keine Rolle.
Die CPU ist noch nicht voll ausgelastet -> da geht noch mehr...
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.12, vom Themenstarter, eingetragen 2022-01-26
|
https://matheplanet.com/matheplanet/nuke/html/uploads/b/47407_Bell50000_40s.PNG
Also
\sourceon BellB[50000]
Mathematica/ YMP =
12847.7 s / 39.97 s = 321 mal schneller als Mathematica !!
\sourceoff
Überprüfen kann man die letzten Stellen mit
\sourceon mathematica
(*
BellB[prime^m + k] mod prime = (m*BellB[k] + BellB[k + 1]) mod prime
Probe: *)
Table[{p = Prime[k], m = 6 - k, Mod[(b = BellB[p^m + k]), 1000],
Mod[b, p], Mod[BellB[1 + k] + BellB[k]*m, p]}, {k, 1, 5}] // Grid
2 5 147 1 1
3 4 474 1 1
5 3 700 0 0
7 2 132 5 5
11 1 147 2 2
\sourceoff
Da hier nur 1 Summe (statt wie bei der Summe aus Stirling2 zwei verschachtelte Summen) berechnet werden muss, wird mit steigenden Argumenten der Geschwindigkeitsgewinn gegenüber Mathematica immer größer & man kann in ganz neue Dimensionen vorstoßen!
Grüße Gerd
|
Profil
|
pzktupel
Aktiv
Dabei seit: 02.09.2017
Mitteilungen: 2180
Wohnort: Thüringen
 | Beitrag No.13, eingetragen 2022-01-27
|
Wäre natürlich schön, wenn die erzielten programmbasierten Leistungen auch mal in anderer Software eingang finden würden.
👍
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.14, vom Themenstarter, eingetragen 2022-01-27
|
Noch größere Argumente: n=100000
\sourceon nameDerSprache
Mathematica/ YMP =
180564 s / 246 s = 734 mal schneller als Mathematica !!
\sourceoff
Die Zeit für Mathematica habe ich diesmal per Näherungsfomeln berechnet, da der Faktor zwischen beiden sehr gut linear ansteigt.
Wenn der RAM ausreichen sollte, muss man also bei Mathematica für BellB[100000] {incl. Speicherung der Dezimalen Ziffern auf schneller RAM-Disc statt der geschummelten Timing-Messung} etwa
2 Tage 2 Stunden 9 min 24 s warten https://matheplanet.com/matheplanet/nuke/html/images/forum/subject/rotate.gif
|
Profil
|
hyperG
Senior
Dabei seit: 03.02.2017
Mitteilungen: 1677
| Beitrag No.15, vom Themenstarter, eingetragen 2022-01-28
|
\quoteon(2022-01-26 17:58 - hyperG in Beitrag No. 12)
...
Überprüfen kann man die letzten Stellen mit
\sourceon mathematica
(*
BellB[prime^m + k] mod prime = (m*BellB[k] + BellB[k + 1]) mod prime
Probe: *)
...
\sourceoff
...
\quoteoff
Überprüfung der letzten Ziffer einer fast 400000stelligen Zahl mit Modulo 5:
\sourceon nameDerSprache
100000 = 5^7 + 21875
BellB[10^5] mod 5 == ... 77020494207143379 mod 5
= (Bell[21875]*7 + Bell[21875 + 1]) mod5
= (... 37731710252239882 * 7 + ...14908875726268245) mod 5
= 4
\sourceoff
Egal welche Stellen noch davor stehen, (10*n+9) mod 5 ist immer 4:
\sourceon mathematica
Table[{Mod[
ToExpression[StringTake[str1, -k]]*7 +
ToExpression[StringTake[str2, -k]], 5],
Mod[ToExpression[StringTake["77020494207143379", -k]], 5]}, {k, 17}]
Out:
{{4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4,
4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}, {4, 4}}
\sourceoff
BellB[21875] könnte man nun wiederum so überprüfen, bis die Argumente klein genug sind, um das Ergebnis im Kopf zu berechnen...
Interessant wären geeignete große Primzahlen, um gleich mehrere letzte Stellen leicht auf Modulo zu überprüfen:
Argument = prime^m+k
Beispiel:
101^3 + 1 = 1030302
also
BellB[1030302] mod 101
=Mod[BellB[1]*3 + BellB[2], 101]
= (1*3+2)mod 101 = 5
Leider enden zu viele Zahlen bei Mod 101 auf 5:
\sourceon nameDerSprache
n n mod 101
101106 5
101207 5
101308 5
101409 5
101510 5
101611 5
101712 5
101813 5
101914 5
102015 5
102116 5
102217 5
102318 5
102419 5
102520 5
102621 5
102722 5
102823 5
102924 5
103025 5
103126 5
103227 5
103328 5
103429 5
...
\sourceoff
|
Profil
|
| hyperG hat die Antworten auf ihre/seine Frage gesehen. | | hyperG wird per Mail über neue Antworten informiert. |
|
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2022 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die
Nutzungsbedingungen,
die Distanzierung,
die Datenschutzerklärung und das Impressum.
[Seitenanfang]
|
Home / Seite ohne Frame
Wie unterstütze ich den Matheplaneten mit Geld?
