Sinus-Eingabe skalieren

05/12/2009 - 10:49 von Markus Wichmann | Report spam
Hi all,
ich versuche hier gerade, den Sinus anzunàhern (wird mal ein Algorithmus
für die dietlibc). OK, bisher hab ich schon interne Abhàngigkeiten zu
floor(), ceil() und fabs(), was ja eigentlich zu vermeiden ist, aber
naja.

So, der Algorithmus, der den Sinus annàhert funktioniert nur mit einem
Eingabewert im Bereich von -Pi bis Pi. Alle anderen Eingaben müssen
entsprechend skaliert werden. Nun hab ich mir folgendes gedacht:

Sei x > 0, dann ist

sin x = sin (x - 2*n*Pi)

mit n = max {n \in N | 2*n*Pi < x}
= max {n \in N| n < x/(2*Pi)}

Das skaliert ein positives x auf den Bereich zwischen 0 und 2*pi. Nun
gut, sollte x immer noch größer als Pi sein, ziehe ich wieder 2*pi ab
und das Problem ist gelöst. So hab ich das ganze auch implementiert:

if (x > 2*M_PI)
x -= 2*M_PI*floor(x / 2*M_PI);

if (x > M_PI)
x -= 2*M_PI;

Die Definition von M_PI hab ich aus der math.h kopiert. Außerdem linke
ich das ganze (momentan) gegen die (g)libm.

Ach ja, für negative x dachte ich mir:

sin x = sin (x + 2*n*Pi)

mit n = min {n \in N | 2*n*Pi > -x}
= min {n \in N | n > -x/(2*Pi)}

Nur sagt mir mein Testprogramm irgendwie, dass das nicht klappt:

sin(6.00) = -0.278463
sin(6.25) = -0.032793
sin(6.50) = 332740.314793

Hier mal der komplette Algorithmus:

#define M_PI 3.14159265358979323846 /* pi */
double ceil(double);
double floor(double);
double fabs(double);
static const double A = -4. / (M_PI * M_PI);
static const double B = 4. / M_PI;
static const double P = 0.775;
static const double Q = 0.225;

double sin(double x)
{
if (x < -2*M_PI)
x += 2*M_PI*ceil(-x / 2*M_PI);

if (x > 2*M_PI)
x -= 2*M_PI*floor(x / 2*M_PI);

if (x > M_PI)
x -= 2*M_PI;

double y = A * x * fabs(x) + B * x;
return P * y + Q * y * fabs(y);
}

#ifdef WITH_MAIN
#include <stdio.h>
int main()
{
for (double d = -2; d < 3; d += 0.0625)
printf("sin(%f * Pi) = %f", d, sin(d*M_PI));
return 0;
}
#endif

Tschö,
Markus

P.S.: Mein Compiler scheint buggy zu sein: Wenn ich die ersten zwei
Zeilen des Sinus durch

if (x < -M_PI)
return -sin(-x);

ersetze, kommt eine Segmentation fault. Ich habe das nachgeforscht: Es
wird an dieser Stelle immer wieder sin(x) aufgerufen, nicht sin(-x).
WTF?

Ach ja: Angaben zu meinem System:
$ uname -a
Linux debbox 2.6.30-2-amd64 #1 SMP Fri Sep 25 22:16:56 UTC 2009 x86_64
GNU/Linux
$ gcc --version |& head -1
gcc (Debian 4.3.4-6) 4.3.4

Nur weil ein Genie nix reißt, muß ja nun nicht gleich jeder Idiot
pausieren... Bully hats ja auch geschafft.
 

Lesen sie die antworten

#1 ram
05/12/2009 - 13:30 | Warnen spam
Markus Wichmann writes:
sin x = sin (x - 2*n*Pi)



Ohne daß ich Deine Frage damit direkt beantworten will,
ist sicher eine interessante Lektüre zu diesem Thema:

http://blogs.sun.com/jag/entry/tran...meditation

Ähnliche fragen