Samstag, 15. Oktober 2011

Eine kleine Anfrage - und die Folgen

Am 22. September erhielt ich per Mail die Anfrage,

...ob Sie die Lösung von folgender Rechnung zufällig wissen:
1 2 3 4 5 6 7 8 9 = 2011

Die Zahlen müssen der Reihenfolge nach gebraucht werden, es darf mit +, -, · und : gerechnet werden und Klammern gesetzt.

Ich wünsche Ihnen einen schönen Abend!

Natürlich wusste ich die Lösung dieser Rechnung nicht. Ich wusste aber, dass Aufgaben dieser Art häufig in den Rätselecken von Zeitschriften zu finden und meist mit ein bisschen Getüftel und Probieren zu bewältigen sind – eine Tätigkeit, die von Laien oft mit Mathematik verwechselt wird, aber praktisch nichts mit ihr zu tun hat. Mathematik wäre es, wenn es einen "Top Down Approach" zur Lösung gäbe, man also aufgrund von systematischen, grundsätzlichen Überlegungen eine Lösung finden könnte. Das schien bei diesem Problem nicht der Fall zu sein.

Wen es als Tüftelei interessiert: Die Fussnote [1] zeigt einen Weg zu einer Lösung. Was man dabei nicht erfährt: gibt es eventuell noch mehr Lösungen? Der Lösungsweg ist nicht systematisch, sondern eine Art zielgerichtetes Probieren. Gibt es auch einen Ansatz, mit dem man mit Papier und Bleistift systematisch alle Lösungen herausfinden kann? (Vermutlich nicht.)

Da mir aber kein mathematischer systematischer Zugang zum Problem einfiel, wechselte ich schliesslich die Rolle: Ich schaute mir das Problem als Informatiker an. Wie praktikabel ist zum Beispiel die Holzhammermethode, alle möglichen Ausdrücke mit einem Computer auszurechnen? Mit wie vielen Möglichkeiten hat man es eigentlich zu tun?

Es fällt auf, dass die Zahl der Möglichkeiten offenbar ein Produkt ist: Die Zahl der möglichen Klammerungen multipliziert mit den Wahlmöglichkeiten für die acht Operationen. Letztere Anzahl ist leicht zu bestimmen: Da nur die vier Grundrechenarten zulässig sind und jede der acht Operationen unabhängig von den anderen gewählt werden kann, gibt es
48 = 65'536
Wahlmöglichkeiten für die Operationen.

Die möglichen Klammerungen führen zu den Baumdarstellungen arithmetischer Ausdrücke und schliesslich zu den sogenannten Catalanschen Zahlen. In diesem Fall kommt man auf
C8 = 1430

verschiedene Möglichkeiten, den Ausdruck zu klammern. Hier fragte ich mich - von diesem Problem ausgehend: Kann man die möglichen Klammerungen in Form eines Algorithmus der Reihe nach durchzählen? Im Blog Ein Algorithmus zum Aufzählen von Ausdrucksbäumen habe ich das Problem und die Lösung skizziert. In allen Details dargelegt, führte das zu einem kleinen Artikel, in dem ich den Weg vom Ausdrucksbaum über ein zugeordnetes Schema, eine Signatur und Zahlenfolgen, die ich "gerichtet" nannte, beschritt. In diesem Artikel beschreibe ich auch den gesuchten Algorithmus.

Das Produkt der Klammerungen mit den möglichen Belegungen der acht Operationen ergibt also 93'716'480 zu prüfende Ausdrücke - eine vom Rechner leicht zu handhabende Grösse, falls man sich für die richtige Programmiersprache entscheidet.[2]

Hier noch einmal die Anforderungen:

  • Es sollte nicht mit den üblichen Gleitkommazahlen, sondern mit rationalen Zahlen, also Paaren {zaehler,nenner} ganzer Zahlen gerechnet werden, damit nicht aufgrund von Genauigkeitsverlusten eine Lösung verlorengeht.

  • Die Auswertung eines Ausdrucks kann natürlich abgebrochen werden, sobald eine Division durch Null auftritt.

  • Nach dem Assoziativgesetz identische Ausdrücke wie a·(b·c) und (a·b)·c sollen nur einmal in der Ergebnisliste vorkommen.

  • Die gefundenen Klammerausdrücke sollten in einer für Menschen lesbaren Form ausgegeben werden – unabhängig davon, welche Datenstruktur intern zur Darstellung der Ausdrücke gewählt wird.

In der Produktbildung
Klammerungen × Mögliche Wahlen der Operationen
ist bereits ein ideales Schema für die Parallelisierung angelegt: Für jedes festes Klammerschema ergibt sich die Aufgabe, alle acht Plätze mit den verschiedenen möglichen Operationen zu belegen, den entstehenden Ausdruck zu evaluieren und zu prüfen, ob er das gewünschte Ergebnis liefert. Der Zeitpunkt, wenn neue Klammerschemata aufgebaut wurden, ist daher gut für einen Schnitt geeignet, bei dem man Threads starten und parallel abarbeiten kann.

Den vollständigen Quelltext und die Binaries des Programms, das ich zur Lösung dieser Aufgabe implementiert habe, findet man bei github unter dem Link

https://github.com/rplantiko/computeExpressions

Hier will ich einige Ideen der Implementierung genauer erörtern.

Eine grundlegende Designentscheidung war, für ein einzelnes Klammerschema eine Klasse ExpressionTemplate einzuführen. Ein Objekt des Type ExpressionTemplate enthält die Angaben über ein konkretes Klammerschema, aber auch einen Vektor für die gefundenen Ergebnisse. Es kann also mit seinen eigenen Daten arbeiten, ohne vom Rest des Programms abzuhängen.

Immer wenn der Algorithmus für Klammerungen ein neues Klammerschema aufgebaut hat, wird eine neue Instanz von ExpressionTemplate gebildet und kommt in einen Vorrat, die globale Variable
std::vector<ExpressionTemplate*> worklist;
Dies ist übrigens die einzige globale Variable, die sich während des Programmlaufs noch ändert.

Über einen Parameter numThreads, den man dem Programm zur Ausführung mitgeben kann, wird gesteuert, wieviele Klammerschemata von der CPU auf einmal verarbeitet werden dürfen. Sobald die worklist diese Anzahl von Einträgen erreicht hat, werden die Threads gestartet. Jeder Thread arbeitet auf einer eigenen Instanz von ExpressionTemplate und kommt damit weder den anderen noch dem aufrufenden Thread ins Gehege. Da die Threads ungefähr alle die gleiche Verarbeitungszeit haben, genügt ein FIFO-Join, um die Beendigung aller Threads abzuwarten: Auf die Beendigung des zuerst gestarteten Threads wird als erstes gewartet, danach auf den zweiten usw.[3]

Hier die Besucherfunktion nextExecutionPlan, die für jedes Klammerschema einmal aufgerufen wird.[4]
void nextExecutionPlan( node* root ) {

if (root != 0) { // There must be a final cleanup call with root = 0
worklist.push_back(
new ExpressionTemplate(root,numbers,expectedResult,size)
);
}

if (worklist.size() >= numThreads || root == 0) {

std::vector<pthread_t> threads;

for (auto templ = worklist.begin(); templ != worklist.end(); ++templ) {
pthread_t t;
pthread_create( &t, 0, evaluate, *templ );
threads.push_back( t );
}

// Waiting for the end of the threads, in the order of their creation
// (not so bad, but could be improved, of course)
for (auto t = threads.begin(); t != threads.end(); ++t ) {
pthread_join( *t, 0 );
}

// The work of this package is done. Print the result & free the ressources
for (auto templ = worklist.begin(); templ != worklist.end(); ++templ) {
(*templ)->printResults();
delete *templ;
}

worklist.resize( 0 );

}

}

Nach Ausführung muss schliesslich noch ein allfällig verbliebener Rest von Threads abgearbeitet werden (wenn nämlich die Gesamtzahl der Klammerschemata nicht durch die gewählte numThreads teilbar ist). Das erledige ich, indem ich am Schluss die Funktion noch einmal mit root = 0 aufrufe.

Wie habe ich die Auswertung eines Ausdrucks implementiert? Das Klammerschema wird im Konstruktor der Klasse ExpressionTemplate in eine "angereichertes" Datenstruktur überführt. Hier ist der Typ eines angereicherten Knotens:
typedef struct valuatedNode {
bool isLeaf;
int* op;
rational value;
struct valuatedNode* left;
struct valuatedNode* right;
} valuatedNode;

Neben den bekannten Elementen isLeaf, left und right gibt es einen Zeiger auf eine Ganzzahl in einem für die ganze Laufzeit der Klasse fixen Array ops = new int[size]; Die Plätze dieses Arrays nehmen nur die Zahlen 0 bis 3 an und werden als Offsets für die Ermittlung des richtigen Funktionszeigers aus der Konstanten OPERATION verwendet.
// Operations with rational numbers  
void plus(rational *first, rational second);
void minus(rational *first, rational second);
void times(rational *first, rational second);
void by(rational *first, rational second);

// Any operation on rationals
typedef void(*operation)(rational*, rational);

const operation OPERATION[4] = { plus, minus, times, by };

So braucht man für jede neue Auswahl von Operationen nur die Ganzzahlen im Array ops zu ändern, und die Evaluierung greift dann automatisch auf die richtige Operation zu.

Um die Kombinationen der Operationen der Reihe nach durchzugehen, verwende ich eine Ganzzahl, die ich einfach von 0 bis 65535 durchzähle. Ihre acht Bitpaare, die ich in den Array ops schreibe, liefern dann alle möglichen Variationen der Operationen. Ein weiterer Trick ist, dass ich die Stelle value des valuatedNode-Typs, die ja streng nur für Blätter, nicht für Operationsknoten nötig ist, für Operationsknoten mit dem aktuellen Zwischenergebnis besetze. Da ich den Array der Knoten von hinten nach vorne durcharbeite, sind alle Referenzen, auf die ein Knoten verweist, bereits mit Werten belegt. So kann der Ausdruck ohne Rekursion in einer ganz gewöhnlichen Schleife abgearbeitet werden.[5]
void ExpressionTemplate::evaluateAllocations() {

// Create all 4**size combinations of operations,
// using a bit field of length 2*size
// If size becomes > 31, think of another implementation of this part
long long int counter = 1 << 2*size;

while (counter) {
counter--;

long long int allocation = counter;
nodes[0].value.denom = 0; // Makes result unequal to everything
int iOp = size-1;

for (valuatedNode* current = nodes + 2*size;
current >= nodes;
current-- ) {

if (!current->isLeaf) {

// Pop next operation from bit field
ops[iOp--] = allocation & 0b11;
allocation >>= 2;

// Skip [ *, [ *, ...], <>[*,...] ] & the like by associativity
if (isDuplicateByAssociativity( current )) break;

// Compute next intermediate result
current->value = current->left->value;
OPERATION[*current->op](&current->value,current->right->value);

// ignore div 0
if (current->value.denom == 0) break;

}

}

if (equals(&expectedResult,nodes[0].value)) {
int* result = new int[size];
for (int i=0;i<size;i++) result[i] = ops[i];
results.push_back( result );
}
}
}


Die Laufzeiten des Programms ergaben noch eine Überraschung. Bei leichten Schwankungen von +/- 0.5 sec ergibt sich folgende Aufstellung [6].

Anz. Threads    Laufzeit
1 : 6.1 sec
2 : 4.1 sec
3 : 4.0 sec
5 : 3.3 sec
10 : 2.5 sec
20 : 2.5 sec
50 : 2.0 sec
100 : 1.9 sec
200 : 1.6 sec
500 : 1.6 sec
1000 : 1.5 sec
1500 : 1.4 sec

Naiv könnte man ja meinen, dass es nichts bringt, eine CPU, die acht Prozessorthreads hat, mit mehr als acht Threads zu behelligen. Die interne Pufferung von auszuführenden Threads ist aber auf modernen CPUs mittlerweile so weit optimiert, dass man offenbar am besten fährt, gleich alles was man hat, auf die CPU zu werfen und ihr selbst die Ausführung zu überlassen! Die Ausführungszeit sinkt kontinuierlich mit der Anzahl der gleichzeitig gestarteten Threads, ohne dass hier ein Minimum erkennbar wird! Wenn alle 1430 Threads gleichzeitig gestartet werden, ist man um den Faktor vier schneller als bei der serieller Abarbeitung.

Ach ja, die Ergebnisse selbst sollen nicht verschwiegen werden! Das Programm gibt die folgenden 41 Lösungen aus. Viele Mehrfachnennungen aufgrund des Assoziativgesetzes werden durch die Implementierung bereits unterdrückt. Manche der Lösungen würde man trotzdem noch als gleichartig betrachten.

2011 = 1-(2·(3+(((4-5)-(6+7))·8·9)))
2011 = 1-((2·3)+(4:((5-6):7)·8·9))
2011 = 1-((2·3)+(4·(5-6)·7·8·9))
2011 = 1+(2:(3:(((4+5)·6·7·8)-9)))
2011 = 1-((2·3)+(4:(((5-6):7):(8·9))))
2011 = 1-((2·3)+(4·(5-6)·7·8·9))
2011 = ((1+(2·(3+(4·5))·6))·7)+(8·9)
2011 = 1-(2·(3+((4-(5+6+7))·8·9)))
2011 = (((1:2)-((3-(4+5))·6·7))·8)-9
2011 = 1-(((2:3)+(4:(5-6)·7·8))·9)
2011 = 1-(((2:3)+(4·(5-6)·7·8))·9)
2011 = 1-(((2:3)+(4:((5-6):(7·8))))·9)
2011 = 1-(((2:3)+(4·(5-6)·7·8))·9)
2011 = 1+(2·3·((4·((5·6)+(7·8)))-9))
2011 = 1-(2:(3:(4+5-(6·7·8·9))))
2011 = (1·2)+((3+4)·((((5·6)+7)·8)-9))
2011 = 1·(2+((3+4)·((((5·6)+7)·8)-9)))
2011 = (1·2)+((3+4)·(5+(6·((7·8)-9))))
2011 = 1·(2+((3+4)·(5+(6·((7·8)-9)))))
2011 = 1+((2+((3-(4·((5:6)-7)))·8))·9)
2011 = 1-((2·3)+(4:(5-6)·7·8·9))
2011 = 1-((2·3)+(4·(5-6)·7·8·9))
2011 = 1+(2:(3:(((4+5)·6·7·8)-9)))
2011 = 1-((2·3)+(4:((5-6):(7·8·9))))
2011 = 1-((2·3)+(4·(5-6)·7·8·9))
2011 = ((1+2+3+4)·((5·6·7)-8))-9
2011 = 1·((((2·3)+4)·((5·6·7)-8))-9)
2011 = (1·((2·3)+4)·((5·6·7)-8))-9
2011 = (((1·2·3)+4)·((5·6·7)-8))-9
2011 = ((1+2+3+4)·((5·6·7)-8))-9
2011 = (1·2)-((3+4)·(((5-(6·7))·8)+9))
2011 = 1·(2-((3+4)·(((5-(6·7))·8)+9)))
2011 = 1-(((2:3)-4)·(5+6+(7·8))·9)
2011 = 1-(2·3·(((4-(5+(6·7)))·8)+9))
2011 = 1+(2:3·(4-(5-(6·7·8)))·9)
2011 = 1+(2:(3:((4-(5-(6·7·8)))·9)))
2011 = (1:(2-3))+(4·(5-(6-(7·8·9))))
2011 = (1·(2-3))+(4·(5-(6-(7·8·9))))
2011 = (1·2)-(3-(4·(5-(6-(7·8·9)))))
2011 = 1·(2-(3-(4·(5-(6-(7·8·9))))))
2011 = 1-(2:(3:(4+5-(6·7·8·9))))


Fussnoten



[1] Mein Kollege Rolf Bertschinger fand die folgende Lösung: Von 2011 eins abgezogen, ergibt 2010, eine durch 6 = 2*3 teilbare Zahl. Man hat also eventuell gute Chancen, wenn man den Ausdruck mit
1 + 2·3·...

beginnen lässt. 2010 ergibt, durch 6 dividiert, 335. Der gepunktete Teil des Ausdrucks sollte also irgendwie 335 ergeben. Nun muss man die Grösse der verbleibenden Möglichkeiten, Produkte zu bilden, ins Spiel bringen: Wir haben noch die Zahlen 4 bis 9 zur Verfügung. Wenn man Produkte von drei Zahlen bildet, ergibt 6·7·8 schon fast das Gewünschte, nämlich 336. Wenn wir die 9 einen Moment aussen vor lassen, hätten wir mit
4 - 5 + 6·7·8 = 335

schon einen brauchbaren Ansatz. Leider stört dann aber die 9. Es scheint hoffnungslos, die 9 nun noch wegzubekommen, wenn doch alle anderen Zahlen schon verbraucht sind! Das Aha ist nun: Wenn man vorne, statt, mit 3 zu multiplizieren, durch 3 dividiert, hat man im Vergleich den nachfolgenden Ausdruck durch 9 dividiert. Das wird gerade korrigiert, wenn man am Schluss noch mit 9 malnimmt! So kommt man auf
1 + ( 2 : 3 ) · ( 4 - 5 + 6·7·8) · 9


[2] Weil ich bei dieser Gelegenheit die Groovy-Bibliothek GPars kennenlernen wollte, machte ich eine erste Implementierung in Groovy. Das war für dieses Problem eine schlechte Wahl. Da Groovy eine hoch dynamische Sprache ist, die alles Nötige erst zur Laufzeit ermittelt, benötigte mein Programm, obwohl es parallelisiert auf 8 Prozessorthreads lief, zehn Stunden, bis es alle Ergebnisse ermittelt hatte! Das ist keine Kritik, sondern nur eine Beobachtung! Groovy ist eine extrem ausdrucksstarke Sprache, die sich dem Ideal "Development by Config" schon sehr stark nähert. Applikationsdesign, bei dem man es mit einem User Interface zu tun hat, kann bei Einsatz von Groovy nur gewinnen – wie etwa das Webframework Grails, oder das auf Swing aufsetzende MVC-Framework Griffon. Aber für rechenintensive Aufgaben wie diese ist Groovy eher nicht zu wählen.

Wen es interessiert: Meine Groovy-Implementierung für diese Aufgabe findet man bei pastebin unter http://pastebin.com/GzcSW0d9.

[3] Wenn die Ausführungszeiten der einzelnen Threads stark streuen, wäre hier ein anderes Verfahren angebracht. Man könnte eine synchronisierte Variable int threadsWorking einführen, die vom Thread beim Start hoch- und bei Beendigung heruntergezählt wird. Der Master-Thread wartet dann darauf, dass threadsWorking == 0 wird.

[4] Ich verwende aus technischen Gründen POSIX-Threads, obwohl ich eigentlich gerne einmal das neu in C++11 eingebaute Multithreading ausprobiert hätte. Leider ist es auf meinem Windows-System noch nicht verfügbar, da die MinGW es noch nicht enthält, und ich hatte bislang noch keine Lust, mir eine zweite Festplatte mit Linux einzurichten. (Alles, was mit blosser Computeradministration zu tun hat, ist unproduktiver, lästiger Zeitfrass.)

[5] Die Rekursion ist für Baumstrukturen im allgemeinen vorzuziehen, da sie lesbarer ist. An dieser Stelle aber, wo es auf höchstmögliche Effizienz ankommt, ist schon der Aufbau neuer Stackebenen pro Operation ein zu vermeidender Overhead.

[6] Ausgeführt auf meinem Rechner, der einen Prozessor vom Typ Intel Core i7 2600K @3.4 GHz hat und somit über acht Prozessorthreads verfügt.

Keine Kommentare :