/*
 * Dominik Erdmann
 * 2004-11-25
 *
 * Programm zur numerischen Loesung von Integralen
 * nach dem Romberg - Verfahren
 *
 * In der Methode start() muessen die Grenzen eingetragen werden
 * und eine beliebige Genauigkeit epsilon
 *
 * In der Methode f(double x) ist die zu berechnende Funktion einzugeben
 *
 */


#include <iostream>
#include <cmath>

using namespace std;

/* CLASS *******************************************************/

class integral {
    public:
        void start();
        ~integral();
    
    private:
        long double grenzeOben,
                    grenzeUnten,
                    schrittweite,
                    epsilon,
                    anz; // Anzahl der Schritte
        long double **erg; // fuer Array mit Pointern auf Array mit long double
        long double **temp; // Pointer, um alte Pointer zu sichern
        
        int genauigkeit; // in Nachkommastellen
        int hoechsterGrad; // der benoetigte Grad wird gespeichert
        // Endergebnis steht in erg[hoechster-1][hoechsterGrad-1]
        
        bool abbruch; // fuer die Abbruchbedingung
        
        long double pi(); // liefert PI zurueck, wenn fuer Grenzen gewuenscht
        void   romberg(); // der Algorithmus
        double f(double x); // zu berechnende Funktion
};

/* MAIN ********************************************************/

main() {
    char c;
    integral i;
    i.start();
    
    cout << "\n\n\nZum Beenden beliebiges Zeichen eingeben! ";
    cin >> c;
}

/* DESTRUKTOR ******************************************************/

integral::~integral() {
    for (int i = 0; i < hoechsterGrad; ++i)    
        delete [] erg[i]; // was man belegt, sollte man auch freigeben
        
    delete [] erg;
}

/* START() *****************************************************/

void integral::start() {

    grenzeUnten = 0;  // hier die Grenzen bestimmen ( PI = pi() )
    grenzeOben  = 10;
    genauigkeit = 5;  // in Nachkommastellen
    epsilon = pow(10.0, -1 * genauigkeit); // und dieFehlertoleranz
    // dran denken: Funktion ganz unten eingeben!
    
    cout.precision(16); // Nachkommastellen der Ausgabe
    
    romberg(); // los gehts
    
    cout << "\n\nErgebnis: " << erg[hoechsterGrad-1][hoechsterGrad-1]
         << "\nbei Grad von: " << hoechsterGrad;
}

/* PI() ********************************************************/

long double integral::pi() {
    return acos(-1.0);
}

/* ROMBERG() **********************************************/

void integral::romberg() {
    
    // Benutze ein Array mit Pointern auf ein Array mit double long
    erg = new long double* [1]; // neuen Pointer fuer Array
    erg[0] = new long double [1]; // neuer Speicher fuer Array - hier kein Array
    
    schrittweite = (grenzeOben - grenzeUnten);
    
    // bedenke: Indexanfang von Array bei 0!!
    // also: Index von Romberg ist erg[1][1]
    erg[0][0] = ( f(grenzeOben) + f(grenzeUnten) ) * schrittweite / 2;
    cout << "Romberg[1][1] = " << erg[0][0] << "\n\n";
    
    abbruch = false;
    for (int k = 2; abbruch == false ; ++k )
    {
        // Um das Array mit den Pointern auf ein Array mit long double
        // zu vergroessern, muss man die alten Pointer sichern, neue Pointer
        // anlegen und die Alten zurueckschreiben.
        
        // temporaeres Pointerarray allokieren
        //long double *temp[k-1];
        temp = new long double*[k-1];
        
        // alte Pointer sichern
        for (int i = 0; i < k - 1; ++i)
            temp[i] = erg[i];

        // Pointerarray loeschen und dann um 1 vergroessern
        delete [] erg;
        erg = new long double*[k];
        
        // alte Pointer zurueckschreiben
        for (int i = 0; i < k - 1; ++i)
            erg[i] = temp[i];

        // temporaere Pointer loeschen
        delete [] temp;

        // neues Array mit long double anlegen und zu den Pointern hinzufuegen
        erg[k-1] = new long double [k];
        
        
        // hier beginnt der eigentliche Algorithmus
        for (int j = 1; j <= k; ++j)
        {
            if ( j == 1 )
            {
                /* Schrittweite = Intervalllaenge / 2^{k-1} */
                anz = pow(2.0, k-1);
                schrittweite = (grenzeOben - grenzeUnten) / anz;
                
                // Initialisierung. Auch hier: Index von Romberg ist erg[k][1]
                erg[k-1][0] = 0.0;
                
                // die Summe der Stuetzstellen
                for (int l = 1; l <= anz/2; ++l)
                    erg[k-1][0] +=  f( grenzeUnten + (2*l-1) * schrittweite );
                
                erg[k-1][0] = erg[k-1][0] * schrittweite * 2;
                erg[k-1][0] += erg[k-2][0];
                erg[k-1][0] = erg[k-1][0] / 2;
            }
            else
            {   // Indexanfang von Array bei 0!!
                erg[k-1][j-1] = erg[k-1][j-2] + // eine Addition
                                ( erg[k-1][j-2] - erg[k-2][j-2] ) // Division
                                / ( pow(4.0, j-1) - 1 );
            }
            
            // Abbruchbedingung ueberpruefen
            // k==j==3 =>> immer hoechster Index der Zeile und min 3. Grad
            if ( k == j && k >= 3 
                 && abs( erg[k-1][k-1] - erg[k-2][k-2] ) < epsilon 
                 && abs( erg[k-3][k-3] - erg[k-2][k-2] ) < epsilon ) 
            {
                abbruch = true;
                hoechsterGrad = k;
            }
            
            // hier gehts zur Ausgabe
            if (k == j)
            {
                cout << "Romberg[" << k << "][" << j << "] = " 
                     << erg[k-1][j-1] << "\n";
            
                if (k >= 3)   // die Abweichung gibt es erst ab Grad 3
                    cout << "Abweichung: "
                         << abs(erg[k-1][j-1]-erg[k-2][j-2]) << "\n\n";
                else
                    cout << "\n";
            }

        } // zweite for-Schleife
    } // erste for-Schleife

}

/* F() *********************************************************/

double integral::f(double x) {
    // Hier die zu berechnende Funktion eingeben!
    return exp(x); // e^x
}