Auf der Suche nach einem Arcsin-Algorithmus

7

Hat jemand einen einfachen Algorithmus zur Berechnung eines einigermaßen genauen Arkussinus? Mit "einfach" meine ich eine Art Polynom, das <= 5 Multiplikationen pro Ausgangsabtastung erfordert. Und mit "ziemlich genau" meine ich ein Algo, dessen Fehler nicht mehr als 10% beträgt, wenn das Eingabeargument nahe bei plus oder minus eins liegt. Ich suchte eine Weile im Internet, fand aber nichts sofort Nützliches.

Richard Lyons
quelle
Dies könnte einige Ideen geben stackoverflow.com/questions/5920467/…
geometrikal
Aber warum nicht einfach eine Nachschlagetabelle?
geometrisch
Ich denke an eine Implementierung, bei der der verfügbare Speicher schmerzlich begrenzt ist. Daher habe ich keine Lösung für Nachschlagetabellen in Betracht gezogen. Danke für deine Gedanken.
Richard Lyons
Erlauben Sie Quadratwurzeln? Aufgrund des Verhaltens der Funktion nahe (unendliche Steigung) funktioniert eine Polynomnäherung dort nicht gut. ±1
Yves Daoust
Was ist mit CORDIC, das nur wenige Additionen und Subtraktionen und keine Multiplikationen benötigt?
Mattatt

Antworten:

5

Hier ist nur eine Polynomversion :

arcsin(x)=x+12x33+1324x55+135246x77
function y = arcsin_test3(x)
    y = x.*(1+x.*x.*(1/6+ x.*x.*(3/(2*4*5) + x.*x.*((1*3*5)/(2*4*6*7)))))
endfunction

Das scheint fünf Multiplikationen (vorausgesetzt, Sie können das Ergebnis von speichern x.*x) und drei Additionen zu haben.

Und die scilabHandlung ist:

Geben Sie hier die Bildbeschreibung ein

Top ist scilab‚s asinvs diese, unten ist der Fehler zwischen den beiden.


Ursprüngliche Antwort

Die Quadratwurzel hier mag ein Ärger sein, aber ich dachte, ich würde es aufschreiben, weil es nach Spaß aussieht. :-)

Diese Seite schlägt vor:

ab Seite 81 des Handbuchs der mathematischen Funktionen von Milton Abramowitz und Irene Stegun:

arcsin(x)=π/21x(a0+a1x+a2x2+a3x3),
wobei
a0=1.5707288a1=0.2121144a2=0.0742610a3=0.0187293

Ich habe dies in implementiert scilabund es funktioniert OK, außer um . Nur die auf reflektieren , eine viel bessere Annäherung.x=10x11x0

Das obere Diagramm zeigt scilabdie asinFunktion gegen die obige Annäherung (rot gestrichelt) gegen meine Änderung in Grün.

Das untere Diagramm zeigt den Fehler für meine Änderung (das Zeichnen des Originals auf derselben Achse bedeutet, dass das Grün überall Null aussieht).

Geben Sie hier die Bildbeschreibung ein

// 25770
function y = arcsin_test(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-x).*(a0 + a1*x + a2.*x.*x + a3.*x.*x.*x)

endfunction

function y = arcsin_test2(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-xx).*(a0 + a1*xx + a2.*xx.*xx + a3.*xx.*xx.*xx)

    y = y.*sign(x); 
endfunction

x = [-1: .0100001 : 1];

clf
subplot(211)
plot(x,arcsin_test2(x),'g.');
plot(x,arcsin_test(x),'r:');
plot(x,asin(x))
subplot(212)
//plot(x,(arcsin_test(x) - asin(x)),'r:')
plot(x,(arcsin_test2(x) - asin(x)),'g.')
Peter K.
quelle
2
"Handbuch der mathematischen Funktionen" lieben dieses Buch
geometrisch
1
Oh schieß. Ich habe diese Wiki-Webseite "Inverse Triggerfunktionen" während meiner Websuche gesehen, aber ich habe nicht weit genug nach unten gescrollt, um das Material "Unendliche Serien" zu sehen! Schande über mich. Peter K., das ist noch einer, den ich dir schulde. (Mein ursprüngliches Problem bestand darin, die Leistung eines digitalen Differenzierers mit zentraler Differenz zu verbessern, was meines Erachtens durch Ausführen einer Arkussinusoperation erreicht werden kann.)
Richard Lyons
Ja, aber Rick, du kannst keine unendliche Serie machen. Wenn Sie es endlich machen wollen, sind die optimalen Koeffizienten nicht genau das, was Sie durch Abschneiden der unendlichen Reihen erhalten. Wenn Sie MATLAB haben (sie haben jetzt eine relativ günstige "Home Use" -Lizenz), kann ich Ihnen MATLAB-Code für den Remez-Austausch über die Funktion Ihres Herzenswunsches senden.
Robert Bristow-Johnson
1
UPDATE: Ich konnte mein ursprüngliches Problem (Erstellen eines sehr einfachen digitalen Unterscheidungsmerkmals, das die Leistung gegenüber einem zentralen Differenzierungsunterscheidungsmerkmal verbessert hat) lösen, ohne eine arcsin () -Funktion zu verwenden. Ich kann auf dsprelated.com einen Blog veröffentlichen, in dem meine Ergebnisse beschrieben werden. Ich danke allen hier für ihre Hilfe!
Richard Lyons
3

Ich habe eine ziemlich gute Implementierung von hier .arctan()

Ich denke, Sie können die Identität verwenden:

arcsin(x)=arctan(x1x2)

um zu bekommen was du willst.

robert bristow-johnson
quelle
Ihr Link zu verschiedenen Funktionen ist interessant, Robert. Nur zum Kichern habe ich versucht, die Funktion sqrt (1 + x) zu implementieren. Anstatt die korrekten Grenzwerte von 0 bis 4 zu verwenden, habe ich es vermasselt und 1 bis 5 verwendet. Natürlich habe ich ein falsches Ergebnis berechnet. Als ich jedoch mein "falsches" Ergebnis mit zwei multiplizierte, erhielt ich das richtige Ergebnis. Interessant, oder?
Richard Lyons
Rick, ich bin mir ziemlich sicher, dass die Funktionen "korrekt" (oder ziemlich genau) sind, da sie mit den angegebenen Grenzen von angegeben sind. Für ist es nur für . Wenn Sie also zwischen und sind, muss dort eine kleine Konstante ( ) gespeichert sein und Sie müssen den Unterschied zwischen einem geraden Exponenten von und einem ungeraden Exponenten von . xx1x224222
Robert Bristow-Johnson
Ich glaube auch, dass Ihre Funktionen korrekt sind. Ich habe lediglich einen dummen Fehler "Summationsgrenzen" kommentiert, den ich gemacht habe, und als ich meinen Fehler gemacht habe, habe ich ein falsches Ergebnis berechnet. Aber ich bemerkte, dass mein falsches Ergebnis genau die Hälfte des richtigen Ergebnisses war. Ich habe nur gesagt, dass mein falsches Ergebnis eine "interessante" Beziehung zum richtigen Ergebnis hat, das ist alles. Entschuldigung für die Verwirrung, Robert.
Richard Lyons
2

Der zentrale Teil der Kurve ist kein wirkliches Problem, da er ziemlich linear ist und die Taylor-Annäherung an zwei oder drei Terme ein guter Ausgangspunkt ist (Polynomanpassung der kleinsten Quadrate etwas besser).

Die Seiten sind wegen der unendlichen Neigung problematischer. Ein Weg, um damit umzugehen, ist über die Transformation

arcsin(x)=π2arcsin(1x2),

das beinhaltet eine Quadratwurzel.


Wenn Ihr Argument mit Gleitkomma dargestellt wird, erhalten Sie eine schnelle Approximation der Quadratwurzel, indem Sie den Exponenten halbieren und eine lineare Transformation auf die Mantisse anwenden.z

Sei mit , dann . Sie können durch .z=m2e1m<2z=m2e/2m(21)(m+2)

  • Nehmen Sie den Exponenten auseinander (wenn Sie ihn löschen, erhalten Sie die Darstellung von ).em
  • wenn ist, berechne ;e(21)(m+2)
  • Wenn ungerade ist, berechne ;e2(21)(m+2)
  • Stellen Sie den Exponenten .e/2

Geben Sie hier die Bildbeschreibung ein

Yves Daoust
quelle
Ich freue mich darauf, mit diesem interessanten Quadratwurzel-Algorithmus zu experimentieren!
Richard Lyons
Bei diesen Koeffizienten ist der Fehler immer positiv. Durch eine leichte Anpassung können wir es symmetrisch machen und den maximalen Fehler halbieren.
Yves Daoust