GNU libc mit Fehlern in Funktionen (sin etc.)

Stefan Kanthak hat mich vor einigen Tagen auf unschöne Fehler in der GNU libc aufmerksam gemacht. Trigonometrische Berechnungen führen zu fehlerhaften Werten, die weit oberhalb dessen liegen, was Mathematiker als Abweichungen gefunden und bestätigt haben.

In einer Mail schrieb Kanthak: "Auf Unix/POSIX-Systemen stellt die 'libc' die von C-Programmen verwendeten Standardfunktionen sowie die Schnittstelle zwischen Benutzerprogrammen und Kernel bereit. Für mathematische (trigonometrische, transzendentale, etc.) Funktionen gibt es noch die 'libm'. Auf GNU/Linux bzw. GNU/*-Systemen stellt die GNU libc beides bereit." Laut Kanthak arbeitet nicht nur das CERN mit der GNU libm und vertraut den Ergebnissen.

Die Abweichungen bei den Ergebnissen

Es gibt eine Dokumentation 19.7 Errors in Math Functions, die sich mit den Fehlern in den mathematischen Funktionen befassen. Dort heißt es u.a.

Errors are measured in "units of the last place". […]
…
Each function with a floating-point result behaves as if it computes an infinite-precision result that is within a few ulp of the mathematically correct value of the function […]
…
the math testsuite only flags results larger than 9ulp (…) as errors; although most of the implementations show errors smaller than the limit.

Kanthak wies darauf hin, dass im Februar 2026 ein paar bekannte Mathematiker (u.a. Dr. Brian Gladman, Prof. Dr. Paul Zimmermann, …), die seit über 30 Jahren mit Gleitpunkt-Arithmetik "spielen", eine Untersuchung zur (Un)genauigkeit verschiedener
"libm"-Implementierungen veröffentlicht. Dort bescheinigen sie sie (fast) allen
mathematischen Funktionen der GNU libm einen maximalen Fehler von 0.5 ULP bescheinigen.

Es gibt Diskrepanzen bei den Ergebnissen

Kanthak schrieb mir, dass er diese obige Untersuchung nicht kannte. Bei eigenen "Spielereien" (Implementation von Payne/Hanek-Argumentreduktion) stieß er binnen Minuten auf das Problem, dass die Ergebnisse größer Abweichungsfehler zeigen.

• cos(214112296674652) weist einen Fehler von 179 ULP auf
• sin(1284673780047912) sowie
• tan(1284673780047912) weisen jeweils einen Fehler von 268 ULP auf.

Hätte noch ein Fehler von Kanthak sein können. Daher hat Kanthak bei <libc-maintainers@gnu.org> nachfragt, ob 179 ULP oder 268 ULP als "a few ulp" gelten. Darauf wird der französische Mathematiker Paul Zimmermann aufmerksam. Wenn ich es richtig gesehen habe, hat Zimmermann einiges zu Abweichungsfehlern in solchen Bibliotheken geschrieben. Im aktuellen Fall bestätigt Zimmermann die von Kanthak ermittelten "groben" Fehler in den trigonometrischen Funktionen.

Weitere Tests von Kanthak

Kanthak hat dann die natürlichen Zahlen 214112296674652 und 1284673780047912 als "fast Vielfaches" von pi, pi/2 oder pi/4 als sogenannte Konvergenten von pi verwendet. Der Zähler/Numerator des Bruchs 6134899525417045/1952799169684491, der 105 Bits bzw. 32 Dezimalstellen von pi liefert, ist ein solcher Konvergent. Zähler/Numerator und Nenner/Denumerator sind kleiner 2**53 und exakt als (doppelt genaue) Gleitpunktzahlen darstellbar.

Der Bruch 5920787228742393/3769290217798865 liefert 32 Dezimalstellen
von pi/2, der Bruch 6027843377079719/7674888557167847 liefert 32 Dezimalstellen
von pi/4.

Kanthak hat (der Einfachheit halber) nur Werte kleiner 2**53 = 0x1.0p+53 =
9007199254740992 verwendet, die als Gleitpunktzahlen EXAKT darstellbar
sind. Es gibt diese auch größer als 2**53. Dabei laufen die ULPs "weg, z.B. liefert:

• sin(0x1.4c96c11134d36p+578) einen Fehler von 102825 ULP,
• cos(0x1.69eab0985179bp+246) einen Fehler von 136921 ULP,
• tan(0x1.69eab0985179bp+246) einen Fehler von 143270 ULP

Das ist also mehr als das 10000-Fache dessen, was GNU als fehlerhaft betrachtet, bzw. das ist das 286000-Fache dessen, was die oben erwähnten Mathematiker als "größten gefundenen Fehler" in der Bibliothek angeben. Klingt nicht gut.

Nachtrag: Die libc-Maintainer haben Kanthak inzwischen mitgeteilt, dass sie diesen seit 2013 vorhandenen Fehler (in einer von IBM eingereichten Routine) in der im Oktober 2026 kommenden Version 2.44 behoben haben (siehe Bug 15267 und Bug 34376).