Numerische Bibliothek in C++

[headmyshoulder]

Hallo,

also die Idee ist denkbar simple:

Eine Bibliothek zur Lösung Standardprobleme der Numerischen Mathematik, wie

• Lösung von gewöhnlichen Differentialgleichungen
• Nullstellensuche
• Numerische Minimierung/Maximierung
• Lineare Regression
• usw.

Das wichtigste dabei ist, daß die Standardtechniken von C/C++ benutzt werden. Insbesondere soll sie ohne eigene Container und Datentypen auskommen und die Methoden der STL, wie vectoren, complex<>, ... nutzen. Ausserdem soll die Bibliothek einfach zu benutzen sein, daß auch unerfahrene C/C++ Benutzer schnell damit zurechtkommen.

Es gibt natürlich schon einige Bibliotheken die soetwas machen, aber entweder sie haben komische Lizenzvereinbarungen oder, sind recht kryptisch zu benutzen oder aber, sie haben ihre eigene (umständlichen) Container.

Das ganze Projekt soll natürlich Open Source sein, also wird es kein Geld dafür geben. Ich stell mir vor, daß daraus eine Community entsteht, welche sich dieser grossen Thematik annimmt und natürlich kann man sowas auch nicht alleine stemmen.

Viele Grüße,

headmyshoulder

[knivil]

Deine Frage oder wie darf ich deinen Beitrag verstehen? Hast du denn schon was?

[nurf]

headmyshoulder schrieb:

also die Idee ist denkbar simple:

Und eine Suchmaschine benutzen ist schwer ?
Wozu das Rad neu erfinden ?

[headmyshoulder]

knivil schrieb:

Deine Frage oder wie darf ich deinen Beitrag verstehen? Hast du denn schon was?

Prinzipiell ist das ein sehr grosses Projekt, welches ich nicht alleine stemmen kann und dafür suche ich hier Leute die daran interessiert sind mitzuarbeiten.

Ein bisschen Code ist schon da, genauer gesagt sind die 3 Solver für gewöhnliche Differentialgleichungen implementiert. Das ganze Projekt liegt auf SF

https://sourceforge.net/projects/yanl/

[headmyshoulder]

nurf schrieb:

headmyshoulder schrieb:

also die Idee ist denkbar simple:

Und eine Suchmaschine benutzen ist schwer ?
Wozu das Rad neu erfinden ?

Nee.

Ich hab mir fast gedacht daß so eine Frage kommt und natürlich ist die auch berechtigt.

Alle Bibliotheken die sowas machen, sind wirklich nicht schön. Wenn man sowas mehrmals benutzt wird man das schnell feststellen. Zum Beispiel hat der Numerical Recipes Code eine die sehr ungenehme Fortranindizierung von arrays. Davon mal davon abgesehen ist der Code auch nicht OS. Die gsl ist von der Benutzung auch nicht sehr schön, und die haben ihre eigenen Vektoren und komplexen Datentypen zwischen man denen man immer hin-und herkonvertiert.

Ich glaub schon, daß im C/C++ Bereich ein grosses Interesse an einer solchen Bibliothek mit einem einfachen Interface besteht.

[headmyshoulder]

1. Projektidee

YaNL - Yet another Numerical Library - ist eine numerische Bibliothek
für C++. Es gibt einige Open Source Bibliotheken die sich diesem Thema
annehmen, beispielsweise die fftw für Fouriertransformation oder die gsl
für den C Bereich. Viele Algorithmen in diesen Bibliotheken werden auch
in den Numerical Reciepes behandelt. Dieses Buch ist für fast alle
Anfänger auf diesem Gebiet der Quasistandard. Es gibt dazu auch
umfangreichen Sourcecode, sowohl in "reinem" C als auch in C++.

Soweit so gut, aber wozu dann eine neue Library und das Rad neu
erfinden? Das Problem mit all den oben genannten Projekten ist, dass
jedes seine eigenen Standards schafft: für Container (wie Vektoren und
Matrizen), Datentypen (complex) und Syntax. Diese Standards sind im
alltäglichen Gebrauch sehr lästig, beispielsweise wird man ständig
Vektoren der verschiedenen libs hin- und herkopieren.

Lange Rede, kurzer Sinn. Wir wollen das ändern. Wir wollen eine
numerische Bibliothek schaffen, die sich soweit wie möglich an den
Standards von C++ hält und nur Standardcontainer verwendet. Soweit ich
weiss, gibt es in der C++ Welt nichts vergleichbares mit der gsl,
Numerical Reciepes und numpy/scipy. Viele Leute würden so ein Projekt
begrüssen und wir haben diese Meinung auch von vielen Kollegen und
Programmierern gehört.

2. Inhalt der Bibliothek

Folgende Probleme sollen von der YaNL gelöst werden können:

* Numerische Lösung von gewöhnlichen Differentialgleichungen
* Numerische Integration
* Nullstellensuche
* Minimierung
* Lineare Regression
* Interpolation
* Matrizen, Matrixoperationen und Eigensysteme
* ...

Einige typische Probleme sind schon quasi standarisiert worden,
beispielsweise

* Spezialfunktionen sind in der STL TR1 enthalten
* Zufallszahlen und Verteilung in der STL TR1 und in der boost

Einige Probleme sind vielleicht auch zu speziel um die hier vernünftig
bearbeitet zu werden, beispielsweise Fast Fourier Transformations und
auch lineare Algebra, mit optimierten Operationen.

3. Was ist schon da?

Bis jetzt haben wir mit den Routinen zur Lösung von gewöhnlichen
Differentialgleichungen angefangen. Der Code dafür existiert, wobei wir
gerade dabei sind verschiedene Designs auf Geschwindigkeit und Usability
zu testen.

Wir sind noch recht weit am Anfang, viel von der Projektorganisation
fehlt auch noch, beispielsweise eine Organisation der Sourcen und
Releaseszyklen.

4. Technisches

Wir programmieren bis jetzt mit den Standardlinuxtools, wie gcc, make,
emacs, ..., Wichtig ist, daß zusätzliche Abhängigkeiten von externen
Bibliotheken nach Möglichkeit vermieden werden. Ausnahmen sind dabei
reine header only libs wie die boost oder die STL.

Der gesamte Code soll einfach zu lesen und zu benutzen sein, so daß auch
C++ Anfänger nach kurzem Studium der Beispiele sich mit YaNL
zurechtfinden. Das heisst, das Templates nur an Stellen die notwendig
sind verwendet werden oder wo Metaprogrammierung sehr grosse
Geschwindigkeitsvorteile bringt.

Das Projekt liegt bis jetzt auf Source Forge:

https://sourceforge.net/projects/yanl/

Unter was für einer Lizenz das ganze stehen soll ist noch offen, klar
ist auf jedenfall, daß YaNL Open Source sein wird.

5. Mitmachen

Um diese Projekt zu realisieren brauchen wir Programmierer die Lust und
Spass auf eim solchen Communityproject haben. Geld wird es ziemlich
sicher nicht geben, wie auch?

Ganz konkrete sollen als nächstes die ODE Solver optimiert werden. Und
das nächste grössere Aufgabengebiet wäre nach meiner Meinung die
Nullstellensuche.

Die Programmierung selber ist wahrscheinlich nichts mehr für Anfänger,
ich denke potentielle Mitprogrammierer sollten ein grundlegendes
Verständnis der gängigen C++ Techniken besitzen und sich mit Numerik und
Mathematik ein bisschen auskennen.

Also, falls jemand Lust und Laune hat sich hier mit einzuklinken, kann
er sich hier melden oder direkt bei SF. Wir freuen und über jeden
Mitstreiter.

Viele Grüße

[Walli]

Was spricht dagegen einfach ein schickes Interface auf stabile Implementierungen zu legen? Man muss solchen Standardquatsch ja nicht immer neu einhacken. Und gerade fuer Matrizen gibt es genug da draussen. Beispiele in annehmbarem C++ waeren MTL4 und boost::ublas.

[rüdiger]

Walli schrieb:

Was spricht dagegen einfach ein schickes Interface auf stabile Implementierungen zu legen? Man muss solchen Standardquatsch ja nicht immer neu einhacken. Und gerade fuer Matrizen gibt es genug da draussen. Beispiele in annehmbarem C++ waeren MTL4 und boost::ublas.

Den besten Eindruck macht imho Eigen.

[headmyshoulder]

Walli schrieb:

Was spricht dagegen einfach ein schickes Interface auf stabile Implementierungen zu legen? Man muss solchen Standardquatsch ja nicht immer neu einhacken. Und gerade fuer Matrizen gibt es genug da draussen. Beispiele in annehmbarem C++ waeren MTL4 und boost::ublas.

Nichts, und ich stimm Dir da voll zu. Wenn Du ein schickes Interface hast, sind 90% getan. Und den Rest dann einhacken kann man dann in einem Abwasch machen^^.

Btw. Lineare Algebra ist ein Fall für sich, der hiermit nicht komplett aufgerollt werden soll. Gerade für dünn besetzte Matrizen gibt es hochoptimierte Pakete, wie lapack, umfpack...

[headmyshoulder]

rüdiger schrieb:

Walli schrieb:

Was spricht dagegen einfach ein schickes Interface auf stabile Implementierungen zu legen? Man muss solchen Standardquatsch ja nicht immer neu einhacken. Und gerade fuer Matrizen gibt es genug da draussen. Beispiele in annehmbarem C++ waeren MTL4 und boost::ublas.

Den besten Eindruck macht imho Eigen.

Sieht interessant aus. Irgendeine Matrixabstrahierung brauchen wir auch und wir wollen da schon irgendetwas fertiges nehmen.

[stsa]

Hallo!

Versteht das Folgende nicht falsch, ich will euer Projekt nicht kaputtreden, aber einige Dinge sind mir doch unklar:

a) Wofür genau sollte ich obige Lib wirklich brauchen? Natürlich nehmt ihr mir Arbeit ab, indem ich mir über Datentypen keine Gedanken machen muss, aber mir stellt sich die Frage, zu welchem Preis. Ich bezweifle einfach, dass ihr eine Lib in der Qualität der fftw schreiben könnt/werdet. Von den Algorithmen aus NumRec zu einer solchen Lib ist es ein sehr sehr weiter Weg.

b) Brauche ich in der heutigen Zeit überhaupt noch eine solche CPU-orientierte Lib? Prinzipiell erschlage ich mittlerweile rechenintensive Probleme einfach mit Cuda/OpenCL, was auch auf Consumer Hardware recht günstig zu realisieren ist. Für alles andere im industriellen Umfeld stünden immernoch Quadro/Tesla Systeme zur Verfügung.

Ich kann mir vorstellen, dass ein solches Projekt unglaublich viel Spaß macht. Man lernt sicher viel dabei und die obige Projektbeschreibung ist in dem Forum ein wirkliches Highlight. Also nicht von meinen Zweifeln entmutigen lassen. Ich werde den Verlauf eures Projektes mit Interesse verfolgen.

[headmyshoulder]

stsa schrieb:

a) Wofür genau sollte ich obige Lib wirklich brauchen? Natürlich nehmt ihr mir Arbeit ab, indem ich mir über Datentypen keine Gedanken machen muss, aber mir stellt sich die Frage, zu welchem Preis. Ich bezweifle einfach, dass ihr eine Lib in der Qualität der fftw schreiben könnt/werdet. Von den Algorithmen aus NumRec zu einer solchen Lib ist es ein sehr sehr weiter Weg.

Stell Dir vor Du bist Doktorand/Diplomand und sollst numerische Probleme in C/C++ lösen, wo verschiedene Algorithmen und Methoden zusammen spielen. Ich bin beispielsweise Physiker und hab viel mit Differentialgleichungen (DGL)s zu tun. Ein Standardproblem ist die Fortsetzung von Lösungen der linearisierten DGLs ins Nichtlineare. Das ist prinzipiell eine Nullstellensuche. Mann kann dann von diesen Lösungen die Stabilität bestimmen, was quasi eine Eigenwertbestimmung ist und diese Lösung dann numerisch in der nichtlinearen DGL integrieren. Das mit NumRec oder der gsl zu machen ist sehr unschön, ständig gibt es irgendwelche nonstandards um auf die Arrays und Strukturen zuzugreifen. Und genau da wollen wir ansetzen und ein Lib schaffen die das mit den Standards macht. Es gibt in C++ schon einen Typ für komplexe Zahlen und zig Container für alles mögliche. Warum muss dann jede Bibliothek das alles nochmal neu schaffen?

In Python gibt es numpy und scipy, die machen genau sowas. Damit sind solche Problem sehr einfach zu lösen, aber in C++ meines Wissens nach nicht.

Die Qualität der Numerical Reciepes kann man locker erreichen. Die Solver, die wir für DGLs haben sind ziemlich genauso stabil und schnell wie NR, mal davon abgesehen haben wir auch Solver mit höheren Ordnungen und symplektische Integratoren. Mit der fftw siehts anders aus, das sind aber auch keine Standardmethoden mehr und sich mit damit zu messen würd ich mir nicht trauen^^.

stsa schrieb:

b) Brauche ich in der heutigen Zeit überhaupt noch eine solche CPU-orientierte Lib? Prinzipiell erschlage ich mittlerweile rechenintensive Probleme einfach mit Cuda/OpenCL, was auch auf Consumer Hardware recht günstig zu realisieren ist. Für alles andere im industriellen Umfeld stünden immernoch Quadro/Tesla Systeme zur Verfügung.

Nunja, wenn man Wert auf Geschwindigkeit legt, wird man natürlich etwas spezielles nehmen. Wir wollen hier auch keinen parallelisierten Code zum Lösen der Navier-Stokes Gleichungen oder zur Berechnung von Molekülorbitalen schreiben.

stsa schrieb:

Ich kann mir vorstellen, dass ein solches Projekt unglaublich viel Spaß macht. Man lernt sicher viel dabei und die obige Projektbeschreibung ist in dem Forum ein wirkliches Highlight. Also nicht von meinen Zweifeln entmutigen lassen. Ich werde den Verlauf eures Projektes mit Interesse verfolgen.

Jo, ich glaub auch, daß das sehr viel Spass machen wird und entmutigen lassen wir uns natürlich nicht, Kritik ist ja auch meistens hilfreich.

[hustbaer]

Hab mir interessehalber den Code etwas angesehen, dabei ist mir was aufgefallen.

// http://yanl.svn.sourceforge.net/viewvc/yanl/trunk/yanl/odeint/solver_runge_kutta4.hpp?revision=19&view=markup :
         std::vector<T> dxdt , dxt , dxm , xt ;

         void resize( size_t new_size ) {
             dxdt.resize( new_size , T(0.0) );
             dxt.resize( new_size , T(0.0)  ); // angenommen hier
             dxm.resize( new_size , T(0.0)  ); // oder hier
             xt.resize( new_size , T(0.0)  ); // oder hier fliegt eine exception
         }

         void next_step( CallMethod deriv, T t , T dt , T* state , size_t n) {
             if( n != dxdt.size() ) resize( n );
             // dann gibt's hier ein problem.
             //    -> lieber xt.size() abfragen

             T dh = dt * T( 0.5 ) , d6 = dt / T(6.0);
             T th = dt * T( 0.5 ) + t;

             deriv( state , dxdt.data() , t );
             for( size_t i=0 ; i<n ; i++ ) xt[i] = state[i] + dh * dxdt[i];

             deriv( xt.data() , dxt.data() , th );
             for( size_t i=0 ; i<n ; i++ ) xt[i] = state[i] + dh * dxt[i];

             deriv( xt.data() , dxm.data() , th );
             for( size_t i=0 ; i<n ; i++ ) { xt[i] = state[i] + dt * dxm[i] ; dxm[i] += dxt[i];}

             deriv( xt.data() , dxt.data() , t+dt );
             for( size_t i=0 ; i<n ; i++ ) state[i] += d6 * ( dxdt[i] + dxt[i] + T(2.0) * dxm[i] );
         }

BTW: wieso ist der code so unübersichtlich formatiert?

[Walli]

rüdiger schrieb:

Walli schrieb:

Was spricht dagegen einfach ein schickes Interface auf stabile Implementierungen zu legen? Man muss solchen Standardquatsch ja nicht immer neu einhacken. Und gerade fuer Matrizen gibt es genug da draussen. Beispiele in annehmbarem C++ waeren MTL4 und boost::ublas.

Den besten Eindruck macht imho Eigen.

Hast Du die mal ausprobiert? Besonders die Sparse-Implementierung interessiert mich. Ich mache manchmal serielle Testprogrämmchen und benutze wegen der Performance immer einen selbstgefrickelten PETSc-Wrapper weil uBLAS grottenlangsam ist und MTL4 noch nicht wirklich stable.

@OP: Wenn Du wirklich viele Leute (inkl. mich) glücklich machen willst, dann schreibe einen gescheiten C++-Wrapper für PETSc oder eine saubere C++-Matrix-Lib, die im Prinzip mit geringem Aufwand alles mögliche interfacen kann . Ich wollte das immer schon mal ernsthaft angehen, aber irgendwie finde ich die Zeit nicht.

http://www.mcs.anl.gov/petsc/petsc-as/

[headmyshoulder]

hustbaer schrieb:

Hab mir interessehalber den Code etwas angesehen, dabei ist mir was aufgefallen.

// http://yanl.svn.sourceforge.net/viewvc/yanl/trunk/yanl/odeint/solver_runge_kutta4.hpp?revision=19&view=markup :
         std::vector<T> dxdt , dxt , dxm , xt ;

         void resize( size_t new_size ) {
             dxdt.resize( new_size , T(0.0) );
             dxt.resize( new_size , T(0.0)  ); // angenommen hier
             dxm.resize( new_size , T(0.0)  ); // oder hier
             xt.resize( new_size , T(0.0)  ); // oder hier fliegt eine exception
         }

         void next_step( CallMethod deriv, T t , T dt , T* state , size_t n) {
             if( n != dxdt.size() ) resize( n );
             // dann gibt's hier ein problem.
             //    -> lieber xt.size() abfragen
             ...
        }

Was genau meinst Du? Falls da ein exception fliegt crasht alles, solange man von aussen nichts fängt, oder irre ich mich?

hustbaer schrieb:

BTW: wieso ist der code so unübersichtlich formatiert?

Über Codeformatierung haben wir uns noch keine Gedanken gemacht, aber so schlecht sieht das nicht aus, oder?

[headmyshoulder]

Walli schrieb:

Hast Du die mal ausprobiert? Besonders die Sparse-Implementierung interessiert mich. Ich mache manchmal serielle Testprogrämmchen und benutze wegen der Performance immer einen selbstgefrickelten PETSc-Wrapper weil uBLAS grottenlangsam ist und MTL4 noch nicht wirklich stable.

@OP: Wenn Du wirklich viele Leute (inkl. mich) glücklich machen willst, dann schreibe einen gescheiten C++-Wrapper für PETSc oder eine saubere C++-Matrix-Lib, die im Prinzip mit geringem Aufwand alles mögliche interfacen kann . Ich wollte das immer schon mal ernsthaft angehen, aber irgendwie finde ich die Zeit nicht.

http://www.mcs.anl.gov/petsc/petsc-as/

Die Hardcore Lineare Algebra wollen wir nicht unbedingt bearbeiten. Das ist echt sehr advanced und man braucht ne Menge Erfahrung, was Compiler können und muss die Mathematik 100%ig verstehen. Obwohl es natürlich sehr interessant ist.

[rüdiger]

@Walli
ich hab Eigen bisher nur für Kleinigkeiten benutzt. Daher kann ich nicht wirklich sagen, ob es alles hält was es verspricht. Aber bisher finde ich es ganz gut.

@headmyshoulder
Willst du das Projekt eher zum lernen machen oder für eine wirkliche Aufgabe?

das solltet ihr euch auf jeden Fall anschauen
http://www.nr.com/

Die ältere 2. Edition gibt es auch kostenlos: http://www.nrbook.com/a/bookcpdf.php

[hustbaer]

headmyshoulder schrieb:

hustbaer schrieb:

Hab mir interessehalber den Code etwas angesehen, dabei ist mir was aufgefallen.

// http://yanl.svn.sourceforge.net/viewvc/yanl/trunk/yanl/odeint/solver_runge_kutta4.hpp?revision=19&view=markup :
         std::vector<T> dxdt , dxt , dxm , xt ;
 
         void resize( size_t new_size ) {
             dxdt.resize( new_size , T(0.0) );
             dxt.resize( new_size , T(0.0)  ); // angenommen hier
             dxm.resize( new_size , T(0.0)  ); // oder hier
             xt.resize( new_size , T(0.0)  ); // oder hier fliegt eine exception
         }

         void next_step( CallMethod deriv, T t , T dt , T* state , size_t n) {
             if( n != dxdt.size() ) resize( n );
             // dann gibt's hier ein problem.
             //    -> lieber xt.size() abfragen
             ...
        }

Was genau meinst Du? Falls da ein exception fliegt crasht alles, solange man von aussen nichts fängt, oder irre ich mich?

Klar crasht alles wenn eine Exception fliegt, die keiner fängt.
Aber es könnte sie ja auch jemand fangen.

Ich meine folgende Situation:

1. next_step() wird aufgerufen, und z.B. dxt.resize() schmeisst eine Exception.
2. die Exception wird gefangen
3. next_step() wird nochmal aufgerufen - mit dem gleichen Wert für "n" wie beim ersten Aufruf
   next_step() vergleicht dann die Grösse von dxdt mit n, und sieht dass sie gleich sind. Und ruft daher resize() nicht nochmals auf.

Wenn "n" nun grösser ist, als die aktuelle Grösse von dxt, dmx oder xt, dann crasht der next_step() wahrscheinlich mit einer Access-Violation statt einer hübschen Exception. Oder noch schlimmer: liefert falsche Ergebnisse.

BTW: die "Lösung" die ich vorgeschlagen habe, ist auch nicht wirklich ausreichend.
Noch besser so in der Art:

template< class T, class CallMethod=dynamical_system<T>& >
     class ode_solver_runge_kutta4 : public ode_solver_base<T, CallMethod>
     {
         std::vector<T> m_dxdt;
         std::vector<T> m_dxt;
         std::vector<T> m_dxm;
         std::vector<T> m_xt;
         size_t m_buffer_size;

         void resize(size_t new_size)
         {
             if (new_size != m_buffer_size)
             {
                 m_buffer_size = 0; // if we're interrupted, force a re-start next time we're called
                 m_dxdt.resize(new_size, T(0.0));
                 m_dxt.resize(new_size, T(0.0));
                 m_dxm.resize(new_size, T(0.0));
                 m_xt.resize(new_size, T(0.0));
                 m_buffer_size = new_size;
             }
         }

     public:

         ode_solver_runge_kutta4() :
             m_buffer_size(0)
         {
         }

         void next_step(CallMethod deriv, T t , T dt , T* state , size_t n)
         {
             resize(n);

             T dh = dt * T(0.5)
             T d6 = dt / T(6.0);
             T th = dt * T(0.5) + t;

             deriv(state , m_dxdt.data() , t);
             for (size_t i = 0 ; i < n ; i++)
                 m_xt[i] = state[i] + dh * m_dxdt[i];

             deriv(m_xt.data() ,m_dxt.data(), th);
             for (size_t i = 0 ; i < n ; i++)
                 m_xt[i] = state[i] + dh * m_dxt[i];

             deriv(m_xt.data(), m_dxm.data(), th);
             for (size_t i = 0 ; i < n ; i++)
             {
                 m_xt[i] = state[i] + dt * m_dxm[i];
                 m_dxm[i] += m_dxt[i];
             }

             deriv(m_xt.data(), m_dxt.data(), t + dt);
             for (size_t i = 0; i < n ; i++)
                 state[i] += d6 * (m_dxdt[i] + m_dxt[i] + T(2.0) * m_dxm[i]);
         }
     };

[headmyshoulder]

hustbaer schrieb:

Ich meine folgende Situation:

1. next_step() wird aufgerufen, und z.B. dxt.resize() schmeisst eine Exception.
2. die Exception wird gefangen
3. next_step() wird nochmal aufgerufen - mit dem gleichen Wert für "n" wie beim ersten Aufruf
   next_step() vergleicht dann die Grösse von dxdt mit n, und sieht dass sie gleich sind. Und ruft daher resize() nicht nochmals auf.

Ja, du hast Recht, das muß man beachten und ändern. Ich glaub dein Vorschlag ist auch eleganter als nach der Größe von xt zu fragen. Das wird geändert.

Ich glaub ich weiss jetzt auch was Du mit unschöner Codeformatierung meinst, obwohl ich members mit m_ zu bezeichnen nicht so mag.

[headmyshoulder]

rüdiger schrieb:

@headmyshoulder
Willst du das Projekt eher zum lernen machen oder für eine wirkliche Aufgabe?

Ich sehe dieses Projekt als richtige Aufgabe, wo was vernünftiges und benutzbares rauskommen soll. Die Idee ist ja nicht aus Spass an der Freude gekommen, sondern weil wir solch eine Bibliothek oft benutzen werden und uns die existierenden Lösungen nicht sehr gut gefallen.

Aber natürlich wird man dabei auch eine Menge lernen, sowohl thematisch als auch programmiertechnisch.

rüdiger schrieb:

das solltet ihr euch auf jeden Fall anschauen
http://www.nr.com/

Die ältere 2. Edition gibt es auch kostenlos: http://www.nrbook.com/a/bookcpdf.php

Ja, die Numerical Reciepes sind das Standardbuch für solche Probleme. Ich kenn das fast auswendig^^, zumindest einige Kapitel.