Mittwoch, 22. Januar 2014

Initialisierung ausschalten

Das Ausschalten der Initialisierung
    if (0) {
    swe_close( );
    const struct swe_data swed0 = {
      .ephe_path_is_set = false,
      .jpl_file_is_open = false,
      .fixfp = NULL,
      .ephepath =  SE_EPHE_PATH,
      .jplfnam = SE_FNAME_DFT,    /* JPL file name, default */
      .geopos_is_set = false,   /* geopos is set, for topocentric */
      .ayana_is_set = false     /* ayanamsa is set */
      };
    memcpy(&swed,&swed0,sizeof(struct swe_data));
    }
treibt die Abdeckung um ein halbes Prozent in die Höhe. Vermutlich ist es anderer Code, der in diesen Fällen durchlaufen wird. Die Tests bleiben aber alle auf OK, was beruhigend ist: Es gibt anscheinend keine weiteren verschleppten Abhängigkeiten mehr, die die Ergebnisse der Tests beeinflussen könnten.

Sonntag, 12. Januar 2014

82 Prozent

Ich habe nun den Testgenerator erweitert, so dass auch Testdaten für siderische Positionen generiert werden können. Ebenso kann eine topozentrische Position in der Fixture angegeben werden, so dass topozentrische Planetenberechnungen in der Fixture zugänglich sind. Für die Fixsterne und ihre Grössen, für die Ayanamsa-Namen, Planetennamen und die Zeitgleichung erlaube ich mir, die Tests und Erwartungen in test_test_calc_reduce.c hart zu codieren. Sie kommen im Testprotokoll immer am Schluss, nach den Tests aus der Fixture

Und so sieht das Protokoll nun aus.

ruediger@herschel:~/workspace/swepar/src$ test_calc_reduce
geocentric positions: 250 OK   0 NOT OK
helio/barycentric positions: 285 OK   0 NOT OK
barycentric Sun: 5 OK   0 NOT OK
astrometric positions J2000: 90 OK   0 NOT OK
topocentric positions: 15 OK   0 NOT OK
sidereal positions with Lahiri Ayanamsha: 25 OK   0 NOT OK
sidereal positions / Suryasiddhanta, strictly computed: 25 OK   0 NOT OK
individual ayanamsa: 25 OK   0 NOT OK
sidereal rectangular geocentric: 40 OK   0 NOT OK
Asteroid Urania: 5 OK   0 NOT OK
Hamburger Poseidon: 5 OK   0 NOT OK
different coordinate systems: 4 OK   0 NOT OK
Obliquity of ecliptic and nutation: 1 OK   0 NOT OK
Moshier, geocentric positions: 140 OK   0 NOT OK
Moshier, heliocentric positions: 100 OK   0 NOT OK
Provoking a "file not found" error: 1 OK   0 NOT OK
swe_calc_ut: 1 OK   0 NOT OK
New mode: Barycentric, directly from file: 1 OK   0 NOT OK
Fixed Stars: 7 OK   0 NOT OK
Get Ayanamsa: 6 OK   0 NOT OK
Get Planet Names: 3 OK   0 NOT OK
Equation of Time for J2000: 1 OK   0 NOT OK
Die damit erreichte Testabdeckung liegt bei 82 Prozent. In dem arg zerklüfteten gcov-File gibt es jetzt keine Hauptfunktion mehr, die nicht abgedeckt wäre. Die vielen kleinen roten Einsprengsel entstehen:
  • durch zahllose nicht durchlaufene Ausnahmesituationen (z.B. "Ephemeride nach Moshier darf nicht benutzt werden für baryzentrische Koordinaten"). Einiges davon hat mit nicht mehr unterstützten Optionen zu tun.
  • durch nicht ausgenutzte Puffer: Da ich zwischen den Tests normalerweise nun swe_close() aufrufe, wird der Code zum Wiederverwenden bereits berechneter Positionen nicht durchlaufen. Dies könnte ich noch mit einer weiteren kleinen Testserie erreichen.

Freitag, 10. Januar 2014

Der Nageltest für die Local History

Gleich nach meinem gestrigen Lobpreis auf die Local History konnte ich nun die Probe aufs Exempel machen. In der Konsole hatte ich versehentlich ein rm test_calc* ausgeführt, was mir ohne Sicherungen einige Backup-Arbeit beschert hätte. Nun konnte ich auf die Local History zurückgreifen: Sie erschien wieder, nachdem ich die Files test_calc_reduce.c und test_calc_reduce.h wieder leer angelegt hatte.

Nach einigen kleineren Macken konnte ich nun fünf Fixsterntests und vier Ayanamsa-Tests hinzufügen, was die Abdeckung zwar nicht weltbewegend hochgetrieben hat (sie liegt stieg von 68% auf 71% an), aber mit swe_fixstar() und swe_deltat() sind nun zwei prominente Funktionen mit im Rennen.

Die folgende Übersicht zeigt,

  • dass im Bereich der Berechnung von Fixsternpositionen und -grössen noch besonders viel Luft ist.
  • lunar_osc_elem könnte mit ein paar Lilith-Berechnungen besser ausgebeutet werden.
  • sweph, die zentrale Funktion für die Planetenberechnungen, ist erst zu drei Vierteln abgedeckt, obwohl sie mit den bestehenden Tests bereits über 6000 Mal aufgerufen wird.
  • Der Code, der die Anfrage aus dem Puffer bedient, wird nicht durchlaufen, da ich bislang keine zwei identische Positionsanfragen nacheinander stelle.
  • Ausserdem werden keine Asteroiden abgefragt, für die eines der ast-Subdirectories abgefragt werden müsste.
  • Die Abdeckung lässt sich deutlich steigern, wenn noch ein paar topozentrische Berechnungen dazukommen (vor allem durch swid_get_observer).
  • Die Zweige, die Planetenberechnungen mit Ayanamsa enthalten, werden noch nicht durchlaufen.

Donnerstag, 9. Januar 2014

Local History

Hätte ich die Local History von Eclipse schon 2012 gekannt, dann hätte ich nicht folgende Frage im Forum der Clean Code Developer abgesetzt:

Delta-Backup für alle Änderungen an Textdateien

Genau die Umkehrung des üblichen Master/Slave-Verhältnisses leistet die Local History, die für mich ohne weiteres Zutun von Eclipse verwaltet wird. Alles, was ich brauche, kann ich damit machen: Versionen vergleichen (und zwar jeden beliebigen Stand, ohne dass ich Commits einchecken und mit einem Text versehen muss), und bei Bedarf auch alte Versionen wieder zurückholen.

Wenn ich will, kann ich sogar einzelne Deltas aus früheren Versionen wieder holen. Alles sehr komfortabel und nützlich.

Mittwoch, 8. Januar 2014

Über die Benutzung vorgesehener Funktionen. Coverage jetzt 68%

Stand gestern: In der bash ausgeführt, lieferte das Programm test_calc_reduce unsinnige Testergebnisse, während es in Eclipse mit dem Pfad Run as->Local C/C++ Application korrekt ausgeführt wurde. Wie ist dieser Unterschied zu erklären?

Zwischen den Tests hatte ich ja die Struktur swed explizit initialisiert.

    const struct swe_data swed0 = {
      .ephe_path_is_set = false,
      .jpl_file_is_open = false,
      .fixfp = NULL,
      .ephepath =  SE_EPHE_PATH,
      .jplfnam = SE_FNAME_DFT,    /* JPL file name, default */
      .geopos_is_set = false,   /* geopos is set, for topocentric */
      .ayana_is_set = false     /* ayanamsa is set */
      };
    memcpy(&swed,&swed0,sizeof(struct swe_data));
Indem ich nun dieser Initialisierung noch einen Aufruf von swe_close() voranstellte, funktionierte das Programm mit beiden Aufrufmethoden wieder gleich. Es lohnt sich vermutlich nicht, der genauen Ursache nachzugehen, warum es zu diesem Fehlverhalten kam. Denn selbst wenn ich es vergessen haben sollte, ein globales Datum zu löschen: Warum verhielt sich das Programm dann beim Aufruf aus Eclipse anders als beim Aufruf aus der Konsole?

Anderes Thema: Ich habe einmal den Coverage Analyzer gcov angeworfen und festgestellt, dass ich mit der Testabdeckung für sweph.c (also mein test_calc_reduce.c) schon ganz gut liege:

Die vertikale Seitenleiste in der Abdeckungs-Quelltextanzeige zeigt, dass es vor allem im unteren Teil des Files noch viel Rot gibt: da sind einige Funktionen, die noch gar nicht aufgerufen werden. Das topozentrische Koordinatensystem habe ich noch mit keinem Testfall angefordert. Die vielen kleinen Splitter zwischendurch sind Ausnahmefälle wie: Die Fehlermeldung bei der Chironberechnung für ein Datum kleiner als CHIRON_START. Oder ein File, das behauptet, Big Endian zu sein, aber es in Wahrheit nicht ist, weil schon seine eigene Grösse in falscher Codierung abgespeichert ist. Hier werde ich irgendwo einen Stopp machen. Ich könnte natürlich in einem Test ein präpariertes File in den Ephemeridenpfad schieben, das dort das entsprechende korrekte File ersetzt. So könnte ich auch viele dieser Splitter noch "erschlagen".

Eine 100%-Abdeckung bekomme ich voraussichtlich nicht hin. Aber ab 95% könnte ich mit gutem Gefühl damit beginnen, auch aggressivere Umbaumassnahmen in meinem "Branch" test_calc_reduce.c vorzunehmen. Und mir dabei ziemlich sicher zu sein, dass die Kernfunktionen der Swiss Ephemeris erhalten bleiben.

Dienstag, 7. Januar 2014

Statische Variablen

Gestern abend räsonierte ich noch über rätselhafte Fixsterne. Das Rätsel der zwei verschiedenen Positionen desselben Fixsterns (Aldebaran) war schnell gelöst.

Es ist so, wie ich gestern schon schrieb:

Das kann nur bedeuten, dass es noch irgendeinen anderen, nicht in swed enthaltenen, Zustand im Programm gibt. Das werde ich morgen weiter analysieren.
Und dieser verborgene Zustand bestand aus zwei in swe_fixstar deklarierten statischen Variablen:
  static char slast_stardata[AS_MAXCH];
  static char slast_starname[AS_MAXCH];
die zur Pufferung des zuletzt gelesenen Satzes im Fixsternkatalog verwendet wurden.

Das erklärt, dass der innere Zustand des Programms beim zweiten, identischen Aufruf ein anderer war als beim ersten Aufruf. Nun waren die beiden statischen Variablen - noch vom ersten Aufruf - gefüllt.

Wenn es aber nur eine Pufferung war - warum ergab der zweite Aufruf dann andere Werte als der erste?

Der Grund ist, dass ich durch das Löschen der swed (z.B. mittels eines swe_close()-Aufrufs, oder - wie in meinem Testprogramm: explizit) auch die Information swed.is_old_starfile = TRUE zerstört hatte, die beim ersten Aufruf ermittelt worden war (in meinem Ephemeridenverzeichnis befand sich nämlich nur das "alte" File fixstars.cat - was aber ja nicht verboten ist).

  if (swed.fixfp == NULL) {
    if ((swed.fixfp = swi_fopen(SEI_FILE_FIXSTAR, SE_STARFILE, swed.ephepath, serr)) == NULL) {
      swed.is_old_starfile = TRUE;
      if ((swed.fixfp = swi_fopen(SEI_FILE_FIXSTAR, SE_STARFILE_OLD, swed.ephepath, NULL)) == NULL) {
        ...
Wenn der Eintrag im Puffer gefunden wurde, springt das Programm direkt zum Label found. Nach Erreichen der Sprungmarke gibt es aber noch eine Stelle, in der mit den Koordinaten abhängig vom Flag swed.is_old_starfile unterschiedlich weitergerechnet wird. Das erklärt's:
  /* speed in ra and de, degrees per century */
  if (swed.is_old_starfile == TRUE) {
    ra_pm = ra_pm * 15 / 3600.0;
    de_pm = de_pm / 3600.0;
  } else {
    ra_pm = ra_pm / 10.0 / 3600.0;
    de_pm = de_pm / 10.0 / 3600.0;
    parall /= 1000.0;
  }
Nach Ausbau dieser Pufferung und Entfernung der beiden statischen Variablen slast_stardata und slast_starname liefen die beiden Fixsterntests nun beide mit OK durch. Insbesondere müssen die beiden sukzessiven Aufrufe nun auch dieselben Ergebnisse gebracht haben.

Nachdem dies erledigt war, habe ich auch das neuere File sefstars.txt ins Ephemeridenverzeichnis hinzugefügt. Im Grunde ist fixstars.cat ja ein alter Zopf, den ich irgendwann kappen kann.

Wie ich ja früher auch schon den Compileparameter TRACE eliminiert habe.

Das Tool nm, das das Symbolverzeichnis eines Objektfiles auflistet, zeigt, dass es noch zwei weitere statische Variablen in sweph.c gibt: epheflag_sv und nutflag:

ruediger@herschel:~/workspace/swepar/src$ nm -aC sweph.o | grep ' [bB] ' 
00000000 b .bss
00000204 b epheflag_sv.3938
00000200 b nutflag.5350
00000000 b slast_stardata.5374
00000100 b slast_starname.5375
Um diese werde ich mich morgen kümmern.

Verwirrend kommt nun hinzu, dass der Start des Programms test_calc_reduce in der Bash nicht mehr funktioniert, wohl aber der Start von Eclipse aus mittels Run As -> Local C/C++ Application

Montag, 6. Januar 2014

Rätselhafte Fixsterne

Das Programm testgen.c funktioniert wunderbar für die Erzeugung von Test Fixtures, die den Aufruf der Funktion swe_calc betreffen. Nun gibt es aber noch eine ganze Reihe weiterer öffentlicher Funktionen, die für eine volle Testabdeckung des Codes aufzurufen wären. Hierzu hätte man testgen.c durch andere Testtypen erweitern müssen.

Leider war testgen.c nur halbherzig für solche Erweiterungen vorbereitet. Es gibt zwar im Header einer Testreihe ein Feld type. Aber damit hat es sich auch schon. Um es wirklich erweiterbar zu halten, müsste für jeden Testtyp ein Set von drei Funktionen plus eine Datenstruktur bereitgestellt werden:

  • Eine Funktion zum Einlesen der textförmigen Fixture in die Datenstruktur;
  • Eine Funktion, die mit gegebenen Werten in dieser Datenstruktur mittels der Originalbibliothek libswe.so Ausgabewerte erzeugt;
  • Eine Funktion, die mit gegebenen Werten in dieser Datenstruktur die entsprechende Funktion im Prototyp test_calc_reduce.c aufruft und die Ergebnisse mit den durch die vorherige Funktion berechneten (und in der binären Fixture abgespeicherten) Ergebnissen des Standards vergleicht.
So etwas wäre machbar. Aber aus Zeitgründen verschiebe ich das nach hinten und schreibe lieber erstmal einige ad hoc-Tests für die noch verbliebenen Funktionen. Allen voran swe_fixstar().

Hier mache ich nun die merkwürdige Erfahrung, dass die swe_fixstar(), zweimal nacheinander mit denselben Werten aufgerufen, zu unterschiedlichen Ergebnissen führt, obwohl ich die globale Datenstruktur swed zwischen den beiden Aufrufen lösche.

static bool test_fixstar() {

  result r = {};
  const double prec = 1e-7;
  const struct swe_data swed0 = {};

  testdata_fixstar td[] = {
      {
        .star ="Aldebaran",
        .jd = 2451544.5,
        .result_exp = {
           .pos = {
            l:69.7903072,
            b:-5.4676328
          }
        }
      },
      {
        .star ="Aldebaran",
        .jd = 2451544.5,
        .result_exp = {
           .pos = {
            l:69.7903072,
            b:-5.4676328
          }
        }
      }
    };
  const int TESTCOUNT = 2;

  for (int i=0;i<TESTCOUNT;i++) {
    memcpy(&swed,&swed0,sizeof(struct swe_data));
    bool ok = test_single_fixstar( td[i], prec );
    if (ok) r.ok++; else r.not_ok++;
    }

  summary("Fixed Stars",r);

  return r.not_ok == 0;
  }

static bool test_single_fixstar( testdata_fixstar td, double prec) {

  double xx_act[6];
  char serr_act[255], sname[41];

  strcpy(sname,td.star);

  int iflag_act = swe_fixstar( sname, td.jd, td.iflag_in, xx_act, serr_act );

  td.result_exp.pos.r = 1;
  td.result_exp.speed.r = 0;
  if (!td.serr_exp) td.serr_exp = "";

  bool ok;
  ok = diff_check_int( "iflag", iflag_act, td.iflag_exp  );
  ok = diff_check_degree( "pos.l", xx_act[0], td.result_exp.pos.l, prec ) && ok;
  ok = diff_check_degree( "pos.b", xx_act[1], td.result_exp.pos.b, prec ) && ok;
  ok = diff_check( "pos.r", xx_act[2], td.result_exp.pos.r, prec ) && ok;
  ok = diff_check_degree( "speed.l", xx_act[3], td.result_exp.speed.l, prec ) && ok;
  ok = diff_check_degree( "speed.b", xx_act[4], td.result_exp.speed.b, prec ) && ok;
  ok = diff_check( "speed.r", xx_act[5], td.result_exp.speed.r, prec ) && ok;
  ok = diff_check_char( "serr",  serr_act,  td.serr_exp ) && ok;

  return ok;

  }
Das kann nur bedeuten, dass es noch irgendeinen anderen, nicht in swed enthaltenen, Zustand im Programm gibt. Das werde ich morgen weiter analysieren.

Sonntag, 5. Januar 2014

Flags mit GNU bc aufschlüsseln

Mehrmals schon hatte ich beim Testen oder Debuggen das Bedürfnis, die Bestandteile eines iflag-Wertes in seine einzelnen Optionen zu zerlegen, ohne ihn erst (z.B. im Kopf) binär umrechnen zu müssen.

Für solche kleinen Rechnereien empfiehlt sich bc(1), ursprünglich ein kleiner Rechner für die Kommandozeile, und von Entwicklern weiter aufgepeppt. bc(1) kann auch Funktionen, hat Kontrollstrukturen und Schleifen, kennt Arrays, symbolische Variable usw.

Hier mein Tool iflag.

#!/usr/bin/bc -q 

scale = 0

print "iflag = "
iflag = read()

define bit( i, p ) {
  return ((i/p) % 2 == 1)
  }

if (bit( iflag, 1)) print "seflg_jpleph\n";
if (bit( iflag, 2)) print "seflg_swieph\n";
if (bit( iflag, 4)) print "seflg_moseph\n" ;
if (bit( iflag, 8)) print "seflg_helctr\n";
if (bit( iflag, 16)) print "seflg_truepos\n";
if (bit( iflag, 32)) print "seflg_j2000\n";
if (bit( iflag, 64)) print "seflg_nonut\n";
if (bit( iflag, 128)) print "seflg_speed3\n";
if (bit( iflag, 256)) print "seflg_speed\n";
if (bit( iflag, 512)) print "seflg_nogdefl\n";
if (bit( iflag, 1024)) print "seflg_noaberr\n";
if (bit( iflag, (2*1024))) print "seflg_equatorial\n";
if (bit( iflag, (4*1024))) print "seflg_xyz\n";
if (bit( iflag, (8*1024))) print "seflg_radians\n";
if (bit( iflag, (16*1024))) print "seflg_baryctr\n";
if (bit( iflag, (32*1024))) print "seflg_topoctr\n";
if (bit( iflag, (64*1024))) print "seflg_sidereal\n";
if (bit( iflag, (128*1024))) print "seflg_icrs\n";
if (bit( iflag, (256*1024))) print "seflg_new_mode\n";

quit
Und hier eine Beispielsitzung:
ruediger@herschel:~/workspace/swepar/src$ iflag 
iflag = 16912
seflg_truepos
seflg_nogdefl
seflg_baryctr

Nachtrag

Gegen eine Scriptsprache, die sich gewaschen hat, wie z.B. Perl, kann bc natürlich nicht ankommen. Das ist auch gar nicht die Absicht. bc will nicht viel mehr als eine Art fortgeschrittener Taschenrechner sein, der in die Unix-Philosophie integriert ist, d.h. sich komfortabel mit anderen Befehlen in der Kommandozeile kombinieren lässt.

Hier das selbe Beispiel in Perl:

#!/usr/bin/perl

my @SEFLG=qw(SEFLG_JPLEPH  SEFLG_SWIEPH    SEFLG_MOSEPH 
             SEFLG_HELCTR  SEFLG_TRUEPOS   SEFLG_J2000   
             SEFLG_NONUT   SEFLG_SPEED3    SEFLG_SPEED   
             SEFLG_NOGDEFL SEFLG_NOABERR   SEFLG_EQUATORIAL 
             SEFLG_XYZ     SEFLG_RADIANS   SEFLG_BARYCTR 
             SEFLG_TOPOCTR SEFLG_SIDEREAL  SEFLG_ICRS   
             SEFLG_NEW_MODE);
my $iflag = $ARGV[0];
foreach my $flag (@SEFLG) {
  print "$flag\n" if $iflag % 2;
  last unless $iflag >>= 1;
  }

Samstag, 4. Januar 2014

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!

Freitag, 3. Januar 2014

Heliozentrisch!

Der Grund dafür, dass meine experimentellen Funktion calcMainPlanetFromFile() für die inneren Planeten fehlschlägt, scheint der zu sein, dass diese im File helio- und nicht baryzentrisch abgespeichert sind. Der Debugger zeigt, dass bei einem Vergleichsaufruf von swetest mit denselben Daten in Zeile 1601ff. von sweph.c in diesem Fall noch die baryzentrische Sonnenposition zur errechneten (heliozentrischen) Position addiert wird:
      /* if planet is heliocentric, it must be transformed to barycentric */
      if (pdp->iflg & SEI_FLG_HELIO) {
 /* now barycentric planet */
 for (i = 0; i <= 2; i++) 
   xp[i] += xps[i];
 if (do_save || (iflag & SEFLG_SPEED))
   for (i = 3; i <= 5; i++) 
     xp[i] += xps[i];
      }
Im Debugger sieht man, dass vor dieser Korrektur die errechneten Werte noch exakt mit denen aus meiner experimentellen Funktion übereinstimmen (im Konsolenfenster im Hintergrund des Screenshots sieht man die Soll/Ist-Abweichungen des Tests).

Die Flagleiste pdp->iflg gehört dabei zur Struktur plan_data, die beim Lesen des Files gefüllt wird. Die Information, ob die Daten heliozentrisch oder baryzentrisch sind, ist also nicht statisch im Programm definiert, sondern mit den Daten zusammen im File abgespeichert.

Habe mich auch wieder daran erinnert, dass gdb-Breakpoints in dynamischen Bibliotheken erst dann gesetzt werden können, wenn die Bibliothek innerhalb der Debuggingsitzung das erste Mal geladen wurde.

pdfposter

Da der Aufrufgraph von sweph.c (siehe letzter Post) mir im Druck immer auf A4 reduziert wurde, so dass selbst ein weniger alterssichtiger Mensch als ich die Texte der einzelnen Knoten nicht mehr lesen kann, suchte ich ein Programm für den Posterdruck und fand Hartmut Goebels pdfposter.

Mit diesem Programm gelang mir die Aufteilung des PDF-Dokuments auf vier A4-Seiten, die zusammen ein A2-Poster im Hochformat ergeben. Beim Ausdruck ergaben sich allerdings Ränder, auf denen die Farbe "schmierte". Ob es am Treiber oder am pdfposter-Programm selbst liegt, war mir zu mühsam herauszufinden. Ich rückte den vier Seiten einfach mit der Schere zu Leibe und habe nun mein sweph.c-Poster.

Aufrufgraph auf vier A4-Seiten

Donnerstag, 2. Januar 2014

Ein Aufrufgraph für sweph.c

Das Kernstück der Swiss Ephemeris ist das Includefile sweph.c. Um die Kohäsion und die Kopplung der einzelnen Funktionen einzuschätzen und das Programm insgesamt besser zu verstehen, hatte ich das Bedürfnis nach einem Aufrufgraphen. Es soll für ein bestimmtes C-Sourcefile (natürlich sweph.c, das fetteste) ein Diagramm erstellt werden, das für jede Funktion zeigt, welche anderen Funktionen sie jeweils aufruft.

Der in Eclipse CDT vorgesehene View Call Graph klingt vom Namen her verheissungsvoll, bricht aber bei Aufruf mit folgendem Stacktrace ab:

org.eclipse.core.runtime.AssertionFailedException: null argument:Action must not be null
 at org.eclipse.core.runtime.Assert.isNotNull(Assert.java:85)
 at org.eclipse.jface.action.ContributionManager.add(ContributionManager.java:75)
 at org.eclipse.linuxtools.internal.callgraph.CallgraphView.createPartControl(CallgraphView.java:532)
 at org.eclipse.ui.internal.e4.compatibility.CompatibilityPart.createPartControl(CompatibilityPart.java:138)
 at org.eclipse.ui.internal.e4.compatibility.CompatibilityView.createPartControl(CompatibilityView.java:155)
 at org.eclipse.ui.internal.e4.compatibility.CompatibilityPart.create(CompatibilityPart.java:313)
        ...
Der Hauptgrund für den Abbruch sind nicht etwa die in meinem System fehlenden Komponenten wie System Tap (ein interessantes Programm, wie mir scheint), sondern dass der View Call Graph nicht statisch, sondern dynamisch gebildet wird. Dieser Call Graph ist ein Werkzeug fürs Profiling, zum Studium des Laufzeitverhaltens – und nicht das statische Analysetool, nach dem ich suche. Also suchte ich nach einer Alternative.

Das Projekt besteht aus einer Reihe von C-Dateien, die wie üblich jede einzeln in eine gleichnamige Objektdatei compiliert wird. Um einen Aufrufgraphen zu erstellen, müssen wir wissen,

  • welche Funktion welche anderen Funktionen aufruft, und
  • welche Funktionen in welchem Include enthalten sind, um die Funktionen als Gruppe markieren zu können.

Listen der in einem Include enthaltenen Funktionen kann man aus dem jeweils generierten Objectfile (*.o) herauslesen, das ja eine Symboltabelle der compilierten Funktionen enthält.

Funktionen, die mittels Compilerschaltern ausgeblendet sind, werden bei diesem Vorgehen allerdings ignoriert, aber das stört mich hier nicht besonders. Falls das in anderen Situationen stören sollte, kann man auch CTags verwenden, um den Quellcode selbst auf Funktionsdefinitionen zu durchforsten.

Mit dem GCC-Tool objdump kann diese Symboltabelle mit der Option -t ausgelesen werden. Für sweph.o beginnt der Symboltabellen-Dump z.B. wie folgt:

ruediger@herschel:~/workspace/swepar/src$ objdump -t sweph.o

sweph.o:     file format elf32-i386

SYMBOL TABLE:
00000000 l    df *ABS* 00000000 sweph.c
00000000 l    d  .text 00000000 .text
00000000 l    d  .data 00000000 .data
00000000 l    d  .bss 00000000 .bss
00000000 l    d  .rodata 00000000 .rodata
00000000 l     O .rodata 000000a8 pla_diam
000000c0 l     O .rodata 000001e0 ayanamsa
00000000 l    d  .data.rel.local 00000000 .data.rel.local
00000000 l     O .data.rel.local 0000006c ayanamsa_name
00000420 l     O .rodata 0000002c pnoint2jpl
00000460 l     O .rodata 00000054 pnoext2int
00000204 l     O .bss 00000004 epheflag_sv.3938
000005b1 l     F .text 00001a56 swecalc
0000f554 l     F .text 00000165 denormalize_positions
0000f6b9 l     F .text 000000d6 calc_speed
0000faa0 l     F .text 000000c4 plaus_iflag
0000395c l     F .text 000004e1 jplplan
00003169 l     F .text 000007f3 sweplan
In der Symboltabelle wird so allerhand Mystisches aufgeführt; aber eben auch alle enthaltenen Funktionen; diese sind dadurch zu erkennen, dass sie in der Liste das Kürzel F haben und ihr Name mit .text eingeleitet wird.

Mit der folgenden Perl-Routine kann ich alle Object Files eines Verzeichnisses auf einmal durchgehen und eine Hashtabelle erstellen, die für jede Funktion die Nummer des Object Files liefert, in der sie vorhanden ist:

sub get_all_symbols {

# Generates a hash containing all symbols 
# of all object files in this directory

  my %symbols;
  my @obj = split /\s+/, `ls *.o`;
  my $count = -1;
  my $exceptions = join "|", @ignore_object_files;
  for my $fname (@obj) {
    next if $fname =~ /$exceptions/; # Skip the exceptions
    $count++;
    foreach (`objdump -t $fname`) {
      chomp;
# This is very GCC specific, expecting the OBJDUMP output format
      next unless /^[0-9a-f]{8}.*?F\s+\.text\s+[0-9a-f]{8}/;
      my $name = $_;
      $name =~ s/^.*?(\w+)$/$1/;
      $symbols{$name} = $count;
      }
    }
  return \%symbols;
  }
Hierbei kann man mit dem Array @ignore_object_files noch eine Reihe von Object Files aufführen, die nicht analysiert werden sollen.

Für die eigentliche Aufrufanalyse gibt es das Programm cflow. Mit der Option -l aufgerufen liefert es die gewünschte Vorwärtsreferenztabelle: Für jede Funktion werden alle von dieser Funktion aufgerufenen Funktionen gelistet. Für sweph.c beginnt die Liste beispielsweise so:

ruediger@herschel:~/workspace/swepar/src$ cflow -l sweph.c | head
{   0} swe_calc() :
{   1}     memset()
{   1}     fopen()
{   1}     fgets()
{   1}     strchr()
{   1}     atol()
{   1}     fclose()
{   1}     swi_open_trace()
{   1}     trace_swe_calc() <void trace_swe_calc (int swtch, double tjd, int ipl, int32 iflag, double *xx, char *serr) at sweph.c:5976>:
{   2}         fputs()
Wieder ist es so, dass bei einem automatischen Auswerten dieser Information bestimmte Aufrufe ignoriert werden sollten, hier z.B. die Aufrufe von Standardfunktionen wie fclose() usw.

Im Grunde sind damit alle Informationen beisammen, um den Graphen zu erstellen, was am besten auf dem Weg über die graphische Sprache dot von graphviz geschieht. Um die GraphViz-Befehle zu erforschen, gibt es übrigens eine schöne Webanwendung, den GraphViz Workspace. Er erzeugt bereits während der Eingabe das resultierende Diagramm.

Ich orientierte mich, was die dot-Syntax und Gestalt der zu erzeugenden Knoten angeht, an dem Script cflow2dot.pl.

Meine Version unterscheidet sich von der Vorlage im wesentlichen nur dadurch, dass es Filtermöglichkeiten für die Objektfiles und die Listen der aufgerufenen Funktionen gibt, und dass die verschiedenen Farbcodes dafür genutzt werden, die Funktionen nach ihrer Zugehörigkeit zu Source-Code-Includes zu gruppieren. Die Formen dagegen entsprechen der jeweiligen Aufruftiefe (so ist es auch im Originalskript).

Hier ist das Ergebnis für sweph.c (noch mit inkscape von SVG ins PDF-Format convertiert):

Die Trace-Funktionen sind gelb markiert, da sie in keinem Objektfile gefunden werden konnten (denn ich kompiliere das Projekt nicht mit dem Traceschalter). Auch konnten die Macro-Funktionen dot_prod und square_sum nicht gefunden werden. Aus Gründen der Compilersicherheit wäre es sowieso besser, diese Funktionen als Inline-Funktionen zu deklarieren (was ich in meinem Prototyp test_calc_reduce.c auch schon umgestellt habe - glaube ich).