Ich muss Rotationstransformationen (und andere) in einer eingebetteten Anwendung durchführen, für die die Funktionen sin(), cos() und tan() erforderlich sind. Ich weiß, dass Sie Nachschlagetabellen verwenden können, und das ist die einzige Lösung, die ich bei meinen eigenen Recherchen finden konnte, aber gibt es da draußen eine gute Festkomma-Trigger-Bibliothek?
Ich denke darüber nach, einen Cortex M3 für die Anwendung zu verwenden, daher möchte ich mich so weit wie möglich von Gleitkommazahlen fernhalten, um Anwendungen flink zu halten.
Ein guter Ansatz für die Trigonometrie in eingebetteten Anwendungen besteht darin, polynomiale Annäherungen an die benötigten Funktionen zu verwenden. Der Code ist kompakt, die Daten bestehen aus wenigen Koeffizienten, und die einzigen erforderlichen Operationen sind Multiplizieren und Addieren/Subtrahieren. Viele eingebettete Systeme verfügen über Hardware-Multiplikatoren, die eine gute Leistung bieten.
_Fractist IMHO ein Stück Mist; Ich hasse die Tatsache, dass es vom C-Komitee "standardisiert" wurde. Es zwingt Sie, Q15 oder Q31 für alles zu verwenden, was in vielen Situationen keinen Sinn ergibt, und lässt Sie in diesen Situationen ohne Hilfe festsitzen.Sind Sie dagegen, die Fixpunkt-Cortex-Bibliotheken dafür zu verwenden?
q31_t arm_sin_q31 (q31_t x)
Schnelle Annäherung an die trigonometrische Sinusfunktion für Q31-Daten.
von:
CMSIS-DSP: DSP-Bibliothekssammlung mit über 60 Funktionen für verschiedene Datentypen: Fixpunkt (fraktional q7, q15, q31) und Gleitkomma mit einfacher Genauigkeit (32-Bit). Die Bibliothek ist für Cortex-M0, Cortex-M3 und Cortex-M4 verfügbar.
Es verwendet eine Nachschlagetabelle mit quadratischer Interpolation, aber es ist ziemlich schnell. Sie könnten es für eine schnellere Geschwindigkeit, aber mehr Fehler an die lineare Interpolation anpassen.
Beachten Sie auch, dass selbst der Cortex M4 nicht unbedingt über eine FPU verfügt. Ich habe gesehen, dass sie "M4F" genannt wurden, wenn sie es tun.
Diese Antwort soll die derzeit akzeptierte Antwort mit einem konkreten Beispiel in zwei Varianten ergänzen und einige spezifische Designhinweise geben.
Die Polynomnäherung ist oft ein überlegener Ansatz, wenn die gewünschte Genauigkeit ziemlich hoch ist und ein Hardware-Multiplikator verfügbar ist. Tabellengrößen neigen dazu, schnell zuzunehmen, selbst wenn Interpolations- (z. B. lineare, quadratische) und Komprimierungsschemata (z. B. zweiteilige Tabellen) verwendet werden, sobald mehr als etwa 16 "gute" Bits erforderlich sind.
Die Verwendung von Minimax-Approximationen für Polynome wird dringend empfohlen, da sie den maximalen Fehler über das Intervall, für das sie generiert werden, minimieren. Dies kann zu einer erheblichen Verringerung der Anzahl der für eine bestimmte Genauigkeit erforderlichen Terme führen, verglichen mit beispielsweise Taylor-Reihenentwicklungen, die nur an dem Punkt, um den sie entwickelt werden, die beste Genauigkeit liefern. Häufig verwendete Tools wie Mathematica, Maple und das Open-Source -Tool Sollya bieten integrierte Methoden zur Generierung von Minimax-Approximationen.
Multiply-high-Operationen sind der grundlegende Rechenbaustein der Polynomauswertung in der Festkommaarithmetik. Sie geben die signifikantere Hälfte des vollen Produkts einer ganzzahligen Multiplikation zurück. Die meisten Architekturen bieten Varianten mit und ohne Vorzeichen, andere bieten Multiplikationen mit Ergebnissen doppelter Breite, die in zwei Registern zurückgegeben werden. Einige Architekturen bieten sogar Multiplizieren-Hoch-Plus-Addieren-Kombinationen, die besonders nützlich sein können. Optimierende Compiler sind typischerweise in der Lage, HLL-Quellcode-Idiome (wie die im ISO-C-Code unten verwendeten) entsprechend diesen Operationen in die entsprechenden Hardwareanweisungen zu übersetzen.
Um die Genauigkeit der Polynomauswertung zu maximieren, möchte man während der Zwischenberechnung jederzeit die maximal mögliche Anzahl von Bits verwenden, indem man ein Festkommaformat mit der maximal möglichen Anzahl von Bruchbits wählt. Aus Effizienzgründen vermeidet ein Skalierungsfaktor, der gleich der Registerbreite ist, die Notwendigkeit, über Verschiebungen neu zu skalieren, wenn er in Verbindung mit den Multiplikations-Hoch-Operationen verwendet wird.
Während das Horner-Schema typischerweise bei Gleitkommaberechnungen verwendet wird, um Polynome mit hoher Genauigkeit auszuwerten, ist dies bei Festkommaberechnungen oft unnötig und kann aufgrund der langen Abhängigkeitskette der Polynomauswertung, die die mehrfache Latenz aufdeckt, die Leistung beeinträchtigen. Parallele Bewertungsschemata, die die beste Ausnutzung von Pipeline-Multiplikatoren mit Multi-Zyklus-Latenz ermöglichen, sind oft vorzuziehen. Im folgenden Code kombiniere ich die Terme jedes Polynoms paarweise und baue daraus die Auswertung des vollständigen Polynoms auf.
Der folgende ISO-C-Code demonstriert die gleichzeitige Berechnung von Sinus und Cosinus gemäß diesen Konstruktionsprinzipien auf dem Intervall [0, π/2], wobei Ein- und Ausgänge im Format S8.23 (Q8.23) vorliegen. Es erzielt im Wesentlichen vollständig genaue Ergebnisse mit einem maximalen Fehler in der Größenordnung von 10 –7 und mehr als 80 % der korrekt gerundeten Ergebnisse.
Die erste Variante, in sincos_fixed_nj(), verwendet einen klassischen Ansatz der Argumentreduktion in [0, π/4] und eine Polynomnäherung an Sinus und Cosinus in diesem Intervall. Die Rekonstruktionsphase bildet dann die Polynomwerte quadrantenbasiert auf Sinus und Cosinus ab. Die zweite Variante, in sincos_fixed_ollyw, basiert auf einem Blogbeitrag von OllyW . Sie schlagen vor, die Transformation a = (2/π)x-1/2 auf das Intervall [-1/2, 1/2] anzuwenden, auf dem man dann sin ((2πa + π)/4 und cos approximieren muss ((2πa + π)/4. Die Reihenentwicklungen dieser ( sin , cos ) sind identisch, außer dass das Vorzeichen für Terme mit ungerader Potenz umgekehrt wird. Das bedeutet, dass man die Terme mit ungerader Potenz und gerader Potenz separat summieren kann und dann Berechnen Sie Sinus und Cosinus als Summe und Differenz der akkumulierten Summen.
Mit Compiler Explorer habe ich mit Clang 11.0 für ein armv7-a32-Bit-ARM-Ziel mit vollständiger Optimierung ( -O3) kompiliert. Beide Varianten wurden in Subroutinen mit 41 Anweisungen kompiliert , wobei jede Subroutine neun gespeicherte 32-Bit-Konstanten verwendet. sincos_fixed_ollyw()verwendet einen Multiplikationsbefehl mehr als sincos_fixed_nj, hat aber einen etwas geringeren Registerdruck. Die Situation scheint ähnlich zu sein, wenn man mit Clang für andere Architekturziele baut, also würde man beide Varianten ausprobieren wollen, um zu sehen, welche auf einer bestimmten Plattform besser abschneidet. Der Tangens könnte berechnet werden, indem das Sinusergebnis durch das Kosinusergebnis dividiert wird.
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#define SINCOS_NJ (1)
#define SINCOS_OLLYW (2)
#define VARIANT (SINCOS_NJ)
/* a single instruction in many 32-bit architectures */
uint32_t umul32hi (uint32_t a, uint32_t b)
{
return (uint32_t)(((uint64_t)a * b) >> 32);
}
/* a single instruction in many 32-bit architectures */
int32_t mul32hi (int32_t a, int32_t b)
{
return (int32_t)(uint32_t)((uint64_t)((int64_t)a * b) >> 32);
}
/*
compute sine and cosine of argument in [0, PI/2]
input and output in S8.23 format
max err sine = 9.86237533e-8 max err cosine = 1.02729891e-7
rms err sine = 4.11141973e-8 rms err cosine = 4.11752018e-8
sin correctly rounded: 10961278 (83.19%)
cos correctly rounded: 11070113 (84.01%)
*/
void sincos_fixed_nj (int32_t x, int32_t *sine, int32_t *cosine)
{
// minimax polynomial approximation for sine on [0, PI/4]
const uint32_t s0 = (uint32_t)(1.9510998390614986e-4 * (1LL << 32) + 0.5);
const uint32_t s1 = (uint32_t)(8.3322080317884684e-3 * (1LL << 32) + 0.5);
const uint32_t s2 = (uint32_t)(1.6666648373939097e-1 * (1LL << 32) + 0.5);
const uint32_t s3 = (uint32_t)(9.9999991734512150e-1 * (1LL << 32) + 0.5);
// minimax polynomial approximation for cosine on [0, PI/4]
const uint32_t c0 = (uint32_t)(1.3578890357166529e-3 * (1LL << 32) + 0.5);
const uint32_t c1 = (uint32_t)(4.1654359549283981e-2 * (1LL << 32) + 0.5);
const uint32_t c2 = (uint32_t)(4.9999838648363948e-1 * (1LL << 32) + 0.5);
const uint32_t c3 = (uint32_t)(9.9999997159466147e-1 * (1LL << 32) + 0.5);
// auxilliary constants
const int32_t hpi_p23 = (int32_t)(3.141592653590 / 2 * (1LL << 23) + 0.5);
const int32_t qpi_p23 = (int32_t)(3.141592653590 / 4 * (1LL << 23) + 0.5);
const int32_t one_p23 = (int32_t)(1.0000000000000e+0 * (1LL << 23) + 0.5);
uint32_t a, s, q, h, l, t, sn, cs;
/* reduce range from [0, PI/2] to [0, PI/4] */
t = (x > qpi_p23) ? (hpi_p23 - x) : x; // S8.23
/* scale up argument for maximum precision in intermediate computation */
a = t << 9; // U0.32
/* pre-compute a**2 and a**4 */
s = umul32hi (a, a); // U0.32
q = umul32hi (s, s); // U0.32
/* approximate sine on [0, PI/4] */
h = s3 - umul32hi (s2, s); // U0.32
l = umul32hi (s1 - umul32hi (s0, s), q); // U0.32
sn = umul32hi (h + l, a); // U0.32
/* approximate cosine on [0, PI/4] */
h = c3 - umul32hi (c2, s); // U0.32
l = umul32hi (c1 - umul32hi (c0, s), q); // U0.32
cs = h + l; // U0.32
/* round results to target precision */
sn = ((sn + 256) >> 9); // S8.23
cs = ((cs + 256) >> 9); // S8.23
/* cosine result overflows U0.32 format for small arguments */
cs = (t < 0xb50) ? one_p23 : cs; // S8.23
/* map sine/cosine approximations based on quadrant */
*sine = (t != x) ? cs : sn; // S8.23
*cosine = (t != x) ? sn : cs; // S8.23
}
/*
compute sine and cosine of argument in [0, PI/2]
input and output in S8.23 format
max err sine = 1.13173883e-7 max err cosine = 1.13158773e-7
rms err sine = 4.30955921e-8 rms err cosine = 4.31472191e-8
sin correctly rounded: 10844170 (82.30%)
cos correctly rounded: 10855609 (82.38%)
Based on an approach by OllyW (http://www.olliw.eu/2014/fast-functions/,
retrieved 10/23/2020). We transform a = 2/PI*x-1/2, then we approximate
sin ((2*PI*a + PI)/4 and cos ((2*PI*a + PI)/4. Except for sign flipping
in the odd-power terms of the expansions the two series expansions match:
https://www.wolframalpha.com/input/?i=series++sin+%28%282*pi*a+%2B+pi%29%2F4%29
https://www.wolframalpha.com/input/?i=series++cos+%28%282*pi*a+%2B+pi%29%2F4%29
This means we can sum the odd-power and the even-power terms seperately,
then compute the sum and difference of those sums giving sine and cosine.
*/
void sincos_fixed_ollyw (int32_t x, int32_t *sine, int32_t *cosine)
{
// minimax polynomial approximation for sin ((2*PI*a + PI)/4 on [-0.5, 0.5]
const uint32_t c0 = (uint32_t)(7.0710676768794656e-1 * (1LL << 32) + 0.5);
const uint32_t c1 = (uint32_t)((1.110721191857 -.25) * (1LL << 32) + 0.5);
const uint32_t c2 = (uint32_t)(8.7235601339489222e-1 * (1LL << 32) + 0.5);
const uint32_t c3 = (uint32_t)(4.5677902549505234e-1 * (1LL << 32) + 0.5);
const uint32_t c4 = (uint32_t)(1.7932640877552330e-1 * (1LL << 32) + 0.5);
const uint32_t c5 = (uint32_t)(5.6449491763487458e-2 * (1LL << 32) + 0.5);
const uint32_t c6 = (uint32_t)(1.4444266213104129e-2 * (1LL << 32) + 0.5);
const uint32_t c7 = (uint32_t)(3.4931597765535116e-3 * (1LL << 32) + 0.5);
// auxiliary constants
const uint32_t twoopi = (uint32_t)(2/3.1415926535898 * (1LL << 32) + 0.5);
const uint32_t half_p31 = (uint32_t)(0.5000000000000 * (1LL << 31) + 0.5);
const uint32_t quarter_p30 = (uint32_t)(0.2500000000 * (1LL << 30) + 0.5);
uint32_t s, t, q, h, l;
int32_t a, o, e, sn, cs;
/* scale up argument for maximum precision in intermediate computation */
t = (uint32_t)x << 8; // U1.31
/* a = 2*PI*x - 0.5 */
a = umul32hi (twoopi, t) - half_p31; // S0.31
/* precompute a**2 and a**4 */
s = (uint32_t)mul32hi (a, a) << 2; // U0.32
q = umul32hi (s, s); // U0.32
/* sum odd power terms; add in second portion of c1 (= 0.25) at the end */
h = c1 - umul32hi (c3, s); // U0.32
l = umul32hi ((c5 - umul32hi (c7, s)), q); // U0.32
o = ((h + l) >> 2) + quarter_p30; // S1.30
o = mul32hi (o, a); // S2.29
/* sum even power terms */
h = c0 - umul32hi (c2, s); // U0.32
l = umul32hi ((c4 - umul32hi (c6, s)), q); // U0.32
e = (h + l) >> 3; // S2.29
/* compute sine and cosine as sum and difference of odd / even terms */
sn = e + o; // S2.29 sum -> sine
cs = e - o; // S2.29 difference -> cosine
/* round results to target precision */
sn = (sn + 32) >> 6; // S8.23
cs = (cs + 32) >> 6; // S8.23
*sine = sn;
*cosine = cs;
}
double s8p23_to_double (int32_t a)
{
return (double)a / (1LL << 23);
}
int32_t double_to_s8p23 (double a)
{
return (int32_t)(a * (1LL << 23) + 0.5);
}
/* exhaustive test of S8.23 fixed-point sincos on [0,PI/2] */
int main (void)
{
double errc, errs, maxerrs, maxerrc, errsqs, errsqc;
int32_t arg, sin_correctly_rounded, cos_correctly_rounded;
#if VARIANT == SINCOS_OLLYW
printf ("S8.23 fixed-point sincos OllyW variant\n");
#elif VARIANT == SINCOS_NJ
printf ("S8.23 fixed-point sincos NJ variant\n");
#else // VARIANT
#error unsupported VARIANT
#endif // VARIANT
maxerrs = 0;
maxerrc = 0;
errsqs = 0;
errsqc = 0;
sin_correctly_rounded = 0;
cos_correctly_rounded = 0;
for (arg = 0; arg <= double_to_s8p23 (3.14159265358979 / 2); arg++) {
double argf, refs, refc;
int32_t sine, cosine, refsi, refci;
#if VARIANT == SINCOS_OLLYW
sincos_fixed_ollyw (arg, &sine, &cosine);
#elif VARIANT == SINCOS_NJ
sincos_fixed_nj (arg, &sine, &cosine);
#endif // VARIANT
argf = s8p23_to_double (arg);
refs = sin (argf);
refc = cos (argf);
refsi = double_to_s8p23 (refs);
refci = double_to_s8p23 (refc);
/* print function values near endpoints of interval */
if ((arg < 5) || (arg > 0xc90fd5)) {
printf ("arg=%08x sin=%08x cos=%08x\n", arg, sine, cosine);
}
if (sine == refsi) sin_correctly_rounded++;
if (cosine == refci) cos_correctly_rounded++;
errs = fabs (s8p23_to_double (sine) - refs);
errc = fabs (s8p23_to_double (cosine) - refc);
errsqs += errs * errs;
errsqc += errc * errc;
if (errs > maxerrs) maxerrs = errs;
if (errc > maxerrc) maxerrc = errc;
}
printf ("max err sine = %15.8e max err cosine = %15.8e\n",
maxerrs, maxerrc);
printf ("rms err sine = %15.8e rms err cosine = %15.8e\n",
sqrt (errsqs / arg), sqrt (errsqc / arg));
printf ("sin correctly rounded: %d (%.2f%%) cos correctly rounded: %d (%.2f%%)\n",
sin_correctly_rounded, 1.0 * sin_correctly_rounded / arg * 100,
cos_correctly_rounded, 1.0 * cos_correctly_rounded / arg * 100);
return EXIT_SUCCESS;
}
Die Ausgabe des beiliegenden Testframeworks sollte im Wesentlichen so aussehen:
S8.23 fixed-point sincos NJ variant
arg=00000000 sin=00000000 cos=00800000
arg=00000001 sin=00000001 cos=00800000
arg=00000002 sin=00000002 cos=00800000
arg=00000003 sin=00000003 cos=00800000
arg=00000004 sin=00000004 cos=00800000
arg=00c90fd6 sin=00800000 cos=00000005
arg=00c90fd7 sin=00800000 cos=00000004
arg=00c90fd8 sin=00800000 cos=00000003
arg=00c90fd9 sin=00800000 cos=00000002
arg=00c90fda sin=00800000 cos=00000001
arg=00c90fdb sin=00800000 cos=00000000
max err sine = 9.86237533e-008 max err cosine = 1.02729891e-007
rms err sine = 4.11141973e-008 rms err cosine = 4.11752018e-008
sin correctly rounded: 10961278 (83.19%) cos correctly rounded: 11070113 (84.01%)
fixed-point sincos OllyW variant
arg=00000000 sin=00000000 cos=00800000
arg=00000001 sin=00000001 cos=00800000
arg=00000002 sin=00000002 cos=00800000
arg=00000003 sin=00000003 cos=00800000
arg=00000004 sin=00000004 cos=00800000
arg=00c90fd6 sin=00800000 cos=00000005
arg=00c90fd7 sin=00800000 cos=00000004
arg=00c90fd8 sin=00800000 cos=00000003
arg=00c90fd9 sin=00800000 cos=00000002
arg=00c90fda sin=00800000 cos=00000001
arg=00c90fdb sin=00800000 cos=00000000
max err sine = 1.13173883e-007 max err cosine = 1.13158773e-007
rms err sine = 4.30955921e-008 rms err cosine = 4.31472191e-008
sin correctly rounded: 10844170 (82.30%) cos correctly rounded: 10855609 (82.38%)
Chris Stratton
Olin Lathrop
Bob
Olin Lathrop
Bob
Chris Stratton