Mythos Top-Down-Schleife

In einem früheren Blog mutmasste ich, dass ich in test_calc_reduce.c die Macros dot_prod und square_sum bereits eliminiert hätte. Diese Mutmassung war falsch. Es passiert mir oft, dass ich in der Erinnerung das intensiv ins Auge gefasste Vorhaben, etwas zu tun, damit verwechsle, ich hätte es bereits getan - wohl ein allgemein-menschliches Problem, das etwas mit der Art zu tun hat, wie unser Gedächtnis funktioniert.[1]

Also führte ich es nun aus: Ich führte zwei neue Inline-Funktionen in test_calc_reduce.h ein:

static inline double vec_dot( void* v, void* w, int n) {
  double *fv = (double *)v,
         *fw = (double *)w,
         s = 0;
  if (v && w) for (int i=n-1;i>=0;--i) s += fv[i]*fw[i];
  return s;
  }

static inline double vec_length( void* v, int n) {
  return v ? sqrt( vec_dot(v,v,n) ) : 0;
  }

Dann ersetzte ich die Aufrufe der Macros dot_prod und square_sum durch die entsprechenden Inline-Funktionen. Während nach der dot_prod-Ersetzung die Tests noch sauber durchliefen, kam es nach der square_sum-Ersetzung (alles schön in einzelnen Versionen der Dropbox nachvollziehbar!) zu einer Verletzung eines Tests. Ein einziger Test schlug fehl - die Berechnung einer Geschwindigkeitskomponente des vermeintlich "wahren" Mondknotens (Swiss Ephemeris Objekt 13) war in einem einzigen Fall zu ungenau:

ruediger@herschel:~/workspace/swepar/src$ test_calc_reduce tcr.bin

          xx[3] :    2.5273368038 <>    2.5273404552
Differences from test case 138
Juldat 2263323.5683292975, planet 13, iflags 258

Swiss Eph, geocentric positions: 249 OK   1 NOT OK

Natürlich prüfte ich zuerst sehr genau, ob ich irgendwelche Fehler bei den Ersetzungen gemacht hatte. Sie waren aber alle korrekt. Als nächstes debuggte ich den Testfall und verglich die Ergebnisse von square_sum und vec_dot bei jedem Aufruf. Dabei kam es zu kleinen Abweichungen. Für die Werte

xpos[2][3]  double  -2.4692434310721423e-05
xpos[2][4]  double  -0.00055779853107697566
xpos[2][5]  double  -4.9608092113099317e-05

ergab die Berechnung

v2 = vec_dot(xpos[i]+3,xpos[i]+3,3);  // 3.1420988038692279e-07
v2 = square_sum((xpos[i]+3));         // 3.1420988038692284e-07 
// Exakter Wert mit Wolfram Alpha:    // 3.14209880386922820362148922449447018e-07

Ärgerlicherweise liefert das Macro sogar den genaueren Wert, obwohl man doch annehmen sollte, dass nach Auflösung der Inline-Funktion praktisch dieselben Maschineninstruktionen entstehen! Die exakte Berechnung führte ich mit den obigen Werten für xpos[2][3..5] mit der Suchmaschine Wolfram Alpha durch.

Wie kommt es zu diesem Unterschied?

Die Erklärung ist banal: Das Rechnen mit Gleitkommazahlen verletzt das Assoziativgesetz aus Rundungsgründen. Es kommt also darauf an, in welcher Reihenfolge man eine Summe mit drei Summanden berechnet. Dagegen brachte Wolfram Alpha, wie zu erwarten, immer dasselbe Ergebnis, unabhängig von der Reihenfolge der Summanden, da sie solche Berechnungen mit arbitrary precision ausführt.

Nun hatte ich aber das Skalarprodukt mit einer Top-Down-Schleife berechnet, wobei ich einem Performance-Mythos aufgesessen war, der möglicherweise in den 80er Jahren noch seine Berechtigung gehabt hatte.[2]

Nach Umstellung auf eine normale Bottom-Up-Schleife,

static inline double vec_dot( void* v, void* w, int n) {
  double *fv = (double *)v,
         *fw = (double *)w,
         s = 0;
  if (v && w) for (int i=0;i<n;i++) s += fv[i]*fw[i];
  return s;
  }

die auch für das Auge angenehmer = lesbarer ausfällt, war wieder alles OK, die Inline-Funktion ergab bis auf die letzte Dezimale denselben Wert wie das Macro.

---

[1] Das ist wohl ein allgemein-menschliches Problem, Ich habe den Eindruck, bei der Erinnerung spielen die Wertungen "wichtig"/"unwichtig" bzw. "interessant"/"uninteressant" eine Rolle, mit der Tendenz, dass Unwichtiges oder Uninteressantes aus dem Gedächtnis herausfällt. Wenn man sich etwas vornimmt, heisst das, dass man einem Thema Aufmerksamkeit schenkt und sich die vorzunehmende Handlung imaginiert. Dadurch kann dieser Vorsatz in der Erinnerung mit der Erinnerung an die tatsächlich vorgenommene Handlung verwechselt werden.
[2] Erst recht ist natürlich das Dekrementieren des Index mit dem Präfix-Operator statt mit der üblichen Postfix-Notation nicht durch Performancegründe zu rechtfertigen. Auch dies beleidigt nur das Auge - weiter nichts!