Dienstag, 6. Januar 2009

Portierung von Java-Code nach Assembler

Ich komme auf das Thema meines früheren Blogs Perl als Portierungshilfe zurück: Dort zeigte ich, wie man mit Hilfe von Perl Code von Java nach Assembler portieren kann.

Das funktionierte zwar gut, hatte jedoch einen Schönheitsfehler, den man bei Codegeneratoren gern vernachlässigt: Die Lesbarkeit des Zielcodes. Auch automatisch generierter Code ist Quellcode und sollte so verfasst sein, dass er für Menschen gut nachvollziehbar bleibt. Die langen trigonometrischen Reihen der Newcombschen Sonnenformeln wurden vom Codegenerator direkt in FPU-Anweisungen übersetzt. Wenn ich diese mit Cut & Paste in mein Assemblerprogramm einfügte, nahmen die Reihen fast ein Drittel des ganzen Listings in Anspruch. Das Auge ermüdete beim Lesen dieser Codestrecken wegen der immer gleichen Anweisungen. Der Blick auf das Wesentliche war erschwert:
; --- perturbations in longitude
; dl2 = 4.838*cos(299.102e0 + g2 - g)
fld FP8(5.22031469929999)
fadd real8 ptr nd.g2+0
fsub real8 ptr nd.g+0
fcos
fmul FP8(0.0844390292114842)

; + 0.116*cos(148.900e0 + 2*g2 - g)
fld FP8(2.59879525621951)
fadd real8 ptr nd.g2+8
fsub real8 ptr nd.g+0
fcos
fmul FP8(0.00202458193231339)
fadd
; + 5.526*cos(148.313e0 + 2*g2 - 2*g)
fld FP8(2.58855017351031)
fadd real8 ptr nd.g2+8
fsub real8 ptr nd.g+8
fcos
fmul FP8(0.096446894465205)
fadd
; + 2.497*cos(315.943e0 + 2*g2 - 3*g)
fld FP8(5.51424559862835)
fadd real8 ptr nd.g2+8
fsub real8 ptr nd.g+16
fcos
fmul FP8(0.0435808714222977)
fadd


Und so weiter, seitenlang. So wird der Quelltext zur Bleiwüste! Da dieser Quellcode für einen menschlichen Leser gedacht ist, sollte er sich auf die wesentliche Information beschränken. Was aber wäre das Wesentliche dieser Anweisungen? Die Koeffizienten und die involvierten Argumente. Alles andere bleibt ja für jeden Term gleich.

Ich hatte mich mittlerweile entschieden, die bogensekundengenaue Sonne aus dem Programm zu entfernen, da die Berechnung selbst auf Assemblerebene mir noch zu lange dauerte. Stattdessen verwende ich für die Sonne nun die Keplerellipse der Erde, was für die insgesamt angestrebte Bogenminutengenauigkeit voll ausreicht. Das Problem mit den trigonometrischen Reihen stellte sich erst wieder beim Mond. Denn die Mondbahn ist so speziell, dass man mit einer wie auch immer zeitlich veränderten Ellipse als Bahnmodell keine bogenminutengenauen Resultate erhält. Hier war also wieder eine Fourier-Entwicklung nötig, wie sie Ernest William Brown (1866-1938) als erster geleistet hat.

Wie bei der Sonne, hatte ich auch diese Reihenentwicklung bereits in der Java-Klasse LowPrecCalculator implementiert. Wieder stellte sich die Aufgabe, den bestehenden Code von Java nach Assembler zu portieren, und wieder half ein Perl-Programm. Aber es generiert nun nicht mehr direkt die FPU-Anweisungen, sondern Macro-Aufrufe. Die trigonometrischen Terme der Vorlagefunktion sehen etwa so aus:
dl = 22639*sin(u)
-4586*sin(u-2d)
+2370*sin(2d)
+769*sin(2u)
-668*sin(v)
-412*sin(2f)
+...

Die Subroutine moon macht etwa folgendes daraus:
; First element: Ecliptical longitude
fld md.u
fsin
fimul INT2(22639)
addTrigoTerm -4586,1,u,-2,d
addTrigoTerm 2370,2,d
addTrigoTerm 769,2,u
addTrigoTerm -668,1,v
addTrigoTerm -412,2,f
...

Die Generierung der FPU-Anweisung ist somit an das Macro addTrigoTerm delegiert. Nur die Argumente, die ja das Wesentliche des Terms enthalten, stehen noch im Quelltext der Routine. Das Macro mit den FPU-Anweisungen steht dagegen nur noch einmal im Code. Das reduziert den benötigten Platz und hat zu einer deutlich lesbareren Version der Mondroutine geführt.

Das Macro addTrigoTerm ist folgendermassen definiert (ich bin mittlerweile für dieses Projekt auf englische Kommentare umgestiegen, da ich plane, es bei sourceforge.net zu veröffentlichen und frei verfügbar zu machen):
addTrigoTerm macro factor, _args:vararg
$cnt = 0
$first = 1
$fact = 0
$off = 0
irp $v, <_args>
$cnt = $cnt + 1
if $cnt eq 1
$fact = $v ; remember arg for next pass
if $v lt 0
$off = ( -$v - 1 ) * 8 ; resulting byte offset
else
$off = ( $v - 1 ) * 8 ; resulting byte offset
endif
else
$cnt = 0 ; reset counter for next pair
if $fact lt 0
if $first eq 1
fldz
$first = 0
endif
fsub real8 ptr md.$v+$off
else
if $first eq 1
fld real8 ptr md.$v+$off
$first = 0
else
fadd real8 ptr md.$v+$off
endif
endif
endif
endm
fsin
fimul INT2(factor)
fadd
endm

Der Vollständigkeit halber sei noch die neue Version des Perl-Programms angefügt. Es generiert nicht nur die trigonometrischen Terme, sondern fügt auch noch eine Statistik an, welche Argumente mit welcher maximalen Multiplizität benötigt wurden.
# Java -> ASM converter for the moon series 
# Output in form of macro calls
use strict;
use warnings;

my ($line,$arg,$coeff,$dep,$fact,$stmt);
my %maxfact = map { $_ => 1 } qw(d f u v);

foreach $line (<DATA>) {
if (($coeff,$arg) = $line =~ /([+-]?\d+)\*sin\((.+)\)/) {

# Eliminate leading '+'
$coeff *= 1;

# Start statement
$stmt = "addTrigoTerm $coeff";

# Progressive matching of the sin args
while ($arg =~ /([+-]?\d*)([dfuv])/g) {
($fact,$dep) = ($1,$2);
# $fact may assume the special values "+", "-" or ""
$fact = $fact."1" unless $fact =~ /\d/;
# Eliminate a leading "+"
$fact *= 1;
# Get statistics: Maximal factor for each dependency
$maxfact{$dep} = abs($fact) if $maxfact{$dep} < abs($fact);
# Next part of the statement
$stmt .= ",$fact,$dep";
}

# Print it out
print "$stmt\n";
}

else {
print "; $line";
next;
}
}

print "; Statistics:\n";
foreach $dep (keys %maxfact) {
print "; $dep:$maxfact{$dep}\n";
}

__DATA__
dl
+22639*sin(u)
-4586*sin(u-2d)
+2370*sin(2d)
+769*sin(2u)
-668*sin(v)
-412*sin(2f)
...

Keine Kommentare :