exp() Funktion selber implementieren



  • Hallo zusammen,

    als Übung soll ich die exp() Funktion selber implementieren.
    Ich sitze schon seit einer gefühlten Ewigkeit an der Aufgabe und weiß absolut nicht wo mein Fehler ist. Das hier ist meine Implementierung:

    long double ehoch(long double x){
        long double erg = 0;
        long long counter = 0;
        long double epsilon = 0.00001;
        long double next = pow(x,counter) / fakultaet(counter);
        counter++;
        while(epsilon < abs(next) && counter < 100){
            erg += next;
            next = pow(x,counter) / fakultaet(counter);
            counter++;
        }
        return erg;
    }
    

    Ich versuche hier mit der Reihenentwicklung der Exponentialfunktion https://de.wikipedia.org/wiki/Exponentialfunktion auf eine bestimmte Genauigkeit (epsilon) den Wert zu errechnen. Eigentlich soll next für counter gegen unendlich gegen 0 gehen. Sprich für ein großes counter immer kleiner werden und schließlich kleiner epsilon abs(next) sein. Allerdings wird next immer größer und ich weiß nicht wieso...

    Kann mir da jemand weiterhelfen?

    MfG
    Asse

    Edit: fakultaet habe ich auch selber implementiert das funktioniert meines Erachtens aber einwandfrei:

    Aus Verzweiflung arbeite ich überall schon mit den größten mir bekannten Datentypen.

    long long fakultaet(long long x){
        if(x < 0){
            //Fehler
            return -1;
        }
        if(x == 0 || x == 1){
            return 1;
        }
        long long erg = 2;
        for(long long i = 3; i <= x; i++)
        {
            erg *= i;
        }
        return erg;    
    }
    


  • @Asse sagte in exp() Funktion selber implementieren:

    Edit: fakultaet habe ich auch selber implementiert das funktioniert meines Erachtens aber einwandfrei:

    Die brauchst du gar nicht.

    Worin unterscheiden sich die Summanden der Reihe?
    Nur diesen Faktor musst du neu berechnen.



  • Zumindestens für kleine x funktioniert deine Funktion (annähernd): Ideone-Code

    Aber so wie @DirkB schrieb, brauchst du nicht immer wieder den kompletten Teilausdruck neu berechnen.
    Worin unterscheiden sich jeweils Zähler und Nenner (vom Vorgänger)?


  • Mod

    Und für große x läufst du in numerische Probleme mit der Fakultät. Merke: Reihenentwicklungen, bei denen Fakultäten von immer größer werdenden Zahlen benötigt werden, sind für numerische Berechnungen ungeeignet. Speziell kommen in der Taylerreihe nicht nur jede Menge Fakultäten vor, die sich meistens nicht wegkürzen, sie konvergiert auch noch denkbar langsam. Lass die Taylorreihe in den theoretischen Beweisen (dafür ist sie super!) und nimm für Numerik etwas anderes. Oder vereinfache die Taylorentwicklung so, dass die Fakultäten wegfallen. Siehe dazu das, was @DirkB geschrieben hat.



  • Ok, stimmt. Ich brauche nicht immer alles neu zu berechnen.
    Aber wenn ich das ändern würde müsste das sich aber dennoch numerisch auslöschen, oder?

    Ich habe im Internet geschaut, wie man das numerisch umsetzt kann https://de.wikipedia.org/wiki/Exponentialfunktion#Numerische_Berechnungsmöglichkeiten
    Allerdings durchblicke ich nicht so ganz was da geschildert ist und wie man das dann umsetzt...



  • Das brauchst du hier nicht.

    Deine Schleife mußt du nur minimal ändern. Wie lautet denn der sich verändernde Faktor je Summand?

    PS: Durch passende Änderung des (Ideone-)Codes kommt dann folgende Ausgabe:

    0: 1 <-> 1
    1: 2.71827877 <-> 2.718281828
    2: 7.389046016 <-> 7.389056099
    3: 20.08553443 <-> 20.08553692
    4: 54.59814722 <-> 54.59815003
    5: 148.4131471 <-> 148.4131591
    6: 403.4287835 <-> 403.4287935
    7: 1096.63315 <-> 1096.633158
    8: 2980.957981 <-> 2980.957987
    9: 8103.083923 <-> 8103.083928
    


  • Ich habe mal umgesetzt, dass es nicht immer neu berechnet wird. Das hat super geholfen:

    long double ehoch(double x){
        if(x == 0)
            return 1.;
        bool invers = false;
        if(x < 0){
            invers = true;
            x *= -1;
        }
        long double erg = 0;
        long long counter = 1;
        double epsilon = 0.00001;
        long double numerator = 1;
        long double denominator = 1;
        while(epsilon < (numerator/denominator)){
            erg += numerator/denominator;
            numerator *=  x;
            denominator *= counter;
            counter++;
        }
        return invers ? 1/erg : erg;
    }
    

    Für negative Exponenten lasse ich einfach Exp(|x|) berechnen und gebe dann die Inverse zurück. Also 1/exp(|x|).
    Jetzt funktioniert es mit der gewünschten Genauigkeit.

    Edit: Ich nutze jetzt noch long double. Dann funktioniert das auch für große Exponenten.



  • Ich frage mich aber weiterhin wieso das komplette Neuberechnen von pow() und fakultaet() zur nummerischen Auslöschung führt und diese Lösung nicht. Letztlich werden bei der Neuberechnung doch einfach alle Schritte "wiederholt" und dann das nächste multipliziert und zurückgegeben. In meinem jetzigen Fall spare ich mich doch nur die Wiederholungen und multipliziere direkt das neue dran. Wo ist denn nun der Trick versteckt?


  • Mod

    Mit dem "Trick" kommt keines deiner Zwischenergebnisse jemals in ungünstige Wertebereiche. Wenn du pow durch die Fakultät teilst, dann sind zwischendurch der Rückgabewert von pow und der Rückgabewert der Fakultät beide sehr groß. Das ist numerisch sehr ungünstig. Auch wenn das Ergebnis am Ende klein sein mag, ist es dann unnötig ungenau. Bei der Fakultät sprengst du sogar ganz schnell den Wertebereich sogar von long long und bekommst dann ein komplett falsches Ergebnis.

    Als Faustregel willst du immer mit möglichst gleich großen und allgemein mit nicht zu großen Zahlen arbeiten (außer bei Subtraktion).

    Damit ergibt sich auch direkt ein weiterer Tipp: Du addierst hier lauter kleine Zahlen zu einer immer größer werdenden Summe. Das ist auch ungünstig, weil diese Addition dadurch mit jedem Schritt numerisch ungenauer wird. Idealerweise möchte man mit dem Addieren der kleinen Zahlen anfangen und dann zu den größeren Zahlen kommen, so dass die laufende Summe und der neue Faktor stets ungefähr gleich groß sind.

    Das ist hier leider nicht so ganz einfach, aber als einfache und effektive Maßnahme betrachte mal den ersten Faktor: Der ist immer 1. Das ist ungünstig, denn 1 plus etwas kleines ist numerisch viel ungenauer als 0 plus etwas kleines. Daher: Lass den ersten Faktor weg und berechne erst exp(x) -1. Beim return am Ende addierst du dann die fehlende 1. So hast du gerade bei kleinen Exponenten durch einfaches Umstellen der Rechnung deutlich an Genauigkeit gewonnen.



  • Hier noch meine Lösung: Ideone-Code
    Bei Reduzierung von epsilon kommen auf 10 Stellen exakt die Werte von exp(x) heraus:

    0: 1 <-> 1
    1: 2.718281828 <-> 2.718281828
    2: 7.389056099 <-> 7.389056099
    3: 20.08553692 <-> 20.08553692
    4: 54.59815003 <-> 54.59815003
    5: 148.4131591 <-> 148.4131591
    6: 403.4287935 <-> 403.4287935
    7: 1096.633158 <-> 1096.633158
    8: 2980.957987 <-> 2980.957987
    9: 8103.083928 <-> 8103.083928
    

  • Mod

    @Th69 sagte in exp() Funktion selber implementieren:

    Hier noch meine Lösung: Ideone-Code
    Bei Reduzierung von epsilon kommen auf 10 Stellen exakt die Werte von exp(x) heraus:

    0: 1 <-> 1
    1: 2.718281828 <-> 2.718281828
    2: 7.389056099 <-> 7.389056099
    3: 20.08553692 <-> 20.08553692
    4: 54.59815003 <-> 54.59815003
    5: 148.4131591 <-> 148.4131591
    6: 403.4287935 <-> 403.4287935
    7: 1096.633158 <-> 1096.633158
    8: 2980.957987 <-> 2980.957987
    9: 8103.083928 <-> 8103.083928
    

    Und mit meiner Verbesserung bekommen wir fast noch eine Stelle mehr geschenkt, ohne das epsilon zu verkleinern zu müssen. Einzig exp(1) kommt durch eine ungünstige Rundung zu einem anderen Ergebnis:
    https://ideone.com/ZRvXHM

    0: 1 <-> 1
    1: 2.7182818284 <-> 2.7182818285
    2: 7.3890560989 <-> 7.3890560989
    3: 20.085536923 <-> 20.085536923
    4: 54.598150033 <-> 54.598150033
    5: 148.4131591 <-> 148.4131591
    6: 403.42879349 <-> 403.42879349
    7: 1096.6331584 <-> 1096.6331584
    8: 2980.957987 <-> 2980.957987
    9: 8103.0839276 <-> 8103.0839276
    


  • @Asse sagte in exp() Funktion selber implementieren:

    Edit: fakultaet habe ich auch selber implementiert das funktioniert meines Erachtens aber einwandfrei:
    Aus Verzweiflung arbeite ich überall schon mit den größten mir bekannten Datentypen.

    long long fakultaet(long long x){  // long long für x ist ganz sicher übertrieben (char würde reicht)
        if(x < 0){                     // ok ein unsigned
            //Fehler
            return -1;
        }
     ....
    
    

    Das Problem liegt auch daran, dass schon ab 22! der Wertebereich vom long long (+/- 63-Bit) überschritten wird.
    unsigned long long macht es auch nicht besser.



  • @SeppJ: Ich habe mal beide Varianten auf 16 Stellen erweitert (sowie epsilon auf 1e-16 gesetzt):

        0: 1 <-> 1
        1: 2.718281828459045 <-> 2.718281828459045
        2: 7.38905609893065 <-> 7.38905609893065
        3: 20.08553692318767 <-> 20.08553692318767
        4: 54.59815003314424 <-> 54.59815003314424
        5: 148.4131591025766 <-> 148.4131591025766
        6: 403.4287934927351 <-> 403.4287934927351
        7: 1096.633158428459 <-> 1096.633158428459
        8: 2980.957987041728 <-> 2980.957987041728
        9: 8103.083927575384 <-> 8103.083927575384
    

    Auch ohne das Ändern von epsilon waren die beiden Ausgaben schon identisch (nur eben noch in den letzten 4-5 Stellen abweichend von exp(x)).



  • Ich danke euch für die Hilfe!!

    Als finale Version nehme ich folgenden Code.

    long double ehoch(double x){
        if(x == 0)
            return 1.;
        bool invers = false;
        if(x < 0){
            invers = true;
            x *= -1;
        }
        long double erg = 0;
        long long counter = 2;
        long double epsilon = 1e-16;
        long double next = x;
        while(epsilon < abs(next)){
            erg += next;
            next *= x / counter;
            counter++;
        }
        return invers ? 1/(erg + 1) : (erg + 1);
    }
    

    Der Funktioniert dann auch für negative Exponenten und ist auch für große Exponenten recht genau.



  • @SeppJ Ok das verstehe ich. Super Tipps, danke!

    Also wäre es nummerisch besser, wenn man erst die Summanden "von klein nach groß" addiert?
    Die Reihe, die wir hier verwenden liefert für die ersten paar Indizes noch relativ kleiner Werte, wächst dann stark an und wird ab einen gewissen Index immer kleiner und geht gegen Null. Müsste es dann helfen die ganze Rechnung "rückwärts" zu machen? Z.B in einer schleife den Summanden next aufzubauen ohne an erg zu addieren und dann ab dem Eintauchindex wo expsilon > abs(next) ist rückwärts anzufangen auf erg zu addieren und dann statt next *= x / counter -> next *= i / x wobei i Index der Schleife vom Eintauchindex bis 0 ist?!

    Ich hoffe die Idee ist verständlich.

    Edit: Eine andere Idee wäre es die einzelnen Summanden in einem Array oder Vector zu speichern, dann der größe nach zu sortieren und dann nacheinander addieren?!


  • Mod

    Hier ist es nicht so ganz wichtig, darauf zu achten, weil mit dem epsilon eine gewisse zu erreichende Genauigkeit vorgegeben ist und genau so lange gerechnet wird, bis diese erreicht ist. Ob die Stellen danach mit gleichem Aufwand auch noch hätten besser sein können, ist dann weniger interessant. Meistens hat man es eher umgekehrt, dass man die Zahl der Reihenglieder vorgibt, und dann eben hofft, dass das Ergebnis so genau wie möglich ist. Dann ist es wichtiger, dass man mit dem vorgegebenen Rechenaufwand ien möglichst gutes Ergebnis zu erzielen.



  • @SeppJ Ok, danke für deine Hilfe!!

    Mit meinem oben genannten Ergebnis bin ich auf jeden Fall zufrieden.


  • Banned

    Falls das hier jemanden interessiert: So haben das die alten Taschenrechner gemacht:
    https://de.wikipedia.org/wiki/CORDIC


Log in to reply