Hallo, für ein Programm muss ich eine Formel implementieren, die x^1.2 rechnet. x ist ein 16 Bit Wert. Ich könnte eine LUT machen, aber da ich schon öfters vor diesem Problem (mit anderen reellen konstanten Exponenten) stand, wollte ich mal versuchen, die mathematische Funktion zu implementieren. Leider habe ich keine Ahnung, wie man das jetzt rechnet. Ich kann mir auch nichts darunter vorstellen. Hat jemand eine Idee?
Das Stichwort lautet Potenzfunktion mit gebrochen rationalem Exponent. y=x^n/m -> kann also als n-te Wurzel aus x^m ausgedrückt werden. Marian
Das ist x * 5teWurzel(x). (x*x^(1/5))=x^(1+1/5)=x^1.2 Am einfachsten (für alle Exponenten) gehts natürlich per Logarithmen (-; Die 5te Wurzel kannst Du ähnlich finden, wie die normale (sqrt). Nur das dafür vermutlich kaum geniale Optimierungen existiren. Ich würde es jetzt ganz einfach anfangen wie bein normalen Wurzeln: y=x/5. Wenn y*y*y*y*y > x, dann y/=5, sonst y*=1.2; Und dann per Excel oder C-Script für die 65000 Werte die Faktoren optimieren.
Marian E. schrieb: > y=x^n/m -> kann also als n-te Wurzel aus x^m ausgedrückt werden. doch wohl eher: m-te Wurzel aus x^n
Thorsten S. schrieb: > doch wohl eher: m-te Wurzel aus x^n Und eben dieses x^n führt bei einer 16-Bit Darstellung schnell zu Problemen durch Überlauf. Deshalb ist es besser die bereits erwähnten Logarithmen zu verwenden.
Wie genau muss die Lösung sein? Kann es auch eine Näherung sein? Ich hatte auch schon solche "Probleme" und da war ein Näherungswert mit z.B.<= 1% Fehler ausreichend und ich konnte eine Reihenentwicklung innerhalb eines Bereichs den ich benötigt nutzen. Das ist manchmal eine gute Lösung, wenn man den Wertebereich und den zu tolerierenden Fehler kennt.
Wie oben schon erwähnt ist x^1.2 = x * x^(1/5). Die 5-te Wurzel läßt sich einfach mittels ein paar Schritten des Newton-Verfahrens ermitteln: Ist x die 5-te Wurzel von a (>0), so ist x eine Nullstelle von f(x)=x^5-a. Darauf Newton anwenden, ein bisschen Bruchrechnung, und schon hat man die Iterationsformel. Ganz analog zur Quadratwurzel, dafür findet sich eine optimierte Assembler-Version z. B. im Linux-Kernel.
Mit 5. Wurzel bzw. Newton ist kein sonderlich praktikabler Ansatz aus 2 Gründen: 1) Es ist keine "einfache" Lösung, die man einfach so hinschreibt. Du wirst Fixed-Point verwenden wollen, und die benötigte Genauigkeit bekommst du nur nach Analyse der Konditionierung. Dito für die benötogte Skalierung / Wertebereicht. Es bringt nix, 16-Bit Fixed zu implementieren, nur um danach zu sehen, dass es interne Überläufe gibt die das Ergebnis schrotten. Zudem braucht's Dvisions- und Multiplikationsalgorithmen für entsprechend skalierte Ein- und Ausgaben. Ob je eine Implementierung ausreicht ist a priori nicht klar. Je größer der Wurzelexponent, desto übler die Situation und desto langsamer konvergiert Newton. 2) Es ist keine allgemeine Lösung, ist also nicht ohne weiteres auf andere Exponenten übertragbar. Zunächst sollte geprüft werden, ob Lookup-Tables verwendet werden können, und ob ein "normaler" Lookup ausreicht oder ob zusätzlich interpoliert werden soll. Ordnung 0: "Normaler" LUT. Man liest den Wert aus einer Tabelle — gegebenenfalls mit Skalierung von Ein- und Ausgabe — und verwendet den so erhaltenen Wert. Einfache Implementierung bei relativ großen Fehlern, d.h. die Dichte der LUT zu erhöhen führt nur zu einer entsprechenden (linearen) Verbesserung der Genauigkeit. Ordnung 1: Wie Ordnung 0 jedoch mit linearer Interpolation der Werte. Zunächst bestimmt man das Intervall, in dem die Eingabe x liegt, z.B. gelte x0 <= x <= x1. Dann ist x = (1-h)*x0 + h*x1 = x0 + h*(x1-x0). Ist y0 = f(x1) und y1 = f(x1) , dann berechnet sich y = f(x) := y0 + h*(y1-y0). Bessere Genauigkeit als Ordnung 0 bei etwas größerem Rechenaufwand. Ordnung 2: Interpolation mit kubischen Splines, entwickelt nach Bernstein-Polynomen und mit Kontrollpunkten bei 1/3 und 2/3 des jeweiligen Intervalls. Hört sich kompliziert an, ist aber recht einfach zu rechnen. Insbesondere sichert die Darstellung in Bernsteins, dass bei jedem Schritt mit Eingabe in [a,b] das Ergebnis wieder in [a,b] ist. Die Wahl der Kontrollpunkte ist einer möglichst einfachen Darstellung und Analyse geschuldet: Die Werte an den Kontrollpunkten ergeben sich, indem man die Tangente durch den nächsten Randpunkt nimmt und am Kontrollpunkt auswertet. Die Auswertung besteht dann aus 3 Schritten à la Ordnung 1. Hier ein Beispiel in Konkreter Implementierung für avr-gcc. Implementiert wurde es für arcsin nach Zerlegung von arcsin in Quadratwurzel und eine Funktion a(), die in [-1,3] analytisch ist. Entwickelt wird a:
1 | #include <stdint.h> |
2 | #include <stdfix.h> |
3 | |
4 | ...
|
5 | |
6 | extern unsigned _Fract __sqrt_ur_floor_mul (uint16_t); |
7 | extern unsigned _Fract __sqrt_ur_round_mul (uint16_t); |
8 | |
9 | typedef unsigned _Fract data_t; |
10 | |
11 | typedef struct |
12 | {
|
13 | // x Start, x End
|
14 | //data_t x_start, x_end;
|
15 | // # Intervals: N
|
16 | uint8_t n; |
17 | // # Intervals: N
|
18 | uint8_t log_n; |
19 | // Step: (end-start) / N
|
20 | //data_t step;
|
21 | // Data for function
|
22 | // # Entries in the table: 2 * (N+1)
|
23 | data_t data[10]; |
24 | } fun_data_t; |
25 | |
26 | // Define _() if you need some transfrom applied to the raw data
|
27 | // like, e.g. convert to some fixed point representation
|
28 | #ifndef _
|
29 | #define _(X) (data_t) X
|
30 | #endif // _
|
31 | |
32 | const __flash unsigned _Fract data_red_asin[] = |
33 | {
|
34 | // sqrt (32) * (a(x)-1)
|
35 | _(0.0000000000000000e+000), _(3.9283710065919304e-002) // _(0.0000000000000000e+000) ... |
36 | , _(1.2501973301494423e-001), _(4.4259597606718604e-002) // _(2.5000000000000000e-001) ... |
37 | , _(2.6698966805210750e-001), _(5.0677394156443352e-002) // _(5.0000000000000000e-001) ... |
38 | , _(4.3126310084303182e-001), _(5.9294454314661967e-002) // _(7.5000000000000000e-001) ... |
39 | , _(6.2633105768720609e-001), _(7.1533945534183907e-002) // _(1.0000000000000000e+000) |
40 | };
|
41 | |
42 | ...
|
43 | |
44 | data_t eval (const __flash uint8_t *data, data_t b, uint8_t interval) |
45 | {
|
46 | const __flash uint8_t *py = & data[interval]; |
47 | |
48 | // Get Control Points' y-Values
|
49 | |
50 | data_t y0_0, y0_1, y0_2, y0_3; |
51 | data_t dy0, dy3; |
52 | asm ("lpm %A0,Z+ $ lpm %B0,Z+" "\n\t" |
53 | "lpm %A1,Z+ $ lpm %B1,Z+" "\n\t" |
54 | "lpm %A2,Z+ $ lpm %B2,Z+" "\n\t" |
55 | "lpm %A3,Z+ $ lpm %B3,Z+"
|
56 | : "=r"(y0_0), "=r"(dy0), "=r"(y0_3), "=r"(dy3) , "+z" (py)); |
57 | |
58 | y0_2 = y0_3 - dy3; |
59 | y0_1 = y0_0 + dy0; |
60 | |
61 | // 3rd Order
|
62 | data_t y1_0 = y0_0 + b * dy0; |
63 | data_t y1_1 = y0_1 + b * (y0_2 - y0_1); |
64 | data_t y1_2 = y0_2 + b * dy3; |
65 | |
66 | data_t y2_1 = y1_1 + b * (y1_2 - y1_1); |
67 | data_t y2_0 = y1_0 + b * (y1_1 - y1_0); |
68 | |
69 | data_t y3 = y2_0 + b * (y2_1 - y2_0); |
70 | |
71 | return y3; |
72 | }
|
73 | |
74 | data_t eval_4 (const __flash data_t *fdata, data_t x) |
75 | {
|
76 | unsigned char interval = bitsur (x) >> 8; |
77 | interval >>= 4; |
78 | interval &= ~3; |
79 | x = urbits (4 * bitsur (x)); |
80 | return eval ((const __flash uint8_t*) fdata, x, interval); |
81 | }
|
82 | |
83 | |
84 | unsigned _Accum asin_ur (data_t z) |
85 | {
|
86 | if (z == 0ur) |
87 | return 0; |
88 | data_t x = -z; |
89 | |
90 | data_t a4sq2_x = eval_4 (data_red_asin, x); |
91 | |
92 | data_t sqrt_x = __sqrt_ur_round_mul (bitsur (x)); |
93 | |
94 | unsigned _Accum x32 = (unsigned _Accum) a4sq2_x + (unsigned _Accum) (4*M_SQRT2); |
95 | x32 = (unsigned _Accum) (2*M_PI) - x32 * sqrt_x; |
96 | return x32 >> 2; |
97 | }
|
98 | |
99 | static inline __attribute__((always_inline)) |
100 | _Accum neg_k (_Accum x) |
101 | {
|
102 | register _Accum rx asm ("22") = x; |
103 | extern uint32_t __negsi2 (uint32_t); |
104 | asm ("%~call %x1" : "+r" (rx) : "i" (__negsi2)); |
105 | return rx; |
106 | //return k_bits_ul (-ul_bits_k (x));
|
107 | }
|
108 | |
109 | _Accum asin_k (_Accum z) |
110 | {
|
111 | uint8_t neg = 0; |
112 | |
113 | if (z < 0k) |
114 | {
|
115 | z = neg_k (z); |
116 | neg = 1; |
117 | }
|
118 | |
119 | if (z > 1k) |
120 | return __ACCUM_MIN__; |
121 | |
122 | if (z == 1k) |
123 | {
|
124 | z = (_Accum) (M_PI/2); |
125 | }
|
126 | else
|
127 | {
|
128 | unsigned _Fract rz = (unsigned _Fract) z; |
129 | z = asin_ur (rz); |
130 | }
|
131 | |
132 | if (neg) |
133 | z = neg_k (z); |
134 | return z; |
135 | }
|
Die Tabelle ist so angelegt, dass die 2 MSBs das Intervall ergeben und die kleineren Bits die Position in Intervall. Die Werte an geraden Positionen sind die Funktionswerte. An ungeraden Stellen das Delta zum nächsten Kontrollpunkt (also dem bei 1/3 des Intervalls). Z.B. y0 = x0, Wert des 1. Kontrollpunkts y0_1 = y0 + dy0. Die Variablen y0_0, y0_1, y0_2, y0_3 bezeichnen ein Intervall: y0_0 und y0_3 Werte an den Endpunkten, y0_1 Wert am Kontrollpunkt bei 1/3 des Intervalls und y0_2 bei 2/3. Zu bemerken ist, das pro Intervall nur 2 Werte gespeichert werden müssen: y0_0 und y0_1 die auch als Tangente an f interpretiert werden können. Das rechte Intervallende gehört dann auch zum nächsten Intervall. Die Berechnung ist zwar in C, sollte aber selbst für einen Assembler-Progrmmieren offensichtlich genug sein :-) Der Vorteil ist, dass dies gut auf andere Funktionen verallgemeinert werden kann. Im Beispiel wird arcsin im Intervall (-1,1) auf die durch Embedded-C vorgeschriebene Genauigkeit berechnet, trotz der Singularität in -1 und 1! (Die Singularität stecht in der Wurzel, nicht in a()). Die Berechnung der Tabelle erfolgt natürlich auf dem Host. Ob eine Aufspaltung von f(x)=x^1.2 als f(x) = exp(1.2*ln(x)) = exp(1.2*ld(x)/ld(e)) besser ist, ist mitr nicht klar. Zwar lässt sich logarithmus dualis recht einfach binär entwickeln, siehe https://de.wikipedia.org/wiki/Logarithmus#Berechnung_einzelner_Bin.C3.A4rziffern Es verbleibt dann aber immer noch Skalierung und exp. Potenzreihenentwicklung ist wegen der auftretenden Divisionen, Auslöschung und möglichen internen Überläufe nur bedingt praktikabel bis unbrauchbar für Fixed-Point; man wird exp also ebenfalls z.B. per LUT darstellen wollen wie oben beschrieben. Vorteil von exp und ld ist deren Funktionalgleichung und daher einfache Skalierung, z.B. ld(2x) = 1+ld(x). Das verwendet man zur Eindämmung der auftretenden Werte(bereiche). ld ergibt im Exponenten übrigens eine andere Konstante, nämlich 1.2/ld(e) = 1.2*ln(2) = 0.8317...
Alternativ die pow() Funktion der C Library aufrufen. Das geht auch von Assembler aus...
Dr. Sommer schrieb: > Alternativ die pow() Funktion der C Library aufrufen. Das geht auch von > Assembler aus... Falls ein Link-fähiger Assembler mit gleichem Binärformat verwendet wird. Mit pow ist's auch nicht getan, es kommen Konvertierungen hinzu.
Zum Wurzelziehen mit Integer-Zahlen kann man auch die Sukzessive Approximation verwenden, wie sie bei Analog-Wandlern eingesetzt wird.
Hans schrieb: > Leider habe ich keine Ahnung, wie man das jetzt rechnet. Wer verbietet Dir denn, das einfach in C hinzuschreiben? Das Atmel Studio kannst Du doch für umme runterladen.
Hai! Sehr hübsch; Du hast das (viel besser) ausgeführt, worüber ich nur im stillen Kämmerlein theoretisiert habe. Johann L. schrieb: > Ordnung 2: Interpolation mit kubischen Splines, entwickelt > nach Bernstein-Polynomen und mit Kontrollpunkten bei 1/3 > und 2/3 des jeweiligen Intervalls. Hmm... das ist doch das, was man landläufig als Bezier-Kurve bezeichnen würde, oder liege ich falsch? > Hört sich kompliziert an, ist aber recht einfach zu rechnen. > Insbesondere sichert die Darstellung in Bernsteins, dass > bei jedem Schritt mit Eingabe in [a,b] das Ergebnis wieder > in [a,b] ist. Hmm. ??? Der atan(x) mit x aus [3;6] liegt aber nicht in [3;6]... ?!? Willst Du darauf hinaus, dass die Funktionswerte der Ersatz- funktion in den Stützstellen exakt sind? > Die Wahl der Kontrollpunkte ist einer möglichst einfachen > Darstellung und Analyse geschuldet: Die Werte an den > Kontrollpunkten ergeben sich, indem man die Tangente durch > den nächsten Randpunkt nimmt und am Kontrollpunkt auswertet. > Die Auswertung besteht dann aus 3 Schritten à la Ordnung 1. Hast Du ein für mich als Laien leicht einsehbares Argument, dass de Casteljau in der Auswertung günstiger ist als das Horner-Schema mit passend gewählten Koeffizienten? -- Mal abgesehen jetzt von der reinen Anzahl der zu speichernden Koeffizienten... # # # Ich hatte noch über zwei Verfeinerunge nachgegrübelt: Ordnung 3: Man muss die Ableitungen in den Stützstellen nicht gleich der Ableitung der Originalfunktion wählen, sondern könnte sehen, ob man das als Freiheitsgrad zur Erhöhung der Genauigkeit ausnutzen kann. Auf die Funktionsauswertung zur Laufzeit hat eine veränderte Lage der Steuerpunkte keinen Einfluss, aber es stellt sich die Frage, wie man in der Analysephase diese Parameter gewinnt. Ordnung 4: Man muss die Intervalle nicht gleich lang wählen. Mit lediglich 2 (bzw. 3) Vergleichen zur Laufzeit kann man den Definitionsbereich in 4 (bzw. 8) Teilintervalle frei wählbarer Länge zerlegen. Insbesondere für Funktionen mit horizontalen Asymptoten (wie dem arctan) ist das günstig, weil die Intervalllänge für große Argumente sehr stark zunehmen darf. Auch hier stellt sich die Frage, wie man die Stützstellen optimal wählt.
Hans schrieb: > wollte ich mal versuchen, die mathematische Funktion zu > implementieren Ansatz mit wenig Arbeit, dafür viel Programmspeicher: Die gewünschte Funktion in C hinschreiben und vom Assembler aus aufrufen. Der Linker wird die benötigten Teile der Gleitkomma-Bibliothek dazuladen. Ausserdem muss man sich mit den Aufruf-Konventionen beschäftigen und mit den nötigen Typumwandlungen, aber das ist auch einfach, wenn man die C-Funktion gleich mit Ein/Ausgabe von Integerwerten definiert, sind ja bloss 2 Parameter und 1 Rückgabewert. Wieviel Code entsteht und wie schnell die Berechnung ist: am besten ausprobieren. Ich habe ja keine Ahnung was dir zur Verfügung steht. Ich habe so (Mixed language programming) schon vor Jahrzehnten mit Z80 gearbeitet. Dass man das selbst besser hinkriegt als ein guter Compiler ist eher bei seltenen Spezialfällen zu erwarten. Georg
Für so etwas fällt mir das CORDIC-Verfahren ein. Schaue dir doch einmal AVR-CORDIC an. (Habe ich selbst aber noch nicht benutzt)
Possetitjel schrieb: > Johann L. schrieb: > >> Ordnung 2: Interpolation mit kubischen Splines, entwickelt >> nach Bernstein-Polynomen und mit Kontrollpunkten bei 1/3 >> und 2/3 des jeweiligen Intervalls. > > Hmm... das ist doch das, was man landläufig als Bezier-Kurve > bezeichnen würde, oder liege ich falsch? Bézier-Kurven sind allgemeiner. Hier geht es ja nicht darum, eine Kurve f(t)=(x(t),y(t)) zu zeichnen, sondern man will eine Abbildung von x auf y. Am einfachsten geht das, wenn man x(t) linear in t wählt. Das ergibt dann effektix eine Funktion y(x). Tatsächlich folgen Kontrollpunkte bei 1/3 bzw. 2/3 aus der Linearität von x(t) mit der Zusatzbedingung übereinstimmender Ableitung an den Randpunkten. >> Hört sich kompliziert an, ist aber recht einfach zu rechnen. >> Insbesondere sichert die Darstellung in Bernsteins, dass >> bei jedem Schritt mit Eingabe in [a,b] das Ergebnis wieder >> in [a,b] ist. Wenn zwei Kontrollpunkte y-Werte ya und yb haben, denn liegt der nächste Zwischenwert im Intervall [ya,yb] weil er sich errechnet zu (1-h)·ya+h·yb mit h in [0,1]. Gleiches gilt übrigens für die x-Koordinaten, aber weil x(t) linear gewählt wurde ist das trivial (hier bedeutet es, dass die x-Koordinaten der Kontrollpunkte für Intervall [xa,xb] in [xa,xb] liegen). > Willst Du darauf hinaus, dass die Funktionswerte der Ersatz- > funktion in den Stützstellen exakt sind? Das ist der einfachste Weg, Kontrolle über die Ordnung des Fehlers zu haben. Bei einer normalen LUT muss man den Abstand zwischen den Stützstellen halbieren bzw. die Anzahl der Stützstellen verdoppeln, um den Fehler zu halbieren. Proportionalitätsfaktor zwischen Stützabstand und Fehler ist i.w. die 1. Ableitung der Funktion. Bei linearer Interpolation wird der Fehler quadratisch kleiner mit dem Abstand der Stützstellen. Zu beachten ist, dass dies der Verfahrensfehler ist und nicht den Fehler durch Rundung beinhaltet. Als Beispiel für eine Analyse für lineare Interpolation siehe z.B. https://www.mikrocontroller.net/articles/AVR_Arithmetik/Sinus_und_Cosinus_(Lineare_Interpolation)#Aufbau_der_Tabelle >> Die Wahl der Kontrollpunkte ist einer möglichst einfachen >> Darstellung und Analyse geschuldet: Die Werte an den >> Kontrollpunkten ergeben sich, indem man die Tangente durch >> den nächsten Randpunkt nimmt und am Kontrollpunkt auswertet. >> Die Auswertung besteht dann aus 3 Schritten à la Ordnung 1. > > Hast Du ein für mich als Laien leicht einsehbares Argument, > dass de Casteljau in der Auswertung günstiger ist als das > Horner-Schema mit passend gewählten Koeffizienten? Ja. Die o.g. Eigenschaft, (1-h)·ya+h·yb im Intervall [ya,yb] liegt. Wichtig bei Fixed-Point. Bei meiner Implementierung ist allerdings darauf zu achten, dass dort nicht die y-Koordinaten der Kontrollpunkte gespeichert sind, sondern die Differenz zum linken y-Wert (Wegen Symmetrie ist nur der 1. Kontrollpunkt bei 1/3 zu speichern). Die Differenzen zu den Kontrollpunkten zu speichern ist also weniger allgemein, war aber angenehmer. > Ich hatte noch über zwei Verfeinerunge nachgegrübelt: > > Ordnung 3: Man muss die Ableitungen in den Stützstellen nicht > gleich der Ableitung der Originalfunktion wählen, sondern > könnte sehen, ob man das als Freiheitsgrad zur Erhöhung > der Genauigkeit ausnutzen kann. Auf die Funktionsauswertung > zur Laufzeit hat eine veränderte Lage der Steuerpunkte keinen > Einfluss, aber es stellt sich die Frage, wie man in der > Analysephase diese Parameter gewinnt. Zum einen das. Zum anderen musst du sicherstellen, dass x(t) linear ist — oder aufwändig x invertieren um von x auf t zu kommen (y ist eine Funktion von t! Dank linearem x wurde nur alles so vereinfacht, dass kein t mehr auftaucht). > Ordnung 4: Man muss die Intervalle nicht gleich lang wählen. Das ändert nicht die Dynamik des Fehlers in Abhängigkeit der Intervallänge. Ohne Kenntnis der Dynamik kannste nur rumprobieren. > Mit lediglich 2 (bzw. 3) Vergleichen zur Laufzeit kann man > den Definitionsbereich in 4 (bzw. 8) Teilintervalle frei > wählbarer Länge zerlegen. Binäre Suche also. Außerdem ist die Darstellung nicht mehr so kompakt, weil Symmetrie verloren geht, d.h. für jedes Intervall muss man beide Stützpunkte speichern. Das macht 50% mehr Daten und Zugriffe. Das Datenvolumen ist vermutlich nicht das Problem, weil aufgrund der Fehlerordnung wenige (im Vergleich zur vanilla LUT) Intervalle gebraucht werden. > Insbesondere für Funktionen mit horizontalen Asymptoten (wie dem > arctan) ist das günstig, Bei atan würd ich stark dazu tendieren, Funktionalgleichung wie atan(x) + atan(1/x) = pi/2 anzuwenden um das zu interpolierende Intervall zu beschränken — trotz blöder Division.
Johann L. schrieb: > Bézier-Kurven sind allgemeiner. Hier geht es ja nicht > darum, eine Kurve f(t)=(x(t),y(t)) zu zeichnen, sondern > man will eine Abbildung von x auf y. Am einfachsten > geht das, wenn man x(t) linear in t wählt. Das ergibt > dann effektix eine Funktion y(x). Mich hatte nur die gemeinsame Erwähnung von "Spline" und "Bernstein-Polynom" aus dem Konzept gebracht, weil ich letztere nur aus dem Kontext von Bezier-Kurven kenne. > Tatsächlich folgen Kontrollpunkte bei 1/3 bzw. 2/3 aus > der Linearität von x(t) mit der Zusatzbedingung > übereinstimmender Ableitung an den Randpunkten. Okay... mal sehen, ob ich Dich verstanden habe: Du gehst von dem Formalismus der Bezier-Kurven (Bernstein-Polynome) aus, schränkst aber die Lage der Steuerpunkte so ein, dass nur Splines mit stetiger 1. Ableitung herauskommen können. >>> Hört sich kompliziert an, ist aber recht einfach zu rechnen. >>> Insbesondere sichert die Darstellung in Bernsteins, dass >>> bei jedem Schritt mit Eingabe in [a,b] das Ergebnis wieder >>> in [a,b] ist. > > Wenn zwei Kontrollpunkte y-Werte ya und yb haben, denn liegt > der nächste Zwischenwert im Intervall [ya,yb] weil er sich > errechnet zu (1-h)·ya+h·yb mit h in [0,1]. Ach so. Ich hatte Dein "in jedem Schritt" nicht verstanden. De Casteljau bestimmt ja, wenn ich mir das richtig gemerkt habe, immer nur gewichtete Mittelwerte; insofern kann der neu erhaltene Wert nie außerhalb liegen. >> Willst Du darauf hinaus, dass die Funktionswerte der Ersatz- >> funktion in den Stützstellen exakt sind? > > Das ist der einfachste Weg, Kontrolle über die Ordnung des > Fehlers zu haben. Bei einer normalen LUT muss man den > Abstand zwischen den Stützstellen halbieren bzw. die Anzahl > der Stützstellen verdoppeln, um den Fehler zu halbieren. > Proportionalitätsfaktor zwischen Stützabstand und Fehler > ist i.w. die 1. Ableitung der Funktion. Ja, soweit einleuchtend. Der Punkt, der meinen Kopf zum Rauchen bringt, ist folgender: Oben wurde ja in den Stützstellen nur Stetigkeit der 1. Ableitung gefordert, der Wert dieser 1. Ableitung aber noch nicht festgelegt. Da ist also noch ein Freiheitsgrad vorhanden. Die Folge davon müsste nach meiner Vorstellung sein, dass es (mindestens) einen weiteren Punkt in jedem Teilintervall gibt, in dem der vorgegebene Funktionswert exakt angenommen wird. (Stimmt diese Vorstellung überhaupt?) Wie muss man die 1. Ableitungen (in den Stützstellen) wählen, damit der Fehler möglichst gleichmäßig verteilt ist, d.h. damit die Maxima der Fehlerbeträge in allen Teilintervallen möglichst gleich und kleiner als ein vorgegebener Grenzwert sind? Das Problem ist ja, dass ich diese Ableitung zwar wählen kann, mit dieser Wahl aber immer zwei Intervalle beeinflusse. Variation dieser Fragestellung: Genügt es, je eine weitere Stützstelle in jedem Intervall vorzuschreiben? Wo müsste diese liegen? In der Mitte jedes Inveralls? >> Hast Du ein für mich als Laien leicht einsehbares Argument, >> dass de Casteljau in der Auswertung günstiger ist als das >> Horner-Schema mit passend gewählten Koeffizienten? > > Ja. Die o.g. Eigenschaft, (1-h)·ya+h·yb im Intervall [ya,yb] > liegt. Wichtig bei Fixed-Point. Okay, das klingt einleuchtend. Muss ich mal in Ruhe überdenken. >> Ordnung 3: Man muss die Ableitungen in den Stützstellen >> nicht gleich der Ableitung der Originalfunktion wählen, >> [...] > Zum anderen musst du sicherstellen, dass x(t) linear ist Hier versagt meine Anschauung. Ich bin ganz naiv von folgender Überlegung ausgegangen: Bei einem kubischen Hermite-Spline ist, wenn ich mir das richtig gemerkt habe, die Kurve durch die Stützstellen und die Forderung nach Stetigkeit der 1. und 2. Ableitung festgelegt. Wenn man die 2. Ableitung ignoriert, muss also noch ein Freiheitsgrad vorhanden sein, den man zur Steigerung der Genauigkeit einsetzen könnte. Natürlich kann ich einfach den Wert der 1. Ableitung der Originalfunktion in den Stützstellen nehmen, aber ist diese Wahl optimal? Wenn nein: Wie geht es besser? >> Ordnung 4: Man muss die Intervalle nicht gleich lang wählen. > > Das ändert nicht die Dynamik des Fehlers in Abhängigkeit der > Intervallänge. Natürlich nicht, nein. Nur sagt diese Dynamik noch nicht unbedingt etwas über den Absolutwert des Fehlers. > Ohne Kenntnis der Dynamik kannste nur rumprobieren. Das ist ja das Problem - mir fehlt der geistige Aufhänger. Wahrscheinlich übersehe ich etwas ganz Elementares. Nützlich wäre ja schon eine Antwort auf folgende Frage: Ich habe im Intervall [a;b] eine Approximation mit Fehler delta. Um wieviel wächst schätzungsweise der Fehler, wenn ich das Intervall auf [a;r*b] mit r>1 vergrößere und die Polynomkoeffizienten anpasse? Funktionswert und 1. Ableitung des approximierenden Polynoms an der Stelle a sollen erhalten bleiben. >> Mit lediglich 2 (bzw. 3) Vergleichen zur Laufzeit kann man >> den Definitionsbereich in 4 (bzw. 8) Teilintervalle frei >> wählbarer Länge zerlegen. > > Binäre Suche also. Ja. > Außerdem ist die Darstellung nicht mehr so kompakt, weil > Symmetrie verloren geht, d.h. für jedes Intervall muss man > beide Stützpunkte speichern. Das macht 50% mehr Daten und > Zugriffe. Das Datenvolumen ist vermutlich nicht das Problem, > weil aufgrund der Fehlerordnung wenige (im Vergleich zur > vanilla LUT) Intervalle gebraucht werden. Ja, das ist für mich der Reiz der Idee: Drei Integer-Vergleiche ist im Verhältnis wenig, und die Tabelle bleibt, verglichen mit linearer oder gar keiner Interpolation winzig. Acht Intervalle mal vier Koeffizienten mal zwei Byte sind 64 Byte für die gesamte Tabelle. Trotzdem bleibt die Methode ziemlich universell. > Bei atan würd ich stark dazu tendieren, Funktionalgleichung > wie atan(x) + atan(1/x) = pi/2 anzuwenden um das zu > interpolierende Intervall zu beschränken — trotz blöder > Division. Ja... den Gedanken, atan(1/t) zu verwenden, hatte ich auch schon, aber so richtig schmeckt mir das nicht. Erstens ist das Verfahren dann nicht mehr universell, und zweitens expodiert der Laufzeitbedarf für kleine Controller.
Leute, ich versteh kein Wort! Ich bin kein Mathematiker. Vielleicht sollte ich jemanden die Funktion schreiben lassen ........ In einem anderen Thema steht folgendes: Josef G. schrieb: > Berechnung mittels Dualdarstellung des Exponenten: > > x^(b3b2b1.a1a2a3) = B3*B2*B1*A1*A2*A3 > > falls b3=0: B3=1 falls b3=1: B3=x^4 > falls b2=0: B2=1 falls b2=1: B2=x^2 > falls b1=0: B1=1 falls b1=1: B1=x^1 > falls a1=0: A1=1 falls a1=1: A1=sqrt(x) > falls a2=0: A2=1 falls a2=1: A2=sqrt(sqrt(x)) > falls a3=0: A3=1 falls a3=1: A3=sqrt(sqrt(sqrt(x))) Das verstehe* ich zwar, aber scheint falsch zu sein. Laut Taschenrechner ist 2^1.1=2.1435469. Nun ist es aber egal, wie ich es rechne, ich komme nicht auf den richtigen Wert. *Ich weiss, wie ich es im Taschenrechner eingeben muss.
Nachtrag: Ich habe im letzten Beitrag den Exponent 1.1 genommen, da dieser das leichteste Beispiel für die angesprochene Methode ist.
Hans schrieb: > Nachtrag: Ich habe im letzten Beitrag den Exponent 1.1 genommen, da > dieser das leichteste Beispiel für die angesprochene Methode ist. Das ist meiner Meinung nach, nicht leicht. Es gibt nämlich keine abbrechende Darstellung von 0.1. Der Nachkommateil sollte sich als Summe aus den Teilen von 1 darstellen lassen, die sich aus der Teilung durch 2-er Potenzen ergibt, damit das Beispiel einfach zu rechnen ist. Also etwa: 0.5, 0.25, 0,125 etc. Das ist bei 0.1 nicht der Fall.
Hans schrieb: > Das verstehe* ich zwar, aber scheint falsch zu sein. Doch, das stimmt schon. > Nun ist es aber egal, wie ich es rechne, ich komme nicht auf den > richtigen Wert. Hast du das berücksichtigt: Josef G. schrieb: > Berechnung mittels Dualdarstellung des Exponenten: ^^^^ Hans schrieb: > Nachtrag: Ich habe im letzten Beitrag den Exponent 1.1 genommen, da > dieser das leichteste Beispiel für die angesprochene Methode ist. Das ist nicht das leichteste Beispiel, da 1,1₁₀ = 1,00011001100110011…₂, also eine eine Zahl mit unendlich vielen Nachkommastellen ist. Probier besser mal 1,25₁₀ = 1,01₂ Edit: Zu lange mit dem Abschicken gewartet :)
:
Bearbeitet durch Moderator
Anscheinend habe ich schon nicht verstanden, was Dualdarstellung ist.
Ich dachte, es sei eine Binärdarstellung gemeint, so dass sich z.B. der
Exponent 4,3 zu 100,11 konvertiert. Das ist aber anscheinend ein grobes
Missverständnis meinerseits.
>1,25₁₀ = 1,01₂
Bahnhof......
Hans schrieb: > Das ist aber anscheinend ein grobes > Missverständnis meinerseits. In der tat... 100,11 ist: 1 Vierer 0 Zweier 0 Einer 1 Halbes 1 Viertel Also 100,11 (binär) = 4.75 (dezimal) Die Nachkommasstellen sind halt nicht zehntel/hundertstel usw, sondern halbe/viertel/achtel ... Und binär lässt sich ein Zehntel, also 0,1 (dez) nicht darstellen, genauso wie sich ein Drittel dezimal auch nicht darstellen lässt! (Mit endlich vielen Ziffern). Hans schrieb: >>1,25₁₀ = 1,01₂ > > Bahnhof...... Na 1,25 (dez) ist binär: 1 Einer 0 Halbe 1 Viertel also 1,01 (bin)
Hmmm so langsam machts vielleicht doch klick. Rechne ich 1,1₁₀ um, indem ich 256/10 teile und das dann binär darstelle?
OK Danke... habe mich nie mit Nachkommastellen im binären System beschäftigt :/
Hans schrieb: > OK Danke... habe mich nie mit Nachkommastellen im binären System > beschäftigt :/ Funktioniert genau wie im Dezimalsystem, aber eben Binär ;-) Eine Zahl der Gestalt
mit den Ziffern z_i (wobei z_0 die Ziffer direkt vor dem Komma ist) im Zahlensystem B (zB 2 oder 10) mit
hat den mathematischen Wert
Das B^k wird hier im Zehnersystem zu Zehnten,Hundertsteln,... im Binärsystem zu Halben, Vierteln...
Hans schrieb: > OK Danke... habe mich nie mit Nachkommastellen im binären System > beschäftigt Das ist auch nicht die Lösung für dein Eingangs-Problem. Du willst ja den Wert von a^x ausrechnen - und das wird ein bissel eklig (jedenfalls, wenn man's allgemein lösen will). 1. Problem: von a nach e, als a^x = e^(x*ln(a)) (wenn ich mich recht erinnere). Das ergibt zwei Berechnungen: ln(a) und e^y (mit y=x*ln(a)) 2. Problem: e^y berechnen. Das geht "relativ" einfach per pseudodivision und Pseudomultiplikation. Bedenke mal: e^(y1+y2+y3) = e^y1 * e^y2 * e^y3 Das sieht doch schon mal gut aus, gelle? Man kann also den Exponenten in additive Stücke zerlegen. Wenn man jetzt den Exponenten in solche Stücke zerlegt, daß e^Stück = 1 + ein Bruchteil wird, dann reduziert sich die Multiplikation in obiger Formel auf ne Verschiebung und Addition, gelle? Also, ein Beispiel (hier im Dezimalen gezeigt): y1 sei so groß, daß e^y1 = 1.1 ist y2 sei so groß, daß e^y2 = 1.01 ist y3 sei so groß, daß e^y3 = 1.001 ist y4 sei so groß, daß e^y4 = 1.0001 ist na und so weiter. Jetzt zerlegst du dein y in Stücke y1, y2, ... yn, indem du vom höchstwertigen y1 beginnend, versuchsweise alle yn abziehst und dir das Ergebnis in einer Zwischenvariablen merkst. Das ist die Pseudo-Division. Die Stücke y1..yn merkst du dir in einer Tabelle im ROM. Anschließend setzt du dein Ergebnis "E" erstmal auf 1.0 (weil ja e^0 = 1.0 ist) und machst dann folgendes: vom höchstwertigen Bit (dem von y1) in deinem Zwischenergebnis beginnend, addierst du (E>>Bitstelle), wenn das Bit 1 ist. Das ist die Pseudo-Multiplikation. Damit hättest du die Funktion e^y ausgerechnet und dabei nur Subtraktionen, Verschiebungen und Additionen gebraucht. Tja, das wäre die Hälfte der Übung, denn das Ausrechnen von ln(a) müßtest du in ähnlicher Weise zuvor erledigen, um von x auf x*ln(a) zu kommen. W.S.
Thread beobachten
Seitenaufteilung abschaltenBitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.