Matroids Matheplanet Forum Index
Moderiert von Spock mire2
Mathematische Software & Apps » Mathematica » Suche schnellsten Algorithmus für Tribonacci(10^7); Folge A000073
Druckversion
Druckversion
Antworten
Antworten
Autor
Universität/Hochschule Suche schnellsten Algorithmus für Tribonacci(10^7); Folge A000073
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Themenstart: 2021-03-08


Hallo zusammen,

eigene Recherchen haben ergeben, dass der schnellste mir bekannte Algorithmus für ganze Tribonacci Zahlen
( mathworld.wolfram.com/TribonacciNumber.html
und oeis.org/A000073 )
etwa 325 mal langsamer als der der Fibonacci Zahlen ist!
Das bekommen wir doch schneller hin - oder?

Zunächst starte ich mit dem Argument n=10^7 (10 Mio.), um nicht zu lange warten zu müssen.
Um objektiv vergleichen zu können, müssen bei der Zeitmessung auch die letzten 30 dezimalen Ziffern mit angegeben werden:
Mod[Funktion[10^7], 10^30] // Timing

Algorithmus 1:
Aus dem Forum hier
gibt es ja den bekannte Weg zur Ermittlung der expliziten Funktion auch bekannt als "BINET'S FORMULA FOR THE TRIBONACCI SEQUENCE"
www.fq.math.ca/Scanned/20-2/spickerman.pdf

Algorithmus 2:
Hypergeometrische Funktion aus den beiden Summen:
Tribonacci(n) = Sum_ {i = 0 .. floor ((n - 2)/4)} ((-1)^i*     binomial (n - 2 - 3*i, i)*2^(n - 2 - 4*i)) - Sum_ {i = 0 .. floor ((n - 3)/4)} ((-1)^i* binomial (n - 3 - 3*i, i)*2^(n - 3 - 4*i));
= 2^(-3 + n) (2 HypergeometricPFQ[...

Algorithmus 3:
Reduzierung der BINET'S FORMEL + Rundung

Ergebnisse:

Berechnungszeiten für Tribonacci(10^7)
Algo 1:        10.8 s BINET'S FORMULA
Algo 2:        79.8 s hypergeometrische Funktion
Algo 3:         5.1 s BINET'S FORMULA reduziert + gerundet
--------- ab hier YMP ------------------------------------
MatrixPo        0.59s siehe Beitrag 2 + 4
TribonacciBinet 0.4 s Beitrag 19!! 
TribonacciTac   0.33 s Beitrag 17
zum Vergleich Fibonacci: 0.015625 s

Jegliche Rekursionen (aus Vorgängern) sind bei solch großen Argumenten unbrauchbar (oder stoßen schnell an RAM-Grenzen).

Kennt jemand schnellere Algorithmen oder Abkürzungen?

Oder kann man die Tribonacci-Zahlen aus den Fibonacci Zahlen berechnen?

Grüße Gerd





Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.1, vom Themenstarter, eingetragen 2021-03-08


Interessant ist der Algorithmus von Pari/GP:
Pari/GP
a(n)=([0, 1, 0; 0, 0, 1; 1, 1, 1]^n)[1, 3]
a(10^7)%10^30 ->  1.54 s      (3.37 mal schneller als Mathematica! Algo 3)
a(10^8)%10^30 -> 16 s         (4.8 mal schneller als Mathematica!)
a(10^9)%10^30 -> 3 min 26,9 s

Mathematicas
TribLinRecur[n_] := LinearRecurrence[{1, 1, 1}, {0, 0, 1}, n] // Last;

erlaubt nur Argumente bis n=10^5 bei 32 GB RAM!

Gibt es da noch einen anderen Syntax?



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.2, vom Themenstarter, eingetragen 2021-03-08


Hab's:
Der Syntax nennt sich MatrixPower:
mathematica
TribMatrixPow[n_]:=MatrixPower[{{0,1,0},{0,0,1},{1,1,1}},n-2][[3]]//Last;
Mod[TribMatrixPow[10^7], 10^30] ->  0.31 s (Timing) -> ganze Zahl RAM-Disk:  0.59 s
Mod[TribMatrixPow[10^8], 10^30] ->  4.59 s (Timing) -> ganze Zahl RAM-Disk:  9 s
Mod[TribMatrixPow[10^9], 10^30] -> 59.66 s (Timing) -> ganze Zahl RAM-Disk:128.8 s

Damit kommen wir in die Geschwindigkeitsregion von Fibonacci!

Interessant: Mathematica nutzt dafür 3 CPU-Kerne! (deshalb schneller als Pari GP)



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Fabi
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.03.2002
Mitteilungen: 4572
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.3, eingetragen 2021-03-08


2021-03-08 19:54 - hyperG in Beitrag No. 2 schreibt:
Hab's:
Der Syntax nennt sich MatrixPower:
mathematica
TribMatrixPow[n_]:=MatrixPower[{{0,1,0},{0,0,1},{1,1,1}},n-2][[3]]//Last;
Mod[TribMatrixPow[10^7], 10^30] ->  0.31 s
Mod[TribMatrixPow[10^8], 10^30] ->  4.59 s
Mod[TribMatrixPow[10^9], 10^30] -> 59.66 s

Damit kommen wir in die Geschwindigkeitsregion von Fibonacci!

Interessant: Mathematica nutzt dafür 3 CPU-Kerne! (deshalb schneller als Pari GP)

Hi,

Willst du die komplette Zahl haben, oder reichen dir von vornherein nur die letzten 30 Stellen? Zweiteres geht um Größenordnungen schneller.

vG,
Fabi


-----------------
"There would be the mathematical equivalent of worldwide rioting." (P.C.)

Willst du Hamburg oben sehen, musst du die Tabelle drehen.



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.4, vom Themenstarter, eingetragen 2021-03-08


Hallo Fabi,

Du hast recht -> ich werde es noch eindeutiger kennzeichnen:
Es geht analog zu Fibonacci(10^9) Benchmarks um das vollständige Speichern der kompletten Dezimalzahl auf einem schnellen Medium.
Da eine RAM-Disk bei diesen Größen weniger als 0.4 % ausmacht, kann man es vernachlässigen.
Was man aber nicht vernachlässigen darf, ist die Hex-> Dez Wandlung.
Ich dachte, dass dies mit Modulo 10 schon halb erledigt sei, aber auch hier haben viele Programme viel optimiert, also der richtige Code zur Geschwindigkeitsmessung:
mathematica
TribMatrixPow[n_]:=MatrixPower[{{0,1,0},{0,0,1},{1,1,1}},n-2][[3]]//Last;
startT0=AbsoluteTime[];
Export["r:\\TribMatrixPow(1Mrd).CSV",TribMatrixPow[10^9],"CSV"]; 
Print[ AbsoluteTime[]-startT0];
128.8026730 s

Das sind also 69,14 s mehr, als wenn man nur die letzten 30 Stellen ausgibt.
Ich werde die Tabelle im Beitrag 2 erweitern.

Damit ist Tribonacci im Vergleich zu Fibonacci unter mathematica bei n=10^9 nur noch 2,274 mal langsamer.
Dafür werden ja auch 55,7 Mio. Stellen mehr berechnet!

Um das mit YMP mit allen 20 Kernen noch zu beschleunigen, fehlt mir noch ein effektiver Algorithmus für die MatrixPower Funktion.
endy hatte da mal einen Mathematica Befehl, wie das "Innere"
abgearbeitet wird...




Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
tactac
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 15.10.2014
Mitteilungen: 2074
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.5, eingetragen 2021-03-08

\(\begingroup\)\(\newcommand{\sem}[1]{[\![#1]\!]} \newcommand{\name}[1]{\ulcorner#1\urcorner} \newcommand{\upamp}{\mathbin {⅋}}\)
$\begin{pmatrix}
0&1&0\\
0&0&1\\
1&1&1
\end{pmatrix}^n$ enthält in der ersten Spalte Tribonacci-Zahlen rund um den Index $n$. Die Potenz lässt sich schnell mit Square-and-Multiply berechnen.
Alternativ kann man auch im Ring $\mathbb Z[X]/(X^3-X^2-X-1)$ die Potenz $X^n$ ausrechnen.

[Die Antwort wurde vor Beitrag No.1 begonnen.]
\(\endgroup\)


Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.6, vom Themenstarter, eingetragen 2021-03-08


Hallo tactac,

genau darum geht es ja gerade, was die "abgekürzte Schreibweise" beim Potenzieren einer quadratischen Matrix eigentlich bedeutet. (Array-Summe aus Feldmultiplikationen!!!)

Hier habe ich das mal in bekannte Array-Operationen für n=2^5-1...2^5+4 aufgeschlüsselt:
mathematica
MatrixPow2[m_]:=Table[Sum[m[[i]][[k]]*m[[k]][[j]],{k,3}],{i,3},{j,3}];
a3={{0,1,0},{0,0,1},{1,1,1}};
m5=Fold[MatrixPow2[#1]&,a3,Range[5]]
m2=MatrixPow2[a3];
m7=Table[Sum[m5[[i]][[k]]*m2[[k]][[j]],{k,3}],{i,3},{j,3}]
Table[MatrixPower[{{0,1,0},{0,0,1},{1,1,1}},k][[1]]//Last,{k,2^5-1,2^5+4,1}]
{m5[[1]]//First,m5[[1]]//Last,m5[[2]]//Last,m5[[3]]//Last,m7[[2]]//Last,m7[[3]]//Last}
 
Out[119]= {{29249425,45152016,53798080},{53798080,83047505,98950096},{98950096,152748176,181997601}}
Out[121]= {{98950096,152748176,181997601},{181997601,280947697,334745777},{334745777,516743378,615693474}}
Out[122]= {29249425,53798080,98950096,181997601,334745777,615693474}
Out[123]= {29249425,53798080,98950096,181997601,334745777,615693474}

Die beiden letzte Zeilen stimmen überein,
sind also Tribonacci(2^5-1)...Tribonacci(2^5+4)

Das mit der abgekürzten Schreibweise "Ring" habe ich noch nie gehört.



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
tactac
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 15.10.2014
Mitteilungen: 2074
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.7, eingetragen 2021-03-09

\(\begingroup\)\(\newcommand{\sem}[1]{[\![#1]\!]} \newcommand{\name}[1]{\ulcorner#1\urcorner} \newcommand{\upamp}{\mathbin {⅋}}\)
2021-03-08 23:36 - hyperG in Beitrag No. 6 schreibt:
Das mit der abgekürzten Schreibweise "Ring" habe ich noch nie gehört.
\(\begingroup\)\(\newcommand{\sem}[1]{[\![#1]\!]} \newcommand{\name}[1]{\ulcorner#1\urcorner} \newcommand{\upamp}{\mathbin {⅋}}\)
Zum Rechnen damit kann man so vorgehen: Die Elemente haben die Form $a+bX+cX^2$ mit $a,b,c \in \IZ$, mit unbekanntem $X$ (sind also eigentlich Tripel ganzer Zahlen), jedoch weiß man von dem $X$, dass $X^3 = 1+X+X^2$, weswegen man Koeffizienten für die höheren Potenzen nicht braucht. Dies ergibt, wie addiert und multipliziert werden muss, und darauf basierend kann potenziert werden. Man rechnet also statt mit 9 Zahlen im Fall der Matrix nur mit 3. Das könnte durchaus schneller sein.
\(\endgroup\)


Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.8, vom Themenstarter, eingetragen 2021-03-09


Vergleich der Geschwindigkeiten:
mathematica
startT0=AbsoluteTime[];zweihoch=26;
MatrixPow2[m_]:=Table[Sum[m[[i]][[k]]*m[[k]][[j]],{k,3}],{i,3},{j,3}];
a3={{0,1,0},{0,0,1},{1,1,1}};
m5=Fold[MatrixPow2[#1]&,a3,Range[zweihoch]];//Timing
Print[AbsoluteTime[]-startT0];
Mod[m5[[1]]//Last,10^30]
Print[AbsoluteTime[]-startT0];
tartT0=AbsoluteTime[];
Mod[MatrixPower[{{0,1,0},{0,0,1},{1,1,1}},2^zweihoch][[1]]//Last,10^30] //Timing
Print[AbsoluteTime[]-startT0];
Test 1 kostenlose wolframcloud: ###################################
{26.612,Null}
26.633396
192181029491212121499150319616
26.70
{30.0164,192181029491212121499150319616}  <- hiermit hat Cloud Probleme!
37.060348
Test 2 Cloud 2. Durchlauf:  ########################################
{26.464,Null}
26.473404
192181029491212121499150319616
26.495252
{28.8386,192181029491212121499150319616} <-- Timing zeigt in der Cloud Müll!
10.006397...12.8  <-- stimmt mit Stoppuhr
 
Test 3 Version 12\MatrixPower nutzt ja beim i9 3 Kerne (offline):#####
Out[4]= {7.578125,Null}
7.8887829
Out[6]= 192181029491212121499150319616
7.9
Out[8]= {0.046875,192181029491212121499150319616} <-- 0.04 s für internen Hex Wert
10.8724600  <-- stimmt mit Stoppuhr

MatrixPower ist gegenüber der "Sum-Nachprogrammierung" etwa 125 mal schneller!
Mit etwas If-Optimierung der Sum-Schleife (doppelt vorkommende Werte kopieren -> statt 27 Mul nur noch 21 Mul) konnte ich die 7.9 s auf 5.9 s drücken, was aber immer noch 95 mal langsamer als der gut optimierte Mathematica Befehl (3 Threads) ist.




Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.9, vom Themenstarter, eingetragen 2021-03-09


2021-03-09 00:00 - tactac in Beitrag No. 7 schreibt:
2021-03-08 23:36 - hyperG in Beitrag No. 6 schreibt:
Das mit der abgekürzten Schreibweise "Ring" habe ich noch nie gehört.
... statt mit 9 Zahlen im Fall der Matrix nur mit 3. ...

Also bei mir hat eine "Matrix-Multiplikation" eines 3*3-Arrays:
Table[Sum[m[[i]][[k]]*m[[k]][[j]],{k,3}],{i,3},{j,3}]
3*3*3= 27 Multiplikationen
Für die beiden hinteren Felder (i>1) kann man je 3 mal kopieren statt multiplizieren, was 27-2*3= 21 Multiplikationen macht. (und meine Geschwindigkeitssteigerung um Faktor 1.34 praktisch zeigt)

Wie man das bei Mio-stelligen Zahlen mit
X[...]³=1+X[...]+X[...]² noch weiter vereinfachen will, verstehe ich nicht.
Gibt's da keine LINKs oder praktisches Zahlenbeispiel statt "abgekürzte Schreibweise" ?



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
tactac
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 15.10.2014
Mitteilungen: 2074
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.10, eingetragen 2021-03-09

\(\begingroup\)\(\newcommand{\sem}[1]{[\![#1]\!]} \newcommand{\name}[1]{\ulcorner#1\urcorner} \newcommand{\upamp}{\mathbin {⅋}}\)
2021-03-09 19:14 - hyperG in Beitrag No. 9 schreibt:
\(\endgroup\)
\(\begingroup\)\(\newcommand{\sem}[1]{[\![#1]\!]} \newcommand{\name}[1]{\ulcorner#1\urcorner} \newcommand{\upamp}{\mathbin {⅋}}\)2021-03-09 00:00 - tactac in Beitrag No. 7 schreibt:
2021-03-08 23:36 - hyperG in Beitrag No. 6 schreibt:
Das mit der abgekürzten Schreibweise "Ring" habe ich noch nie gehört.
\(\begingroup\)\(\newcommand{\sem}[1]{[\![#1]\!]} \newcommand{\name}[1]{\ulcorner#1\urcorner} \newcommand{\upamp}{\mathbin {⅋}}\)
... statt mit 9 Zahlen im Fall der Matrix nur mit 3. ...
\(\endgroup\)

Also bei mir hat eine "Matrix-Multiplikation" eines 3*3-Arrays:
Table[Sum[m[[i]][[k]]*m[[k]][[j]],{k,3}],{i,3},{j,3}]
3*3*3= 27 Multiplikationen
Für die beiden hinteren Felder (i>1) kann man je 3 mal kopieren statt multiplizieren, was 27-2*3= 21 Multiplikationen macht. (und meine Geschwindigkeitssteigerung um Faktor 1.34 praktisch zeigt)

Wie man das bei Mio-stelligen Zahlen mit
X[...]³=1+X[...]+X[...]² noch weiter vereinfachen will, verstehe ich nicht.
Gibt's da keine LINKs oder praktisches Zahlenbeispiel statt "abgekürzte Schreibweise" ?

\(\begingroup\)\(\newcommand{\sem}[1]{[\![#1]\!]} \newcommand{\name}[1]{\ulcorner#1\urcorner} \newcommand{\upamp}{\mathbin {⅋}}\)

Ich habe doch alles beschrieben: Man rechnet mit Tripeln $(a,b,c)$, die für $a+bX+cX^2$ stehen. Zusammen mit der Festlegung, dass $X^3 = 1+X+X^2$ sein soll, ergibt sich, wie zwei solcher Tripel zu multiplizieren sind. (Das ist so ähnlich wie mit komplexen Zahlen, nur, dass dort mit Paaren der Form $a+bX$ und der Festlegung $X^2 = -1$ gerechnet wird.)
Nämlich: $(a,b,c) \cdot (d,e,f) = (k_0,k_1,k_2)$ mit:
* $k_0 = ad + k_3$,
* $k_1 = ae + bd + k_3 + k_4$,
* $k_2 = af + be + cd + k_3 + k_4$,
* $k_3 = bf + ce + k_4$,
* $k_4 = cf$.

So. $X$ selbst ist das Tripel $(0,1,0)$. Jedes $X^n$ ist ein Tripel $(a_n,b_n,c_n)$, das man effizient mit square-and-multiply auf der Basis der Multiplikation der Tripel (in der 9 Multiplikationen in $\IZ$ vorkommen) ausrechnen kann.
Die so entstehenden Folgen $(a_n)_{n\in\IN}$ und $(c_n)_{n\in\IN}$ sind Folgen von Tribonacci-Zahlen (nicht beim Index 0 beginnend, sondern leicht verschoben). Z.B. ist $T_n = c_{n+1}$, wenn $T_0 = 0, T_1 = 1, T_2 = 1$.
\(\endgroup\)


Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.11, vom Themenstarter, eingetragen 2021-03-10



Input n= 2^28
Mathematica 3 + 1 Kern: 28.9 s

YMP  1 Kern & 21-Mul-Algorithmus bin ich bei:

also 17.5 s

Das geht noch schneller...


Also tactac,
immer noch zu viele abgekürzte Symbole,
zu viele Sprünge,
unfertig umgestellte Formeln (X²=-1),
aus dem Nichts auftauchende k3 , k4...
weder sehe ich Input noch Output, wie es sich für einen Algorithmus mit "roten Faden" gehört.

Ich zeige das mal morgen am Beispiel von
Input: n=2^3
Algo: MatrixPower[{{0, 1, 0}, {0, 0, 1}, {1, 1, 1}}, 2^3][[1]] // Last
Output: 24



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Delastelle
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 17.11.2006
Mitteilungen: 1670
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.12, eingetragen 2021-03-10


Hallo,

...
um tactac zu verstehen sollte man erst mal einen Schluck
Zaubertrank nehmen. (Avatar von tactac)
...

Viele Grüße
Ronald



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
DerEinfaeltige
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 11.02.2015
Mitteilungen: 2851
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.13, eingetragen 2021-03-10


tactac's "Mathematisch" übersetzt in "Python für Arme" sähe ganz naiv so aus:
Python
def mult(A,B):
    a1,a2,a3 = A
    b1,b2,b3 = B
    k4 = a3*b3
    k3 = a2*b3 + a3*b2 + k4
    k2 = a1*b3 + a2*b2 + a3*b1 + k3 + k4
    k1 = a1*b2 + a2*b1 + k3 + k4
    k0 = a1*b1 + k3 
    return (k0,k1,k2)
 
def pow(T,n):
    if n==1:
        return T
    if n%2==0:
        return pow(mult(T,T),n//2)
    return mult(pow(T,n-1),T)
 
def Tribonacci(n):
    if n == 0:
        return 0
    T = (0,1,0)
    a,b,c = pow(T,n+1)
    return c




-----------------
Why waste time learning when ignorance is instantaneous?
- Bill Watterson -



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Kitaktus
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 11.09.2008
Mitteilungen: 6850
Wohnort: Niedersachsen
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.14, eingetragen 2021-03-10


Ein Bemerkung:

Möchte man T(n) von vornherein nur modulo einer bestimmten Zahl M berechnen, so wächst die notwendige Rechenzeit nur logarithmisch mit n, also linear mit der Stellenzahl(!) von n(*).

Ist man an einem exakten Ergebnis interessiert (also ohne modulo), dann wächst die Zahl der Stellen, die dafür berechnet werden müssen linear mit n, was zu einem Gesamtaufwand von $O(n\log n)$ führt.

Algorithmisch ist das Problem sehr leicht, für die Berechnung von T(109) reichen ein paar hundert Multiplikationen und Additionen aus. Das einzige Problem ist die Länge der beteiligten Zahlen.

(*) Modulo 106 rechnet Matlab T(109) in unter einer Millisekunde aus. Größere Module sind problematisch, weil Matlab (zumindest in meiner Version) kein Zahlenformat für lange Ganzzahlen hat.

[Die Antwort wurde nach Beitrag No.12 begonnen.]



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
tactac
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 15.10.2014
Mitteilungen: 2074
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.15, eingetragen 2021-03-10


Für Leute, die Python nicht so gut lesen können, vll. noch eine Haskell-Version:
Haskell
data T n = T n n n deriving (Eq, Ord, Show)
 
instance Num n => Num (T n) where
    -- der Vollständigkeit halber auch mit + und -, auch wenn
    -- das für's vorliegende Problem nicht benötigt wird
    fromInteger n = T (fromInteger n) 0 0
    T a b c + T d e f = T (a+d) (b+e) (c+f)
    T a b c - T d e f = T (a-d) (b-e) (c-f)
    T a b c * T d e f = T k0 k1 k2 where
        k0 =       a*d       + k3
        k1 =    a*e + b*d    + k3 + k4
        k2 = a*f + b*e + c*d + k3 + k4
        k3 =    b*f + c*e         + k4
        k4 =       c*f
    abs = undefined
    signum = undefined
 
tribo :: Int -> Integer
tribo n = let T _ _ c = T 0 1 0 ^ (n+1) in c


[Die Antwort wurde nach Beitrag No.13 begonnen.]



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.16, vom Themenstarter, eingetragen 2021-03-10


2021-03-10 09:24 - DerEinfaeltige in Beitrag No. 13 schreibt:
tactac's "Mathematisch" übersetzt in "Python für Arme" sähe ganz naiv so aus:
Python
def mult(A,B):
    a1,a2,a3 = A
    b1,b2,b3 = B
    k4 = a3*b3
    k3 = a2*b3 + a3*b2 + k4
    k2 = a1*b3 + a2*b2 + a3*b1 + k3 + k4
    k1 = a1*b2 + a2*b1 + k3 + k4
    k0 = a1*b1 + k3 
    return (k0,k1,k2)
...


Oh, danke DerEinfaeltige,
jetzt sieht man doch genau, dass der tactac-Array-Mul-Algorithmus mit 9 "echten Multiplikationen" auskommt. Danke tactac.

In Mathematica sind das 3 Zeilen:
mathematica
MulTac[a_,b_]:=Module[{k4=a[[3]]*b[[3]]},k3=a[[2]]*b[[3]]+a[[3]]*b[[2]]+k4;ret={a[[1]]*b[[1]]+k3,a[[1]]*b[[2]]+a[[2]]*b[[1]]+k3+k4,a[[1]]*b[[3]]+a[[2]]*b[[2]]+a[[3]]*b[[1]]+k3+k4};ret];
PowTac[T_,n_]:=If[n==1,T,If[Mod[n,2]==0,PowTac[MulTac[T,T],Quotient[n,2]],MulTac[PowTac[T,n-1],T]]];
TribonacciTac[n_]:=Module[{aT},aT=PowTac[{0,1,0},n];aT[[3]]];
Probe:
Table[{TribonacciTac[k],MatrixPower[{{0,1,0},{0,0,1},{1,1,1}},k][[1]]//Last},{k,10}]//Grid
nach n-Offset Korrektur stimmen beide:
0	0
1	1
1	1
2	2
4	4
7	7
13	13
24	24
44	44
81	81

Nun kann ich das nach YMP wandeln (sehr Hardware-nah: FFT-Mul, AVX2, Multitasking) & mit Mio.-stelligen Zahlen testen...

Bei meinem gestrigen Programm (21 Multiplikationen) konnte ich nun auch die 20 Kerne zuschalten & die Berechnungszeit der 71 Mio. stelligen Zahl von 17.5 s auf 4.4 s reduzieren!


Mit 9 Mul sollte das noch zu halbieren sein...



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.17, vom Themenstarter, eingetragen 2021-03-10


Tatsächlich, YMP kann mit dem 9-Mul-Algorithmus die Rechenzeit halbieren:


Tribonacci(2^28)
Mathematica MatrixPow      28.9 s  3 Threads
Mathematica TribonacciTac  26.6 s  1 Thread
YMP Tribonacci2h 21 Mul     4.4 s 20 Threads
YMP TribonacciTac 9 Mul     2.3 s 20 Threads

YMP TribonacciTac(10^7) braucht nur 0.33 s -> Startbeitrag.

YMP TribonacciTac(10^9) -> 35 s  264649443 Bytes


Super!

Vielen Dank an alle Mitwirkenden!

Grüße Gerd



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.18, vom Themenstarter, eingetragen 2021-03-10

Benchmark Tribonacci(10^9) incl. RAM-Disk
Mathematica 9 Mul       297.9 s 1 Thread
Mathematica MatrixPower 128.8 s 3 Threads
YMP 9-Mul                35.0 s 20 Threads
YMP BINET's Formula      33.48s 20 Threads Beitrag 19

Während bei 2^28 der tactac's 9 Mul Algo noch vorn lag, holte Mathematicas 3 Thread Befehl MatrixPower bei 10^9 auf.

Es kann auch sein, dass MatrixPower durch folgende Optimierung schneller wurde: en.wikipedia.org/wiki/Addition-chain_exponentiation



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.19, vom Themenstarter, eingetragen 2021-03-12


Nur weil Mathematica bei BINET's Formel so langsam ist,
muss das bei YMP nicht auch so sein!

Also habe ich ALLES von Ganzzahl auf Gleitkommazahl umgestellt.
Wurzel-Iteration gab es schon für 1/Sqrt(x) -> also Sqrt(x)=x*1/Sqrt(x)
(es geht bei Rekorden möglichst darum, die langsame Division zu vermeiden)

Für 3. Wurzel = x^(1/3) habe ich aber keine schnelle Iteration gefunden -> also die altbekannte Newton-Iteration.

Für n=10^7 kam tatsächlich 0,4 s heraus!

Das macht optimistisch -> also n=10^9



Rekord! 2 s schneller!

Der Wert stimmt exakt: 5574879143...871128684941792256.00

UNGLAUBLICH !!!


Falls sich noch eine schnellere Iteration für die 3. Wurzel findet als
r0=x/2.5;
r:=2*r/3+x/(3*r*r)
lim r = x^(1/3)

...oder jemand eine Iteration ohne die beiden Divisionen kennt...

...kann sich ja mal melden.

Grüße Gerd



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Slash
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 23.03.2005
Mitteilungen: 8420
Wohnort: Sahlenburg (Cuxhaven)
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.20, eingetragen 2021-03-13


Kenne jetzt den Thread nicht... Hattest du schon die OEIS gecheckt?

A000073

Da stehen viele Formeln und auch Programme. PARI ist doch auch sehr flott bei solchen Geschichten.

Hier noch etwas zur n-ten Zahl: Simon Plouff

Gruß, Slash


-----------------
Science can't see what it doesn't have the language to describe.



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 03.02.2017
Mitteilungen: 1454
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.21, vom Themenstarter, eingetragen 2021-03-13


2021-03-13 00:20 - Slash in Beitrag No. 20 schreibt:
Kenne jetzt den Thread nicht... Hattest du schon die OEIS gecheckt?

A000073

Da stehen viele Formeln und auch Programme. PARI ist doch auch sehr flott bei solchen Geschichten.

Hier noch etwas zur n-ten Zahl: Simon Plouff

Gruß, Slash

Man merkt, dass Du fast nichts gelesen hast:
- A000073  Startbeitrag
- PARI -> Beitrag Nr. 1
- daraus entwickelte sich die schnellere MatrixPower
- plouffe.fr -> hat noch 4 3. Wurzeln 2000 Stellen -> ich 2 und gehe gegen 1 Mrd Stellen

Meine letzte Frage ging Richtung: 3. Wurzel Iteration,
denn da hat der Weltmeister bereits tolle Ideen zur 2. Wurzel. (Beitrag 19)

Solche 1/x * Faktor Iteration habe ich aber noch nicht gefunden.



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
hyperG hat die Antworten auf ihre/seine Frage gesehen.
Neues Thema [Neues Thema] Antworten [Antworten]    Druckversion [Druckversion]

 


Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2021 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]