Wurzefunktion für uint32_t in C

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

typedef struct

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

    unsigned long wurzel;

14

    signed int shift;

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

    int exact;

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_t sqrt_32 (unsigned long r)

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

    const unsigned long MAXWERT = 1 << (8*sizeof (long) -2);

36

37

    // Ergebnis mit 0 vorbesetzen

38

    wurzel_t s = { .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 { unsigned long frac, r; };

46

        unsigned long long r_frac;

47

    }

48

    f = {{.r = r, .frac = 0}};

49

50

    // Spezialfall: r == 0

51

    {

52

        s.exact = 1;

53

        return s;

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

        unsigned long q2 = (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

    return s;

124

}

125

126

// Wandelt die oben erklärte Darstellung nach double

127

double wurzel_to_double (wurzel_t w)

128

{

129

    // Die Wurzel so hinschieben, daß sie 32 binäre 

130

    // Nach- und Vorkommastellen hat

131

    unsigned long long ll = 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

vielen dank für die schnelle antwort.

hab den code von alex's seite mal abgetippt (variante 4)

funktioiert einwandfrei!

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

Arne schrieb:
> siehe Anhang

Schnell ist das auch nicht, das ist einfach nur

Und Divisionen sind teuer (zumindest auf AVR).

Johann

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

Quadriert ergibt das rund 1.999

Johann

>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

1

#include <stdint.h>            

2

extern uint16_t sqrt32_avr (uint32_t);

S

1

__zero_reg__ = 1                     ; avr-gcc ABI

2

;; #include <stdint.h>            

3

;; extern uint16_t sqrt32_avr (uint32_t);

4

5

.global sqrt32_avr

6

.text

7

;--------------------------------------------------------------------

8

;  R25:R24 = SQRT (R25:R24:R23:R22) rounded to the nearest 

9

;  integer (0.5 rounds up)

10

;  Destroys R0, R18-R23, Clears R1

11

;  This is avr-gcc ABI conform

12

;--------------------------------------------------------------------

13

sqrt32_avr: ldi   R21,0xc0

14

            clr   R20               ; rotation mask in R19:R18:R0

15

            ldi   R19,0x40

16

            clr   R18               ; developing sqrt in R17:R16:R1

17

            mul   R18,R18           ; R0=R1=0, C=0

18

1:          brcs  2f                ; C --> Bit is always 1

19

            cp    R23,R1            ; Extra instruction

20

            cpc   R24,R18

21

            cpc   R25,R19           ; Does test value fit?

22

            brcs  3f                ; C --> nope, bit is 0

23

2:          sub   R23,R1            ; Extra instruction

24

            sbc   R24,R18

25

            sbc   R25,R19           ; Adjust argument for next bit

26

            or    R1,R0             ; Extra instruction

27

            or    R18,R20

28

            or    R19,R21           ; Set bit to 1

29

3:          lsr   R21

30

            ror   R20               ; Shift right mask

31

            ror   R0                ; Extra instruction

32

            eor   R1,R0             ; Extra instruction

33

            eor   R19,R21

34

            eor   R18,R20           ; Shift right only test bit in result

35

            lsl   R22

36

            rol   R23

37

            rol   R24

38

            rol   R25               ; Shift left remaining argument (C used at _sq32_1)

39

            sbrs  R0,5              ; Skip if 17 bits developed

40

            rjmp  1b                ; Develop 17 bits of the sqrt

41

            lsl   R1                ; C if round up

42

            clr   __zero_reg__      ; avr-gcc ABI

43

            adc   R18,__zero_reg__  ; Round up if C

44

            adc   R19,__zero_reg__

45

            movw  R24, R18          ; avr-gcc ABI

46

            ret

Falls noch jemand eine schnelle 24x24Bit Signed Multiplikationsroutine 
mit 48Bit Ergebnis braucht:

;*********************************************************************** 
***********************
;
;  muls24x24_48:
;  Signed Multiply von 2 24Bit breiten Zahlen mit 48Bit ergebnis
;
;  x                       = a           * b
;  R21:R20:R19:R18:R17:R16 = R27:R26:R25 * R24:R23:R22
;  hi                  lo    hi      lo    hi      lo
;
;*********************************************************************** 
***********************

muls24x24_48:
    clr    r2          ;Zero Register
    muls  r27,r24        ; (1) signed Multiply a(MSB) * b(MSB)
    mov    r21,r1
    mov    r20,r0

    mul    r26,r23        ; (2) unsigned
    mov    r18,r0
    mov    r19,r1

    mul    r25,r22        ; (3) unsigned multiply a(LSB) * b(LSB)
    mov    r17,r1
    mov    r16,r0

    mul    r26,r22        ;(4) unsigned
    add    r17,r0
    adc    r18,r1
    adc    r19,r2
    adc    r20,r2
    adc    r21,r2

    mul    r25,r23        ;(5) unsigned
    add    r17,r0
    adc    r18,r1
    adc    r19,r2
    adc    r20,r2
    adc    r21,r2

    push  r16
    push  r17

    mov    r16,r27
    mulsu  r16,r22        ;(6) unsigned * signed
    sbc    r20,r2
    sbc    r21,r2
    add    r18,r0
    adc    r19,r1
    adc    r20,r2
    adc    r21,r2

    mulsu  r16,r23        ;(7) unsigned * signed
    sbc    r21,r2
    add    r19,r0
    adc    r20,r1
    adc    r21,r2

    mov    r16,r24
    mov    r17,r25
    mulsu  r16,r17        ;(8) unsigned * signed
    sbc    r20,r2
    sbc    r21,r2
    add    r18,r0
    adc    r19,r1
    adc    r20,r2
    adc    r21,r2

    mov    r17,r26
    mulsu  r16,r17        ;(9) unsigned * signed
    sbc    r21,r2
    add    r19,r0
    adc    r20,r1
    adc    r21,r2

    pop    r17
    pop    r16
    ret


Gruss Helmi

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.

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

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 :)

dafür gibt's nen Artikel :)

http://www.mikrocontroller.net/articles/AVR-GCC-Tutorial#Assembler_und_Inline-Assembler

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