Forum: Mikrocontroller und Digitale Elektronik [AVR Assembler] Wie rechnet man x^Kommazahl?


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
von Hans (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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?

von Marian E. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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

von A. S. (achs)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Thorsten S. (thosch)


Bewertung
2 lesenswert
nicht lesenswert
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

von Hp M. (nachtmix)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Sven S. (boldie)


Bewertung
0 lesenswert
nicht lesenswert
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.

von A. B. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Johann L. (gjlayde) Benutzerseite


Bewertung
4 lesenswert
nicht lesenswert
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...

von Dr. Sommer (Gast)


Bewertung
-1 lesenswert
nicht lesenswert
Alternativ die pow() Funktion der C Library aufrufen. Das geht auch von 
Assembler aus...

von Johann L. (gjlayde) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
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.

von Jobst Q. (joquis)


Bewertung
1 lesenswert
nicht lesenswert
Zum Wurzelziehen mit Integer-Zahlen kann man auch die Sukzessive 
Approximation verwenden, wie sie bei Analog-Wandlern eingesetzt wird.

von bastler (Gast)


Bewertung
0 lesenswert
nicht lesenswert

von Waldo (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Exp(Ln(x)*a) = x^a

von Peter D. (peda)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Possetitjel (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Georg (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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

von Detlev T. (detlevt)


Bewertung
0 lesenswert
nicht lesenswert
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)

von Johann L. (gjlayde) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
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.

von Possetitjel (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Hans (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Hans (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Nachtrag: Ich habe im letzten Beitrag den Exponent 1.1 genommen, da 
dieser das leichteste Beispiel für die angesprochene Methode ist.

von Theor (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

von Yalu X. (yalu) (Moderator)


Bewertung
1 lesenswert
nicht lesenswert
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
von Hans (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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......

von Dr. Sommer (Gast)


Bewertung
1 lesenswert
nicht lesenswert
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)

von Hans (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Hmmm so langsam machts vielleicht doch klick. Rechne ich 1,1₁₀ um, indem 
ich 256/10 teile und das dann binär darstelle?

von Hans (Gast)


Bewertung
0 lesenswert
nicht lesenswert
OK Danke... habe mich nie mit Nachkommastellen im binären System 
beschäftigt :/

von Dr. Sommer (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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...

von W.S. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
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.

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.