Schnelle 32Bit-Integer Sinusberechnung

Schnelle 32Bit-Integer Sinusberechnung

Kürzlich fand ich [hier] eine genial einfache Methode, wie man mit einer Tabelle aus 21 Werten und einem kurzen Programm eine genial einfache Sinusberechnung mit erstaunlicher Genauigkeit von 0,1% erzielen kann. Die gesamte Berechnung benötigt nur <200Bytes Code und Daten. Zum Vergleich: Wenn man die Floating point Arithmetik des GCC benutzt, ist man gleich mit 7,5 kByte dabei.

Mathematische Grundlagen

Wir nehmen uns das bekannte Sinus-Cosinus-Additionstheorem:

[math]\displaystyle{ \sin(a+b) = \sin(a) \cos(b) + \sin(b) \cos(a) }[/math]

Wir zerlegen unseren Winkel a+b in Zehnerstufen: a=n·10° + b=Rest

sin(a)	bestimmen wir durch eine überschaubare Tabelle mit 10 Einträgen
cos(a)	bestimmen wir mit derselben Tabelle
cos(b)	nähern wir durch eine überschaubare Tabelle mit 10 Einträgen
sin(b)	nähern wir durch eine Gerade.

Letzteres geht deshalb, weil für kleine x gilt (x in Bogenmaß):

[math]\displaystyle{ \sin(x) \approx x }[/math]

Berechnet wird im Programm natürlich nur der erste Quadrant, die anderen werden einfach daraus abgeleitet.

Programmcode

#define MAX16BIT    0x7FFF

// benoetigt fuer die Umrechnung Grad in Bogenmass bei der Annaeherung sin(b)=b
int hollyConstant = MAX16BIT*0.017453292519943295769236907684886; // (Pi/2)/90°

int sinTable[] = {
MAX16BIT*0.0,                                    //sin(0)
MAX16BIT*0.17364817766693034885171662676931 ,    //sin(10)
MAX16BIT*0.34202014332566873304409961468226 ,    //sin(20)
MAX16BIT*0.5 ,                                   //sin(30)
MAX16BIT*0.64278760968653932632264340990726 ,    //sin(40)
MAX16BIT*0.76604444311897803520239265055542 ,    //sin(50)
MAX16BIT*0.86602540378443864676372317075294 ,    //sin(60)
MAX16BIT*0.93969262078590838405410927732473 ,    //sin(70)
MAX16BIT*0.98480775301220805936674302458952 ,    //sin(80)
MAX16BIT*1.0                                     //sin(90)
};

int cosTable[] = {
MAX16BIT*1.0 ,                                   //cos(0)
MAX16BIT*0.99984769515639123915701155881391 ,    //cos(1)
MAX16BIT*0.99939082701909573000624344004393 ,    //cos(2)
MAX16BIT*0.99862953475457387378449205843944 ,    //cos(3)
MAX16BIT*0.99756405025982424761316268064426 ,    //cos(4)
MAX16BIT*0.99619469809174553229501040247389 ,    //cos(5)
MAX16BIT*0.99452189536827333692269194498057 ,    //cos(6)
MAX16BIT*0.99254615164132203498006158933058 ,    //cos(7)
MAX16BIT*0.99026806874157031508377486734485 ,    //cos(8)
MAX16BIT*0.98768834059513772619004024769344      //cos(9)
};

/* Integer Sinus-Funktion
--------------------------------------
Prototype: int Sinus ( int angle );

Example: i = Sinus (30);
Result:	 INT32_MAX/2 * 0.5    (sin(30°)=0,5)    
*/
int Sinus ( int angle ){

   int a, b, quadrant, ret;

   quadrant = (angle % 360)/90;
   angle = (angle % 90);     // modulo 90
   if ((quadrant%2)!=0) angle = (90 - angle);

   a = angle / 10;
   b = angle - 10 * a;
   ret = sinTable[a] * cosTable[b] + b * hollyConstant * sinTable[9-a];

   if (quadrant>=2) ret=-ret;
   return ret;
}

Erweiterter Source code

Der gleiche Code mit gleicher Genauigkeit funtioniert natürlich auch für Winkel mit 1-2 Nachkommastellen:

Source code: [Sinus.zip]

Quellen und Links

• Fast sine function
• Sinusberechnung für 32Bit-Integer
• Diskussion über die Taylor-Reihe, Beitrag von waldo
• Alternative Interpolation, Artikel: Sinus berechnung
• Genaue Iteration für avr-gcc, Artikel: Sinus und Cosinus (CORDIC)
• Lineare Interpolation für avr-gcc: Sinus und Cosinus (Lineare Interpolation)

Autor: Marten Petschke