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).

Donnerstag, 26. Dezember 2013

Backup und kurzfristige Versionierung mit der Dropbox

Da ich dieses Projekt nicht im Team bearbeite, sondern solo, ist ein "richtiges" Repository ein bisschen mit Kanonen auf Spatzen schiessen. Eigentlich reicht vorerst die auf dreissig Tage befristete Versionierung, wie sie die Dropbox-Gratis-Accounts anbieten.

In Linux gibt es die schönen symbolischen Links. Ich habe einen Ordner angelegt, der symbolische Links auf die Dateien enthält, die von der SwissEphemeris-Distribution abweichen:

#!/usr/bin/perl -w

use strict;

my @syncfiles = qw(
  test_calc_reduce.c
  test_calc_reduce.h
  test_test_calc_reduce.c
  test_test_calc_reduce.h
  testgen.c
  testgen.h
  tcr.fix
  tests.cpp
  tests.hpp
  makefile
  );

foreach my $file (@syncfiles) {
  `ln -s ~/workspace/swepar/src/$file $file`;
  }
... und hier ist der

Link auf den Dropbox-Ordner