Generische Sinus und Cosinus Implementierung

Normalerweise verlässt man sich auf die C Standard Bibliothek, wenn man etwas Trigonometrisches braucht.

Aber was tut man, wenn man diese nicht parat hat, wie in einer nativen EFI App, und auch kein Inline-Assembler vom Compiler unterstützt wird.

Oft ist es einfach wichtig, den richtigen Begriff zu kennen und den Rest erledigt die Suchmaschine des Vertrauens.
Und da war doch was in der Mathe-Vorlesung zur Analysis:

Die Taylor-Reihe

Wenn wir von einer Funktion immer Ableitungen bilden und diese nach dem Taylorschein Lehrsatz aufsummieren, erhält man abstrakt:

f(x) = f(x0) 
     + f1(x0) * (x - x0) ^ 1 / 1! 
     + f2(x0) * (x - x0) ^ 2 / 2!
     + ...
     + fn(x0) * (x - x0) ^ n / n!

Und nachdem die Ableitungen von Sinus und Kosinus sich mit alternierendem Vorzeichen abwechseln (cos, -sin, -cos, sin, cos, …) und sin(0) == 0 und cos(0) == 1 ist, haben wir in Schule folgende alternierende Annäherung kennen gelernt:

sin(x) = 0 + x^1 / 1! - x^3 / 3! + x^5 / 5! - x^7 / 7! + ...
cos(x) = 1 - x^2 / 2! + x^4 / 4! - x^6 / 6! + x^8 / 8! - ...

Zumindest, wenn x ein Winkel in Radianten zwischen -PI und +PI ist.

Die Taylor-Reihe wird um so genauer, je länger sie ist. Mathematiker können zeigen, dass wir nur ein paar Iterationen brauchen, um ausreichend genau zu sein, damit wir damit ein bisschen rechnen können. Siehe dotancohen.com/eng/taylor-sine.php.

Nachdem ich bei meinen Unit-Tests Genauigkeiten von mindestens 4 Nachkommastellen haben wollte, musste ich die Iterationen bis zu x^13 erweitern.

Sinus und Cosinus selbst implementieren

Wenn in Embedded Plattformen was fehlt, dann fast alles. Und nachdem man sein eigenes malloc und free nachgebaut hat, folgen die mathematischen Funktionen gänzlich ohne math.h.
Aber, ein Code sagt mehr als 1000 Worte:

 static double const pi = 3.14159265358979323846;
 
 static double normalize(double angle)
 {
   while(angle <= pi) angle += pi;
   while(angle > pi) angle -= pi;
   return angle;
 }
 
static double const n2 = 2.0;
static double const n3 = 6.0;
static double const n4 = 24.0;
static double const n5 = 120.0;
static double const n6 = 720.0;
static double const n7 = 5040.0;
static double const n8 = 40320.0;
static double const n9 = 362880.0;
static double const n10 = 3628800.0;
static double const n11 = 39916800.0;
static double const n12 = 479001600.0;
static double const n13 = 6227020800.0;

double approx_cos(double value)
{
  double const x = normalize(value);
  double const x3 = x * x * x;
  double const x5 = x3 * x * x;
  double const x7 = x5 * x * x;
  double const x9 = x7 * x * x;
  double const x11 = x9 * x * x;
  double const x13 = x11 * x * x;
  return x - x3 / n3 + x5 / n5 - x7 / n7 
       + x9 / n9 - x11 / n11 + x13 / n13;
}

double approx_sin(double value)
{
  double const x = normalize(value);
  double const x2 = x * x;
  double const x4 = x2 * x * x;
  double const x6 = x4 * x * x;
  double const x8 = x6 * x * x;
  double const x10 = x8 * x * x;
  double const x12 = x10 * x * x;
  return 1 - x2 / n2 + x4 / n4 - x6 / n6 
       + x8 / n8 - x10 / n10 + x12 / n12;
}

double approx_tan(double value)
{
  return approx_sin(value) / approx_cos(value);
}

So in etwa sieht der Kompatibilitäts-Layer im GATE Projekt aus, der genau für solche Spezialfälle aktiviert werden kann. (z.B. bei EFI Apps.)

Für den Arkussinus gibt es ebenfalls eine Reihenentwicklung und der Arkuskosinus sollte mit acos(x) = PI / 2.0 - asin(x) erreicht werden können.

Wenn nötig, kann man dann auch die “echten” C Funktionen so nachbauen, und genau das musste ich tun, um diverse Fremdbibliotheken durch den Linker zu bekommen.

Fazit

Wenn man “normale” Bibliotheken in diverse embedded Umgebungen einpflanzen möchte, passiert es, dass es keine C-Lib Implementierung gibt. Und dann tauchen zahlreiche “undefinierte Symbole” beim Linken auf.
Und das kann auch sin(), cos() und ihre Geschwister treffen.

Genau dafür ist es gut, dass man eine primitive Implementierung parat hat, damit man fortsetzen kann.
Selbstverständlich sollte immer eine ordentliche Mathematik-Bibliothek im finalen Kompilat landen, die auf die Plattform optimiert ist.

… doch leider ist das eben nicht immer möglich.