Hallo,
ich benötige für mein Programm eine Funktion, die mir die Wurzel aus
einer unsigned long variablen zieht. Bei der Suche im Forum bin ich nur
auf Assembler-Programme gestoßen.
Hat jemand so eine Funktion schon mal in C umgesetzt? Wenn ja wäre es
super, wenn ihr sie mir einmal schicken könntet? Würdet mir damit echt
weiterhelfen!
Hier wird ihnen geholfen, falls du englisch beherrschst kannst du dabei
sogar noch ein paar Sachen lernen, ansonsten belasse es bei Copy&Paste
...
http://www.embedded.com/98/9802fe2.htm
Samarium Kobald schrieb:
> Hallo,>> ich benötige für mein Programm eine Funktion, die mir die Wurzel aus> einer unsigned long variablen zieht. Bei der Suche im Forum bin ich nur> auf Assembler-Programme gestoßen.>> Hat jemand so eine Funktion schon mal in C umgesetzt? Wenn ja wäre es> super, wenn ihr sie mir einmal schicken könntet? Würdet mir damit echt> weiterhelfen!
Ich hab das mal in C gemacht exemplarisch:
1
// Für das Ergebnis der Wurzelberechnung
2
typedefstruct
3
{
4
// Das Endergebnis ist
5
// .wurzel / (2 ** .shift)
6
// d.h. .wurzel muss um .shift Binärstellen nach rechts geschoben werden,
7
// um die Quadratwurzel zu bekommen.
8
// .shift ist immer >= 0
9
// Durch diese Darstellung geht keine der berechneten Stellen verloren.
10
// Soll das Ergebnis die Quadratwurzel im Q-Format
11
// mit 2N Nachkommastellen sein, dann ist das Ergebnis
12
// .wurzel * (2 ** (N - .shift))
13
unsignedlongwurzel;
14
signedintshift;
15
16
// Falls .exact == 1 ist, dann ging die Wurzel-Berechnung ohne Rest auf.
17
// Die Eingabe — als nicht-negative ganze Zahl betrachtet — war eine Quadratzahl.
18
// Für .exact == 0 gilt:
19
// .wurzel führt zu einem Wert kleiner als die exakte Quadratwurzel
20
// .wurzel+1 führt zu einem Wert größer als die exakte Quadratwurzel
21
intexact;
22
}wurzel_t;
23
24
// Berechnet die Quadratwurzel mit 32 Bit Genauigkeit in der
25
// eben beschriebenen Darstellung.
26
// Der Algorithmus ist eine direkte Übertragung des in
27
// http://www.diaware.de/html/wurzel.html
28
// erklärten schriftlichen Wurzelziehens vom Zehner- ins Zweiersystem.
29
// Sprach-Standard ist GNU-C (gcc -std=gnu99 ...)
30
wurzel_tsqrt_32(unsignedlongr)
31
{
32
// MAXWERT hat das zweit-oberste Bit gesetzt
33
// Werte >= MAXWERT mit 4 zu multiplizieren führt zu einem Überlauf,
34
// was für Werte < MAXWERT nicht der Fall ist.
35
constunsignedlongMAXWERT=1<<(8*sizeof(long)-2);
36
37
// Ergebnis mit 0 vorbesetzen
38
wurzel_ts={.wurzel=0,.shift=0,.exact=0};
39
40
// r wird rechts um 32 binäre Nachkommastellen erweitert.
41
// Durch die Union f kann einfach nach rechts/links geschoben werden,
42
// da Vor- und Nachkommastellen zu 64 Bits zusammengefasst sind.
43
union
44
{
45
struct{unsignedlongfrac,r;};
46
unsignedlonglongr_frac;
47
}
48
f={{.r=r,.frac=0}};
49
50
// Spezialfall: r == 0
51
{
52
s.exact=1;
53
returns;
54
}
55
56
// r solange um 2 Stellen nach rechts schieben, bis die
57
// Wurzel einstellig (also 1) ist
58
while(f.r>=4)
59
{
60
// Anzahl Rechtsshifts und hinausgeschobene
61
// Ziffern werden für später gemerkt
62
s.shift--;
63
f.r_frac>>=2;
64
}
65
66
// Die erste Ziffer der Wurzel ist 1
67
s.wurzel=1;
68
f.r--;
69
70
// Schleifen, bis die Wurzel aufgeht oder die maximale
71
// Genauigkeit erreicht wurde. Mit jedem Durchlauf gewinnt
72
// das Ergebnis eine Stelle hinzu.
73
while(1)
74
{
75
// Falls der Rest gleich 0 ist, dann geht die Wurzel auf
76
if(f.r_frac==0)
77
{
78
s.exact=1;
79
break;
80
}
81
82
if(s.wurzel>=MAXWERT)
83
// Mehr Stellen passen nicht ins Ergebnis:
84
// q2 unten würde überlaufen
85
break;
86
87
// Eine weitere Stelle zu .wurzel hinzufügen
88
s.wurzel<<=1;
89
s.shift++;
90
91
if(f.r>=MAXWERT)
92
{
93
// Der Rest würde unten überlaufen; wir sind fertig
94
// q2 passt aber noch 1x in r hinein
95
s.wurzel++;
96
break;
97
}
98
99
// 2*.wurzel + 1 ergibt sich aus der Binomischen Formel
100
unsignedlongq2=(s.wurzel<<1)+1;
101
102
// Zwei (Binär)-Ziffern aus f.frac zu r hinzuholen
103
f.r_frac<<=2;
104
105
// passt q2 in r hinein?
106
if(q2<=f.r)
107
{
108
// Ja: Diese Stelle der Wurzel ist 1 (ansonsten 0)
109
// q2 von r abziehen
110
s.wurzel++;
111
f.r-=q2;
112
}
113
}
114
115
// .shift nicht-negativ machen
116
if(s.shift<0)
117
{
118
s.wurzel<<=-s.shift;
119
s.shift=0;
120
}
121
122
// Ergebnis zurückliefern
123
returns;
124
}
125
126
// Wandelt die oben erklärte Darstellung nach double
127
doublewurzel_to_double(wurzel_tw)
128
{
129
// Die Wurzel so hinschieben, daß sie 32 binäre
130
// Nach- und Vorkommastellen hat
131
unsignedlonglongll=w.wurzel;
132
ll<<=32-w.shift;
133
134
// Die Multiplikation mit 2**32 rückgängig machen.
135
return(double)ll/(1ULL<<32);
136
}
Brauchbar ist auch
http://www.restena.lu/convict/Jeunes/Math/square_root_CORDIC.htm
wobei da multipliziert werden muss, und man muss auf Überläufe
aufpassen. Vielleicht lässt sich aber die Verwendung von 64-Bit
Arithmetik vermeiden, wenn man die Überläufe im Griff hat (hab mir noch
keine Gedanken drum gemacht).
Johann
Johann
>ich benötige für mein Programm eine Funktion, die mir die Wurzel aus>einer unsigned long variablen zieht. Bei der Suche im Forum bin ich nur>auf Assembler-Programme gestoßen.
Du kannst doch Assembler Programme in C aufrufen. Nimms du dir die
Assemblerfunktion und baust dir ein C Aufruf drumrum.
Gruss Helmi
Schade, daß man das nicht so einfach für Q-Format übernehemn kann :-(
Ich bräucht auch ne schnellere Wurzel für Fixpunkt-Zahlen mit 4 Vor- und
12 Nachkommastellen.
Johann
Ersetze testweise einfach Root <<= 8 durch Root <<= 6 für dein
Festkommaformat. Es ist schon etwas her, seit ich das gemacht habe, aber
mit etwas Glück tut es ...
Hier die Wurzelfunktion in Assembler. Stammt aus einer APP Note von
Atmel.
Ich verwende sie um einen Effektivwert auszurechnen (und das bis zu 5000
mal in der Sekunde). Dabei ist der AVR immer noch nicht voll ausgelastet
und kann noch ein paar Filter mitberechnen. Der AVR (ATMEGA88) ist mit
16MHz getaktet. Du brauchst dir nur noch fuer deinen C Compiler einen
Aufruf zu schreiben und das wars. Fuer eventuelle Fixpoint Formate muss
du nur am Ergebnis noch ein paar Bits schieben.
Gruss Helmi
;-----------------------------------------------------------------------
---
;
; R17:R16=SQRT(R23:R22:R21:R20) rounded to the nearest integer (0.5
rounds up)
; Destroys the argument in R23:R22:R21:R20 and R0, R1
;
;-----------------------------------------------------------------------
---
Sqrt32: ldi R19,0xc0
clr R18 ; rotation mask in R19:R18:R0
ldi R17,0x40
clr R16 ; developing sqrt in R17:R16:R1
mul R16,R16 ; R0=R1=0, C=0
_sq32_1: brcs _sq32_2 ; C --> Bit is always 1
cp R21,R1 ; Extra instruction
cpc R22,R16
cpc R23,R17 ; Does test value fit?
brcs _sq32_3 ; C --> nope, bit is 0
_sq32_2: sub R21,R1 ; Extra instruction
sbc R22,R16
sbc R23,R17 ; Adjust argument for next bit
or R1,R0 ; Extra instruction
or R16,R18
or R17,R19 ; Set bit to 1
_sq32_3: lsr R19
ror R18 ; Shift right mask
ror R0 ; Extra instruction
eor R1,R0 ; Extra instruction
eor R17,R19
eor R16,R18 ; Shift right only test bit in
result
lsl R20
rol R21
rol R22
rol R23 ; Shift left remaining argument (C
used at _sq32_1)
sbrs R0,5 ; Skip if 17 bits developed
rjmp _sq32_1 ; Develop 17 bits of the sqrt
lsl R1 ; C if round up
adc R16,R19
adc R17,R19 ; Round up if C (R19=0)
ret
Was bedeutet denn der Kommentar "Extra Instruction"?
Leider funzt das auch nicht für Q-Format. Wurzel aus 2 liefert nämlich
1. Und so simpel ist das auch nicht anzupassen.
Mein Algorithmus von oben liefert als Ergebnis hingegen 46340 und sagt,
daß das noch um 15 Stellen nach rechts zu schieben ist um Wurzel 2 zu
bekommen. Wie weit man wirklich nach rechts schiebt hängt dann davon ab,
wo man das Komma haben will. Für 12 Nachkommastellen würde man als um 3
nach rechts schieben und hätte als Ergebnis
>Was bedeutet denn der Kommentar "Extra Instruction"?
Wie gesagt stammt aus der APP Note von Atmel
>Leider funzt das auch nicht für Q-Format. Wurzel aus 2 liefert nämlich>1. Und so simpel ist das auch nicht anzupassen.
Das ist Integermaessig auch richtig.
Du must vorher dein Argument entsprechend vorher nach links shiften um
mehr stellen zu bekommen. Du bekommst die Wurzel nur Integermaessig
berechnet.
Beim Q Format nimmst du ja auch nicht die niederigsten Stellen sondern
die hoechsten und justierst (shiftest) das Ergebnis wieder nach deinem Q
Format.
Von daher ist deine Eingabe von 2 auch kein Q Format.
Johann L. schrieb:
> Was bedeutet denn der Kommentar "Extra Instruction"?
Ich hab eine Theorie.
Dieser Kommentar taucht nur dann auf, wenn eine
'Mehr-Register-Operation' sich über mehrere Opcodes hinzieht.
cp R21,R1 ; Extra instruction
cpc R22,R16
cpc R23,R17 ; Does test value fit?
In diesem Abschnitt wird also geprüft 'Does test value fit?' und die
Sequenz beginnt bei 'Extra instruction'.
Dieses 'Extra instruction' ist also wie ein 'Achtung: Wenn du das
verstehen willst, musst du auch noch die nächsten par Anweisungen mit
dazunehmen' zu sehen.
Ich werd das bei Gelegenheit mal testen, ob's schneller ist als meine
Impementierung (kosten beide in der Größenordnung von 2 16/16
Divisionen)
Derweil hab ich's schon mal so umgeschrieben daß es zum avr-gcc ABI
passt. Und hoffe ich hab keine Tippfehler drinne :-)
Johann
h
Ist das nur bei mir so, dass in math.h die sowieso auf 32 Bit getrimmt
sind?
Benutze GCC für ARM7. Funzt einwandfrei bis 32 Bit
math.h:
float sqrtf(float);
double sqrt(double);
Helmut Lenzen schrieb:
> Falls noch jemand eine schnelle 24x24Bit Signed Multiplikationsroutine> mit 48Bit Ergebnis braucht:
Liebt von Dir.
Eine Bitte hab ich (ist jetzt nicht nur auf Helmut bezogen)
Wer solche Routinen der Allgemeinheit zur Verfügung stellen will, wärt
ihr so lieb und fügt sie auch ins Wiki ein?
Es gibt da eine Seite (die irgendwann mal Bestandteil des Turials wird)
auf der solche Basisfunktionen gesammelt werden.
Ich such die auch immer. Der für mich einfachste Weg, die Seite zu
finden sieht so aus:
Im Assembler Tutorial gibt es eine Seite: Arithmetik - Einfache
Grundoperationen
Die ist zunächst nur für 8-Bit Operationen gedacht und dafür ein
gewisses Grundverständnis aufzubauen, wie eigentlich Multiplikation,
Division, 2-er Komplement etc. funktionieren.
Aber: Ganz unten auf dieser Seite gibt es einen kurzen Abschnitt und
einen Link auf eine Folgeseite die sich mit Mehrbyteoperationen
beschäftigt. Und ich wäre sehr dankbar, wenn sich solche Codestückchen
dort sammeln würden. Es braucht nichts dokumentiert werden, es reicht
völlig, wenn das Codestückchen inhaltlich an die richtige Stelle
eingefügt wird.
Danke
und auch Danke, dass wieder ein Codestückchen mehr dort steht.
Irgendwann bleibt es nicht aus, die Codestückchen dort etwas zu
vereinheitlichen (Registerbenutzung, Aufrufkonventionen, etc) und die
Wirkungsweise zu dokumentieren. Dann wird daraus ein neuer Abschnitt
fürs Tutorial. Aber bis dahin bin ich schon dankbar, wenn die Sammlung
wächst.
Ich hab die AVR-GCC-ABI-sqrt32_avr-Funktion gerade mal versucht zu
testen.
Hab die sqrt32.c und sqrt32.S von Johann L. einfach mal so übernommen.
AVR-Studio hat das Makefile von selbst entsprechend ergänzt.
Fehlermeldung:
/main.c:12: undefined reference to `sqrt32_avr(unsigned long)'
Code sieht etwa so aus:
Bei mir tut's. Allerdings mit selbst geklöppeltem Makefile und ohne
AStudio.
Sieht aus als wäre der Header nicht so wie er sein soll, da es eine
Compilerfehler ist und kein Linkerfehler.
Im Header fehlt natürlich noch ein #ifndef SQRT32_H etc, was bei dem
einfachen Beispiel hier aber nicht stört.
Was heißt eigentlich das "a" an zero-reg? Wundert mich, daß das
überhaupt im Object landet.
Nochwas: Wenn ich ne foo.c Datei hab, mach ich in's Projekt keine foo.S,
weil Windoofs Huddel macht wenn es eine foo.s gibt (zB gcc -save-temps).
Johann
Hab das Makefile mal in das Stammverzeichnis des Projektes geschoben,
nachdem ich mal das Makefile-Template von WinAVR genommen hab
(standardmäßig liegt's in default/). Geht jetzt. Danke trotzdem.
Das mit foo.c geht schon klar - Gibt ja keine sqrt32.c.
Und soweit ich das sehe war das schon ein Linkerfehler. "... has not
been declared" würde der Compiler wohl sagen.
Gruß
Dennis
Nochwas: Hat mal jemand ne Speed-Messung gemacht für den
Atmel-Algorithmus bzw. die gcc-Variante davon?
Im Wiki wär ne ungefähre Obergrenze für die Anzahl der Ticks interessant
zu lesen.
Beim Überfliegen komm ich grob auf 450 Ticks WCET überschlagsmässig.
Johann
Das hier sieht relativ ähnlich der Funktion von Atmel aus:
http://elm-chan.org/docs/avrlib/sqrt32.S
Ich frage mich nur, wieso diese Funktion größer ist als die Lösung von
Atmel, obwohl sie anscheinend den gleichen Algorithmus verwendet.
Rechnet der vielleicht exakter?
Benedikt K. schrieb:
> Rechnet der vielleicht exakter?
Keine Ahnung...
Für die Wurzel hat man prinzipiell zwei denkbare Ergebnisse: Abrunden
oder Rundung zum Nächsten (Atmel). Die beiden Möglichkeiten sollten sich
nicht allzusehr in der Implementierung unterscheiden, schon garnicht in
deren Komplexität. Beim Atmels Wurzel fallen wahrscheinlich nur die
letzten paar Instruktionen weg bei Abrunden .
Die Atmel-Lösung hab ich verwendet um Julia-Mengen in Echtzeit mit den
AVR auf ne Oszi-Röhre zu pinseln. Sieht soweit gut aus und auch der
Graph der komplexen Wurzel stellt der ATmega168 wie erwartet dar.
Ich würd eher vermuten daß die Atmel-Jungs was mehr Hirnschmalz
investiert haben.
Johann
Helmut Lenzen schrieb:
> ;----------------------------------------------------------------------- ---> ;> ; R17:R16=SQRT(R23:R22:R21:R20) rounded to the nearest integer (0.5> rounds up)> ; Destroys the argument in R23:R22:R21:R20 and R0, R1> ;> ;----------------------------------------------------------------------- ---
Hallo,
wo gibt's ne Beschreibung, wie man Assembler-Code in C-Files aufruft?
Muss man die Register(R23:R22:R21:R20) in C dann erst selbst befüllen
und dann den ASM-Codeteil aufrufen? Muss man sonst noch irgendwelche
Register sichern?
Danke & Gruß
Micha :)
Michael Förster schrieb:> Hallo,> wo gibt's ne Beschreibung, wie man Assembler-Code in C-Files aufruft?> Muss man die Register(R23:R22:R21:R20) in C dann erst selbst befüllen> und dann den ASM-Codeteil aufrufen? Muss man sonst noch irgendwelche> Register sichern?
Konkret für deinen Fall gibt's eine Implementierung in
http://www.mikrocontroller.net/articles/AVR_Arithmetik#avr-gcc_Implementierung_.2832_Bit.29
Natürlich muss eine solche asm-Funktion so geschrieben sein, daß sie das
gleiche Interface (ABI) verwendet, wie es avr-gcc erwartet. Den
Rückgabewert also in R25:R24 anstatt in R17:R16.
Die Funktion wird über einen C-Header deklariert und wie gewohnt von C
aus aufgerufen, ohne Inline-Asm-Zauberei.
Falls die Funktion nicht ABI-konform ist, müssen Register umkopiert bzw.
gesichert werden, und zwar so, daß der Code ABI-konform wird. Ansonsten
fliegt dir eher früher als später alles um die Ohren.
Die Anpassung ans ABI erfolgt in Assembler, wobei das in Inline-Asm
ziemliches Gefrickel wäre weil es keine Constraints für einzelne
Register gibt.
Teilweise ist das ABI beschrieben in
http://www.nongnu.org/avr-libc/user-manual/FAQ.html#faq_reg_usage