Sinus berechnung

Die laufenden Diskussionen um eine schnelle und kompakte Berechnung der Sinus- und Cosinusfunktion hat mich angeregt, die einfache Taylor-Entwicklung um Stützstellen in einem Quadranten effektiv zu programmieren. Aus dem Analysis-Lehrbuch bekommen wir:

sin(x_0+dx)=sin(x_0)+cos(x_0)*dx-1/2*sin(x_0)*dx^2-1/6*cos(x_0)*dx^3

Ich teile nun das Intervall [0,pi/2] in n Intervalle. An den (n-1) Teilintervallgrenzen berechne ich nur den Sinus mit z. B. 16 Bit Genauigkeit, so dass die (n-1) Werte je in 2 Byte abgespeichert werden. Im Folgenden verwende ich 16 Bit und n=5, so dass 4 16 bit Werte vorab berechnet werden müssen.

Für einen beliebigen Wert x im Intervall [0,1] oder äquivalent im Bereich xb=[0,pi/2] bestimme ich die nächstgelegene Stützstelle k und den Abstand von der Stützstelle DX bzw. dx.

Den Cosinus Wert bekomme ich aus der Sinustabelle durch Ersetzung k->n-k. Folgendes Python-Programm bestimmt für 1 Mio zufallswerte im Intervall [0,pi/2] den Sinus-wert und vergleicht mit dem der c-Bibliothek:

---

 #! /usr/bin/python
# Interpolationsformel fuer den sinus im Intervall 0-pi/2
#Zahl der nicht trivialen Stuetzstellen n-1
#Genauigkeit der abgespeicherten Stuetzwerte in bitlaenge=nbit
#es wird von rechts oder links von der Stützstelle mit dem Ausdruck
#sin(x0+dx)=sin(x0)*(1-dx**2/2)+cos(x0)*dx*(1-1/6*dx**2)
#genaehert.
#Die cosinuswerte stehen in der gleichen Tabelle rueckwaerts
#bei 4 nichttrivialen 16 bit Zahlen ist die Genauigkeit 3.3*10^-5
import math,random
nbit=16
div=2**nbit
n=5
sin=(n+1)*[0]
sin[0]=0
for i in range(n):
  msin=math.sin(math.pi/2.0/n*(i+1))
  sin[i+1]=int(msin*div)
  print i+1,msin,sin[i+1]
maxdif=0
def sinp(x):
  k=int(x*n+0.5)
  DX=x*n-k
  dx=DX*math.pi/2.0/n
  dx2=dx**2
  si=sin[k]*(1.0-0.5*dx2)
  co=sin[n-k]*dx*(1.0-1.0/6.0*dx2)
  return (si+co)/float(div)

for i in range(1000000):
  x=random.random()
  xb=x*math.pi/2.0
#inline version ersetzt Funktionsaufruf sinp(x)
#  k=int(x*n+0.5)
#  DX=x*n-k
#  dx=DX*math.pi/2.0/n
#  dx2=dx**2
#  si=sin[k]*(1.0-0.5*dx2)      
#  co=sin[n-k]*dx*(1.0-1.0/6.0*dx2)
#  isin=(si+co)/float(div)
  isin=sinp(x) #(si+co)/float(div)
  masin=math.sin(xb)
  if abs(isin-masin)>maxdif:
    maxdif=abs(isin-masin)
print "maxdif=",maxdif

---

Das Ergebnis ist mit n=5 und nbit=16:

1 0.309016994375 20251
2 0.587785252292 38521
3 0.809016994375 53019
4 0.951056516295 62328
5 1.0 65536
maxdif= 3.27405091449e-05