Forum: FPGA, VHDL & Co. Kubikwurzel bilden


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 BiBi (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Für die Quadratwurzel gibt es Cores und das Heronverfahren (siehe auch 
Quake-code), das im allgemeinen Fall ein Newtonverfahren ist und im 
Prinzip für jede Wurzel taugt.

Leider muss man je nach Güte des Schätzwertes mehrere mal iterieren, bis 
man zu einer guten Wurzel kommt. Da bei jeder Iteration eine Division 
auftaucht, taugt das nicht un bedingt für eine Implementierung in VHDL.

Wie berechnet man am schnellsten eine Kubikwurzel?

von Matthias N. (nippey)


Bewertung
0 lesenswert
nicht lesenswert
Da man das alles als Potenz sehen kann..
(Die 2. Wurzel aus x ist x hoch 1/2)
.. könntest du das alles mit Reihenentwicklungen machen.
Die Methode funktioniert in C gut, aber ich vermute, dass sie nicht 
wirklich für HDL geeignet ist..

[Ge-Guttenbergt von Wikipedia: 
http://de.wikipedia.org/wiki/Potenz_(Mathematik); 
http://de.wikipedia.org/wiki/Logarithmus]

Damit hättest du NUR noch Divisionen durch gerade Zahlen ;)
Viel Erfolg noch!

von BiBi (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Das mit der Potenz und Logarithmus hatte ich auch schon angedacht, 
allerdings hatte ich nicht an eine Reihenentwicklung gedacht. Das geht 
in VHDL durchaus einfach, aber die Nenner steigen nicht sonderlich 
schnell, weswegen das nicht gut konvergiert.

von Lothar M. (lkmiller) (Moderator) Benutzerseite


Bewertung
2 lesenswert
nicht lesenswert
BiBi schrieb:
> Wie berechnet man am schnellsten eine Kubikwurzel?
Weil du nur nach der Geschwindigkeit und nicht nach dem 
Ressourcenverbrauch gefragt hast: über eine Tabelle.

von Lustiger Denkermolch (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Mach doch einfach eine binäre Suche, oder auch eine "intelligente" 
interpolierende Suche. Gerade bei Ganzzahlen geht das schön flott; je 
höher die Ordnung der Wurzel, desto flotter...

von BiBi (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Ich würde gerne Kubikwurzeln aus Zahl ziehen, die wenigsten 32 Bit 
Breite haben - eher wahrscheinlich 48. Da wird es mit den Tabellen etwas 
schwer.

von D. I. (Gast)


Bewertung
0 lesenswert
nicht lesenswert
und wie genau brauchst du die wurzel? Also nur Ganzzahlanteil oder 
wieviele Nachkommastellen?

von Jürgen S. (engineer) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
In VHDL als tempoorientierte Lösung geht das am Besten über exp/log im 
Binärformat, also ld = logarithmus dualis, statt mit dem natürlichen ln.

Wenn es per SW sein soll, dann würde ich lieber gleich eine 
Potenzreihenentwicklung der dritten Wurzel nehmen.

von Dr. Mathematicus (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Hast Du vlt. auch Mathematica oder Maple? Dann könntest du über ein 
Funtionswindow auch Berechnungen dieser Programme durch führen lassen, 
da gibt es auch für C, Pascal, VDHL, Cobol, Fortran... fertige Patches. 
Was für ein Basisbetriebssystem hast Du?

von BaBa (Gast)


Bewertung
-2 lesenswert
nicht lesenswert
Ich programmiere mit CAN-Bus, geht dass damit auch?

von Dr. Mathematicus (Gast)


Bewertung
0 lesenswert
nicht lesenswert
BaBa schrieb:
> Ich programmiere mit CAN-Bus, geht dass damit auch?

Ja klar, darum heisst es ja cann-bus, weil man damit alles machen "can".

Damit kannst Du sogar Lichterketten mit ansteuern und aus der Anzahl der 
verwendeten Leuchtmittel die Kubikwurzel ziehen, dass ist doch toll, 
oder?

von ... (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Das Hauptproblem bei der Reihenentwicklung ist die Geschwidnigkeit, mit 
der der Wert gegen das Ergebnis konvergiert. Aus dem Grund normiert man 
den Wert und berechnet davon dann die eigentliche Funktion. Das ist 
alles sehr gut hier beschrieben:

yacas.sourceforge.net/Algo.book.pdf

Meistens funktioniert auch die Newton-Iteration sehr gut. Also einfach 
mit Newton-Iteration die Nullstelle von x^(1/3)-y berechnen. Einen guten 
Startwert kann man mit einem Polynom berechnen.

von Dr. Mathematicus (Gast)


Bewertung
0 lesenswert
nicht lesenswert
... schrieb:
> Meistens funktioniert auch die Newton-Iteration sehr gut. Also einfach
>
> mit Newton-Iteration die Nullstelle von x^(1/3)-y berechnen. Einen guten
>
> Startwert kann man mit einem Polynom berechnen.

Da muss ich Dir recht geben. Ich sehe das auch so. Gut gemacht.

von Yalu X. (yalu) (Moderator)


Bewertung
1 lesenswert
nicht lesenswert
Ich habe hier mal einen Algorithmus zusammengebastelt, der aus einer
ganzen Zahl x den ganzzahligen Anteil von x**(1/3) berechnet. Er ver-
wendet nur Additionen und Schiebeoperation mit konstanter Schiebeweite
<= 3, so dass er gut in Hardware implementierbar sein sollte. Da ich
kein VHDL kann und in Verilog nicht besonders fit bin, habe ich den
Algorithmus in Python aufgeschrieben:
n = 16

def cbrt(x):
  y3 = y2 = y1 = 0
  y0 = 1 << 3*(n-1)
  while True:
    s = y3 + (y2+y1 << 1) + y2 + y1 + y0
    if s <= x:
      y3 = s
      y2 += (y1<<1) + y0
      y1 += y0
    y0 >>= 3
    if y0 == 0:
      return y1
    y1 >>= 2
    y2 >>= 1

print(cbrt(280137072301568))
# Ergebnis: 65432

n ist ein konstanter Parameter und gibt die Bitbreite des Ergebnisses
an. Der Radikand x hat somit 3·n Bits. y0, y1, y2 und y3 sind Register,
die ebenfalls 3·n Bits breit sind. Die Schleife wird n-mal durchlaufen.
Ordnet man die Rechenoperationen innerhalb der Schleife geeignet an, so
dass pro Durchlauf die vier Register genau einmal beschrieben werden,
sollte das Ergebnis also in n Taktzyklen verfügbar sein.

Der Algorithmus dürfte in ähnlicher Form auch für Festkommazahlen imple-
mentierbar sein.

Edit: Für die Register y0, y1 und y2 genügen 3·n-2 Bits.

von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
Interessantes Verfahren. Kannst Du noch einen Satz dazu sagen, woher das 
stammt und worauf das beruht?

Ansonsten liefert Newton für die Kubikwurzel die Iteration
x(n+1)=2/3*x(n)+S/(3*x(n)^2)

Das benötigt ne integer-Division, das will man nicht.

Cheers
Detlef

von Uwe Bonnes (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Mit aktueller Hardware wie Cortex M4 braucht man die Divisionen auch 
nicht mehr so zu fuerchten...

von BiBi (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Uwe Bonnes schrieb:
> Mit aktueller Hardware wie Cortex M4
Der cortex ist aber keine Lösung für das Thema (VHDL) - oder meinst Du 
mit dem Vorschlag einen Core im FPGA?

von Yalu X. (yalu) (Moderator)


Bewertung
0 lesenswert
nicht lesenswert
Detlef _a schrieb:
> Interessantes Verfahren. Kannst Du noch einen Satz dazu sagen, woher das
> stammt

Das habe ich mir selber ausgedacht, ist aber garantiert nicht neu.

> und worauf das beruht?

Die Grundidee ist ganz einfach: Im Ergebniswort werden — beginnend mit
dem höchstwertigen — nacheinander die einzelnen Bits testweise auf 1
gesetzt und geprüft, ob die dritte Potenz des Ergebnisses kleiner oder
gleich dem Radikanden ist. Ist dies der Fall, wird das 1-Bit in das
Ergebnis übernommen, ansonsten bleibt das Bit 0. So wird das Ergebnis
schrittweise verbessert, bis nach dem Bearbeiten des letzten Bits der
Fehler kleiner als 1 ist.

Für die Berechnung der dritten Potenz braucht man normalerweise 2
Multiplikationen. Wenn man schnelle Multiplizierer zur Verfügung hat,
kann man die Potenz direkt berechnen, sonst geht es auch mit Additionen
und Schiebeoperationen, wenn man immer die zuletzt berechnete Potenz
mitbenutzt. Folgende Variablen werden dazu verwendet:

  y  = aktuelles Ergebnis
  z  = die dem aktuell gesetzten Bit entsprechende Zweierpotenz
  y₃ = y³
  y₂ = y²z
  y₁ = yz²
  y₀ = z³

Dabei tauchen y und z im Programm nicht explizit auf.

Wird nun ein zusätzliche Bit im Ergebnis y gesetzt, ist das neue Ergeb-
nis y+z. Die dritte Potenz davon ist

  s = (y + z)³ = y³ + 3(y²z + yz²) + z³ = y₃ + 3(y₂ + y₁) + y₀

Das ist die Formel in der ersten Zeile der While-Schleife. Die Werte y₃,
y₂, y₁ und y₀ stammen aus dem letzten Schleifendurchlauf bzw. der Initi-
alisierung vor der Schleife.

Ist s>x, d.h. das Ergebnis ist zu groß, wird das alte y beibehalten,
ansonsten wird z zu y addiert und die davon abhängigen Variablen wie
folgt aktualisiert:

  y₃ := (y + z)³  = s
  y₂ := (y + z)²z = y²z + 2yz² + z³ = y₂ + 2y₁ + y₀
  y₁ := (y + z)z² = yz² + z³ = y₁ + y₀


Das gleiche Spiel setzt sich nun mit dem nächsten Bit fort. Die Zweier-
potenz z wird also durch 2 dividiert und die Variablen abhängig davon,
wie oft der Faktor z darin enthalten ist, ebenfalls dividiert:

  y₀ := y₀ / 2³
  y₁ := y₁ / 2²
  y₂ := y₂ / 2

Bevor y₁ und y₂ dividiert werden, wird aber noch die Abbruchbedingung
z=0 bzw. y₀=0 überprüft. Ist diese erfüllt, sind alle Bits bearbeitet
worden, und das Gesamtergebnis steht in y₁.

Im Programm sind die Multiplikationen mit 2 und 3 und die Divisionen
durch 2, 4 und 8 mit Schiebeoperationen realisiert.

von MJF (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Hallo BiBi,

es gibt eine einfache Möglichkeit, die dritte Wurzel zu berechnen:

Zunächst eine Annahme, damit es sich leichter erklärt:

y = x^1/3   mit 0<=x<1

Der Eingangswert x kann auch wie folgt beschrieben werden:

x = 2^(-3*n) * x'   mit  1/8 <= x' < 1  und  n \in IN

In HW kann "n" leicht ermittelt werden, indem die Anzahl der führenden 
Nullen detektiert und "weggeshiftet" werden (Multiplexer). Die Anzahl 
der
entfernten führenden Nullen muss ein Vielfaches von 3 sein.

Für die Wurzel gilt dann:

y = x^1/3 = ( 2^(-3*n)*x' )^1/3 = 2^-n * x'^1/3

Nun muss die Wurzel nur noch für einen Eingangswertebereich zwischen 1/8 
und 1 berechnet werden. in diesem Bereich ist die Wurzel recht gutmütig,
d. h. sie kann mit einer liniearen Interpolation sehr gut genähert 
werden.
Für die Lineare Interpolation benötigt man ein RAM für Stützstellen, ein
RAM für die Steigungen einen Multiplizierer und einen Addierer.

Die Multiplikation mit 2^-n ist wieder eine Shift-Operation 
(Multiplexer).

Durch das fehlerfreie Shiften entstehen Fehler nur durch die Lineare
Interpolation. Dieser Fehler tritt immer nur relativ zum Eingangwert
auf, d.h.

\Delta y / y

ist ungefähr konstant.

Man kann die Dimension des Fehlers auch grob abschätzen. Bei einer 
Linearen Interpolation ist Fehler bei N Stützstellen ungefähr 1/N^2
(grobe Daumenregel).

Verwendest du z. B. 1024 Stützstellen, erreichst du einen relativen 
Fehler
von ca. 1/1.000.000.

Um die tatsächlichen Fehler zu bestimmen, muss man schon etwas länger
rechnen.

Dies ist keine iterative Methode, d. h. mit jedem Takt kann ein 
Ausgangswert berechnet werden.

Man kann es mit etwas mehr Aufwand auch für beliebige Wurzeln erweitern.

Gruss

Markus

von Jürgen S. (engineer) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
Ich komme mit dem Python Konstrukt nicht ganz klar:
s = y3 + (y2+y1 << 1) + y2 + y1 + y0

Bezieht sich der Schiebeoperator (y2+y1 << 1) auf y2 und Y1 oder nur Y1?

Ich habe mal etwas gekramt: Mit Newton geht es in der Tat recht gut, 
wenn man einen guten Startwert. Den kann man zur Basis 2 einfach 
bestimmen, wenn man die Bits etwas schiebt. Das höchstwertige einfach 
auf ein Drittel der Position, den Rest als Multiplikator oben drauf. 
Dann kann man für ein eventuell existentes Bit eines darunter noch 
feststellen, dass der Wert im Bereich von 50% liegen muss und daher die 
dritte Wurzel aus 1,5 aufmultiplizieren -> rund 1+1/8

von Yalu X. (yalu) (Moderator)


Bewertung
0 lesenswert
nicht lesenswert
J. S. schrieb:
> Ich komme mit dem Python Konstrukt nicht ganz klar:
> s = y3 + (y2+y1 << 1) + y2 + y1 + y0

Der Vorrang der Operatoren ist gleich wie in C:

s = y3 + ((y2+y1) << 1) + y2 + y1 + y0

Das ist wiederum das Gleiche wie

s = y3 + 3 * (y2 + y1) + y0

Ich wollte mit der Aktion nur die Multiplikation wegoptimieren. Es kann
aber gut sein, dass das der VHDL-/Verilog-Compiler schon tut.

von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
Hallo,

ich habe Routinen für Quadrat- und Kubikwurzel für Integerrechnung 
implementiert und in die Codesammlung eingestellt.

Beitrag "square root and cubic root for integer and FPGA implementation"

Diese Routinen benötigen nur Shift und Addition und eignen sich damit 
gut für die Implementierung in Hardware.

Yalu, Dein Verfahren habe ich in ähnlicher Form für Quadratwurzeln bei
http://en.wikipedia.org/wiki/Methods_of_computing_square_roots
und
http://www.cc.utah.edu/~nahaj/factoring/isqrt.legalize.c.html
gefunden.

Für Kubikwurzeln habe ich es in C implementiert und auch in eine Form 
modifiziert, die mehr Richtung der genannten Referenzen geht: mit einem 
remainder anstelle der eigentlichen Werte.

War großer Spaß gewesen, ordentlich was bei gelernt, danke für die 
Anregung. Die Idee, was 'Bit für Bit' zu iterieren hat auch woanders 
Potential, glaube ich.

Cheers
Detlef

von Jan (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Hast Du das so gedacht, dass dann in HW eine FSM laufen soll oder ein 
soft core soll das so in der Synthese aufgezogen werden?

von Detlef _. (detlef_a)


Bewertung
1 lesenswert
nicht lesenswert
Das kann man im FPGA auf verschiedene Weisen ausführen, ich bin da nicht 
besonders kundig. Kubikwurzel dauert für 30Bit Zahl 10 Runden, schneller 
kann man nicht werden. Ansonsten läßt sich die Verarbeitung beliebig 
serialisieren, parallelisieren und pipelinen.

Cheers
Detlef

von Mohan (Gast)


Bewertung
-2 lesenswert
nicht lesenswert
Kannst Du mir helfen mit Matlab Simulink schematic fur die dritten 
wurzel in  FPGA, sogar mit einem Beischpiel zu bilden ?

von Jürgen S. (engineer) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
Ich glaube, an der Stelle ist Simulink Kanonen auf Spatzen. Das geht mit 
einer einfachen Ausgleichsrechnung oder händischen Potenz mit "hoch 
1/3", im Zweifel über Logarithmus und CORDIC. Kommt drauf an, ob es 
schnell, schnell programmiert oder einfach oder einfach programmiert 
oder klein sein soll. Am kleinsten und Stromsparendsten ist sicher ein 
iterativer Schätzer, mit Newton und der ersten Ableitung.

von frk (Gast)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Ich habe heute mal fix das ganze in VHDL gehackt. Es wird wohl keinen 
Schönheitspreis gewinnen. Auf dem FPGA habe ich es nicht getestet und 
auch in ModelSim nur ganz rudimentär getestet.

von Jürgen S. (engineer) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
Da müsste noch pipelining rein, meine ich.

von Markus F. (mfro)


Bewertung
0 lesenswert
nicht lesenswert
Jürgen S. schrieb:
> Da müsste noch pipelining rein, meine ich.

Ist doch schon drin?

von frk (Gast)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Jürgen S. schrieb:
> Da müsste noch pipelining rein, meine ich.

Ja, das sollte gehen. Bzw.: Ja, das geht (siehe Anhang).

Ich hatte noch ein Problem mit dem cubedbit, weil Integer in VHDL auf 32 
Bit begrenzt sind und ich das Bit als Integer berechne bevor ich es in 
ein unsigned umwandele. Ich habe das jetzt nachgebessert und noch eine 
Rundung des Ergebnisses hinzugefügt.

Bitte entschuldigt, wenn ich Fragen nur mit großer Verzögerung 
beantworte. Ich bin aktuell sehr beschäftigt und wollte eigentlich nur 
schnell den Code hier zur Verfügung stellen, den ich sowieso schon für 
mich geschrieben habe, um Leuten zu helfen, denen VHDL nicht so schnell 
von der Hand geht.

von Jürgen S. (engineer) Benutzerseite


Bewertung
0 lesenswert
nicht lesenswert
Markus F. schrieb:
> Jürgen S. schrieb:
>> Da müsste noch pipelining rein, meine ich.
>
> Ist doch schon drin?
Nicht so, wie ich mir das denke :-)

frk schrieb:
> Ja, das sollte gehen.
Die Art des Rundens würde ich mir nochmal ansehen.

Warum wird das überhaupt gemacht? Ich würde noch die Nachkommastellen 
prozessieren.

: Bearbeitet durch User

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]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [vhdl]VHDL-Code[/vhdl]
  • [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.

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