3. Problem mit Gleitkommazahlen

Problem mit Gleitkommazahlen

  • ?

    Ich habe ein Problem, welches mich doch etwas zum verzweifeln bringt.
    Im folgenden ist ein Stück meines Codes angegeben. Die Funktion "Intgr_log"
    wird jedesmal garantiert mit den gleichen Argumenten aufgerufen und liefert
    auch bis zur 19ten Nachkommastelle den gleichen Wert zurück:

    double fff= Intgr_log(&MeanExcitationEn, E1)-Intgr_log(&MeanExcitationEn, E1),
         d1= Intgr_log(&MeanExcitationEn, E1),
         d2= Intgr_log(&MeanExcitationEn, E1),
         d3= d2-d1;

  cout<<"  "<<d2<<endl<<" -"<<d1<<endl<<" ="<<d3<<endl<<"!="<<fff<<endl;

  printf("fff: %Lx \r\n",fff);
  printf(" d1: %Lx \r\n",d1);
  printf(" d2: %Lx \r\n",d2);

    Das merkwürdige Ergebnis ist:

    0.9223072223880005582
    -0.9223072223880005582
    =0
    !=-4.163336342344337026e-17

    fff: bc88000000000000
    d1: 3fed838a6fa0ac89
    d2: 3fed838a6fa0ac89

    Wie kann das sein, dass unterschiedliche Ergebnisse herauskommen, wenn ich die
    zurückgegebenen Werte zwischenspeichere und subtrahiere (sinnvoll: =0) oder direkt subtrahiere (sinnlos: !=0)!
    Die Rückgabewerte sind bitweise gleich.
    Hat jemand eine Idee was dort passieren könnte?
    Danke

    0
  • M

    der computer erkennt ja bei deiner subtraktion nicht das die beiden gleich sind. ich denke des er den einen wert vieleicht negiert und es dann addiert, dann könnte sich der fehler bei dem negieren einschleichen. (nur spekulation)

    mfg Mirauder Mo

    0
  • R

    direkt seh ich keinen Fehler

    Aber was liefert Intgr_log zurück? Vielleicht ist das ja kein double sondern nur ein Typ, der implizit nach double konvertiert werden kann und der - operator kommt nicht mit dem Typ klar.

    und bist du dir sicher, dass du Intgr_log immer mit den gleichen Wertan aufrufst? Weil der erste Parameter wird ja als Pointer übergeben und könnte ja intern verändert werden.

    0
  • ?

    Beim negieren wird beim double-Typ ja eigentlich nur das Vorzeichenbit
    invertiert. Die Mantisse der Zahlendarstellung bleibt aber ja die gleiche,
    somit sollte auch das Ergebnis 0 sein ...

    Grüße Uwe

    0
  • ?

    Die "Intgr_log"-funktion ist durch Vererbung ein Mitgliedsfunktion der Klasse aus der der Code stammt geworden. Ihr Rückgabetyp ist wirklich double. Alle Befehle in dieser Funktion verwenden den Zeiger nur als Argument. Er wird nicht für Schreibzugriffe benutzt.

    0
  • ?

    Gerade hab ich noch etwas mehr heraus bekommen.
    Verwende ich die Funktion so wie hier:

    double CInterpol::Intgr_log( vector<S_Wertepaar>* d,
                             double x) {

...

   return I_log((*d)[idx].x,   (*d)[idx].y,(*d)[idx+1].x, (*d)[idx+1].y, x)
          + Offset;
}

    ist das Ergebnis ungleich Null.

    Schreibe ich die Funktion um zu:

    double CInterpol::Intgr_log( vector<S_Wertepaar>* d,
                             double x) {

...

   double ddd = I_log((*d)[idx].x, (*d)[idx].y, (*d)[idx+1].x, (*d)[idx+1].y,x)
               +Offset;
   return ddd;

}

    erhalte ich Null als Ergebnis für 'fff'.
    Weiß vielleicht jemand wieso oder was daran falsch ist?

    Danke allen Bisherigen

    Grüße Uwe

    0
  • P

    Daran ist erst erstmal gar nichts falsch.
    Du arbeitest hier mit Zahlen doppelter Genuaigkeit d.h nach IEEE 64 Bit damit gelten folgende
    Grenzen entnommen float.h MSVC 6.0

    #define DBL_DIG         15                      /* # of decimal digits of precision */
#define DBL_EPSILON     2.2204460492503131e-016 /* smallest such that 1.0+DBL_EPSILON != 1.0 */
#define DBL_MANT_DIG    53                      /* # of bits in mantissa */
#define DBL_MAX         1.7976931348623158e+308 /* max value */
#define DBL_MAX_10_EXP  308                     /* max decimal exponent */
#define DBL_MAX_EXP     1024                    /* max binary exponent */
#define DBL_MIN         2.2250738585072014e-308 /* min positive value */
#define DBL_MIN_10_EXP  (-307)                  /* min decimal exponent */
#define DBL_MIN_EXP     (-1021)                 /* min binary exponent */
#define _DBL_RADIX      2                       /* exponent radix */
#define _DBL_ROUNDS     1                       /* addition rounding: near */

    Damit sind deine Zahlen nicht auf 19 Stellen gleich
    sondern nur auf 15 alles was darunter ist ist Zufall bzw Ungenauigkeiten bei der Formatierung.

    Schade das du nicht auch d3 als Hexwert angegebn hasst.

    @ ur1: Ich bin mir nicht sicher ob wirklich nur das Sign Bit der Mantisse sich ändert
    bei int ist es ja deutlich anders
    +1 == 0x0001
    -1 == 0xFFFF
    Aber ich habe im Moment keine Möglichkeit das nachzuschauen.

    Floats werden üblicherweise nicht mit == verglichen
    sondern es wird überprüft ob die Differenz der Beiden Zahlen innerhalb eines vorgegebenen Intgervalls ist hier bietet sich typischerweise das Interval ]-2*DBL_MIN..+2*DBL_MIN[ an.
    kann man auch ausdrücken der Betrag der Differenz soll kleiner sein als 2*DBL_MIN.

    💡 Zur Problematik der Rechengenauigkeit mit Gleitkommazahlen solltest du vielleicht mal in Mathe forum gehn richtung numerische Mathematik oder nach den Stichworten googlen.

    0
  • P

    ich glaube ein auf DBL_EPSILON basierendendes intervall ist deutlich sinnvoller als ein auf DBL_MIN basierendes.
    weil die fehler treten ja ca. in der genauigkeit auf, in der man rechnet, und meist ist man eher um 1 als bei e-308.
    noch sinnvoller ist vielleicht eine relative differenz. also wenn die differenz kleiner als e-10*(summe der zu vergleichenden werte) wird gleichheit angenommen. dann ist es egal, in welcher groessenordnung man sich befindet.

    0
  • P

    @PeterTheMaster

    💡 Danke für die Korrektur hab die falsche der beiden erwischt.

    😕 Was ist der Unterschied zwischen relativer differenz und dem was ich beschrieben habe

    |Zahl_B - Zahl_A| < 2*DBL_EPSILON => die beiden Gleitkommazahlen sind gleich

    Gleitkomma deshalb weil das Verfahren für float, double, long double funktioniert, natürlich mit anderen Epsilons

    0
  • P

    na statt
    |Zahl_B - Zahl_A| < 2DBL_EPSILON
    lieber
    |z1-z2|<1e-10    (|z1|+|z2|)
    damit passt man sich automatisch an genauigkeit und groessenordnung der argumente an. man weiss ja nicht bevor man diese vergleichsfunktion schreibt, in welchen groessenordungen man sich mal rumtreibt.
    und bei 1e-200 - 1e-201 wird man nach der ersten variante niemals einen unterschied feststellen koennen.
    man koennte natuerlich bei der vergleichsfunktion auch eine konstante fuer die intervallgroesse mit angeben lassen. aber man will ja moeglichst wenig schreiben und das moeglichst automatisch haben.

    0
  • P

    Ich werd langsam alt, habe gerade was von mir gefunden, ist heute noch im Einsatz.

    /***********************************************************************.FA*
 .FUNCTION [ floatcomp                                    ]
--------------------------------------------------------------------------
 .GROUP    [ GenUtil                      ]
 .AUTHOR   [ PAD                          ]
--------------------------------------------------------------------------
 .DESCRIPTION
    comparing two float values for equality
    two float values are equal if they are in an intervall
    of [(-0.25e-12 * max(a,b))..(+0.25e-12(max(a,b))]

--------------------------------------------------------------------------
 .INDEX
   Utilities
--------------------------------------------------------------------------
 .PARAMETER
   IN   long double a  these two variables get compared
        long double b

--------------------------------------------------------------------------
 .RETURNVALUE
   0  equal within limits
   1  notequal within limits
--------------------------------------------------------------------------
 .VARIABLE_REFERENCES
--------------------------------------------------------------------------
 .HISTORY
  Date      Author          Comment
  09.05.89  PAD        comments adapted for SOFTDOC
  16.07.89  PAD        Empfindlichkeit von e-14 nach
                            e-12 verringert
**********************************************************************.HE**/

//#define FEQU(a,b) floatcomp((long double)(a),(long double)(b))
//#define FNEQU(a,b) (!floatcomp((long double)(a),(long double)(b)))
//#define TESTMAIN

int  floatcomp(double a, double b)
{
  double fuzz = 0.5e-12L;
  double c,d;
  c=fabs((double)(a-b));
  d=(fuzz * max(fabs((double) a),fabs((double)b)));
  if (c <= d)
    return PASS;
  else return FAIL;

}

#ifdef TESTMAIN
main()
{
  long double f = 0L;

  while (FEQU(1,f))
  {
    printf ("%.30Lf, %.30Lf\n",f,(long double)f);
    f+= 0.1L;
  }

  f=0.0L;
  while (FNEQU(0,f))
  {
    printf ("%.30Lf, %.30Lf\n",f,(long double)f);
    f+= 0.1L;
  }
  return 0;
}
#endif

    🙂

    0
Anmelden zum Antworten