====== Quadratwurzel ziehen ======
Das Quadratwurzel zu ziehen ist ein Standard-Problem.
Es gibt verschiedene Verfahren dafür.
===== signed integer square root (Heron) =====
Hierbei wird mit dem Ausdruck
x' = {{n\over x} + x\over 2}
also Näherungswert=(Schätzwert+Eingangswert/Schätzwert)/2 so lange an das Ergebnis angenähert, bis eine vorher festgelegte Genauigkeit erreicht worden ist. Das Verfahren konvergiert rasch. Algebraisch formuliert lautet die Vorschrift:
xn+1 = ( xn + a/xn ) / 2
==== Beispiel ====
a=25, x0=a ; Integer rechnen.
x1 = (25+25/25)/2 = (25+1)/2 = 13
x2 = (13+25/13)/2 = (13+1)/2 = 7
x3 = (7+25/7)/2 = (7+3)/2 = 5
x4 = (5+25/5)/2 = (5+5)/2 = 5
x5 = (5+25/5)/2 = (5+5)/2 = 5 ...
Der Wert beleibt nun auf 5 stehen, dh, die Näherung ist abgeschlossen, wenn xn+1=xn wird. Anders ausgedrückt, Schluß ist wenn gilt:
| xn+1 - xn | <= 0
==== sqrt ====
: sqrt ( n -- x )
dup BEGIN 2dup / over + 2/ swap over - abs 1 <= UNTIL nip ;
Um diesen Algorithmus zu entwickeln setzen wir die Näherungsformel zunächst in RPN um.
a xn / xn + 2/ ==> xn+1
Dann erfolgt die Umsetzung in Forth. Dazu erstellen wir ein Stackbild der Operationen. Ausgangspunkt sind der Anfangswert n und der Schätrzwert x0.
( -- n x0 )
Diese müssen auf dem Stack für die Berechnung in der Schleife erhalten werden. Um daraus den Bruch n/xo zubilden, wird also zunächst verdoppelt und dann geteilt.
2dup ( n x0 -- n x0 n x0 )
/ ( -- n x0 a )
Nun kann x0 dazu addiert werden, anschließend wird durch 2 geteilt um den nächsten Schätzwert zu erhalten.
over ( -- n x0 a x0 )
+ ( -- n x0 b )
2/ ( -- n x0 x1 )
Der neue Schätzwert ersetzt nun den Alten.
swap ( -- n x1 x0 )
Damit kann nun geprüft werden, ob die Annäherung an das Ergebnis schon erreicht worden ist.
over ( -- n x1 x0 x1 )
- abs 1 <= ( -- n x1 flag )
Damit wird die Schleife solange durchlaufen, bis flag=true geworden ist, und wir haben x1=x als Ergebnis. Der Eingangswert n wird fallen gelassen.
( -- n x1 )
swap drop ( -- x )
Der Term swap drop ist, weil oft gebraucht, in manchen Forthsystemen vorgefertigt und wird nip genannt.
: nip swap drop ;
Achtung
In diesem Verfahren führen negative Eingangswerte zu einer Endlosschleife. Je nachdem welches Verhalten man für diesen Falle haben möchte, kann der Eingangswert zunächst abgefragt werden. So käme -24 throw in Frage: invalid numeric argument.
==== sqrt-2 ====
: sqrt-2 ( n -- x )
dup 0< IF -24 throw THEN
dup BEGIN 2dup / over + 2/ swap over - abs 1 <= UNTIL nip ;
==== sqrt-3 ====
: sqrt-3 ( u1 -- u2 )
dup 0< IF -24 throw THEN
dup begin >r dup r@ / r@ + 2/ dup r> - abs 2 < until nip ;
===== unsigned integer square root (Wil Baden) =====
Soll hingegen ohne Vorzeichen gerechnet werden, wird gern das folgende Verfahren genommen, weil es so einfach ist. Es hat allerdings den Nachteil für größere Eingangswerte relativ langsam zu sein, da n-te wurzel + 1 Näherungsschritte erforderlich sind.
[[proof of sqrt baden|proof of square root (W.Baden)]]
\ Paul E Bennet comp.arch.embedded 4 May 2008
\ Wil Baden; comp.lang.forth; "Square Root the Slow Way"
\ square root unsignd integer
: sqrt -1 swap over do 2 + dup +loop 2/ ;
"It is a standard method of calculating square root by summing
the odd integers. It works for all unsigned integers from
0 to FFFFFFFF (or FFFF).
I discovered the Forth code when implementing the sum of odd
integers algorithm. It is neat how the Forth primitives work
to implement this algorithm."
zb gforth (cell=4; 32bit Sytem):
-1 sqrt-pb u. 65535 ok
Hierbei wird also vorzeichenlos eine Näherung ermittelt.
===== Weitere Verfahren =====
==== isqrt und sqrt-tab (bushmils) ====
: ISQRT ( n -- isqrt )
DUP 0= ?EXIT
DUP >R
BEGIN R@ OVER / OVER + 2/
DUP ROT - ABS 1 <=
UNTIL R> DROP ;
Das "Problem" hierbei ist, dass es 52 Stellen gibt, wo der ISQRT einen anderen Wert liefert als FSQRT FROUND F>S .
Mit etwas mehr data space geht es schneller, wenn man z.B. eine 64K Tabelle anlegt:
\ gforth compatibility for testing:
: -R postpone r> postpone drop ; immediate
256 dup * constant #256 ( 64K)
CREATE sqrdtab #256 CELLS ALLOT
MARKER -nonce :NONAME #256 0 DO I I * sqrdtab I CELLS + ! LOOP ; EXECUTE -nonce
: SQRT-TAB ( x -- sqrt[x] )
>R #256 -1 \ initialize lower and upper limits
BEGIN 2DUP - 1 > \ not yet done?
WHILE 2DUP + 2/ \ compute midpoint ( hi low mid )
DUP CELLS sqrdtab + @
R@ U<= IF NIP \ replace lower limit
ELSE ROT DROP SWAP \ replace upper limit
ENDIF ( cr .s )
REPEAT -R NIP
DUP -1 = IF 1+ EXIT ENDIF ;
==== 32bit sqrt auf 16bit Maschine (Klaus Kohl) ====
\ Quelle: 4d1992_3.pdf; Fouriertransformation, Klaus Kohl. Screen#5 23.05.89/KK
( SQRT)
( Quadratwurzel einer 32-Bit-Integerzahl)
( Zulässiger Bereich: $3FFF.0001 -> 0..32767 ( 15 Bit )
( : under swap over ; )
: u2/ ( u - u/2 ) 2/ $7FFF and ; ( lower 15 bit erhalten )
: d2* ( d - d*2 ) over 2* -rot 2* swap 0< - ; ( Bit 15 übertragen )
: d2/ ( d - d/2 ) under 2/ swap u2/ rot 1 and IF $8000 or THEN swap ;
: d- dnegate d+ ;
: SQRT-kk ( ud - u )
0 0
15
FOR >r 2* over 0< - 2* over 2* 0< - ( 2 bits in S )
>r d2* d2* r> r> ( Quadrat*4)
2* 1+ ( Wert*2+1)
2dup u< ( Summe kleiner Wurzel?)
IF 1- ( dann nur Wert*2 )
ELSE dup >r - r> 1+ ( oder Summe=Differenz)
THEN
NEXT
nip nip nip U2/ ;
===== Aufgaben =====
* Welcher Algorithmus steckt jeweils dahinter?
* Wieso terminiert die Schleife immer? Gibt es kein Oszillieren?
* Funktioniert das mit vorzeichenlosen und vorzeichenbehafteten Zahlen?
* Wie sieht es mit dem Zeitverhalten aus? Wie oft wird die Schleife durchlaufen?
( gforth, 32bit System)
: test ( von bis -- )
swap DO cr i .
i sqrt .
i s>d d>f fsqrt fdup FROUND F>d d>s . f.
LOOP ;
===== Dank & Links =====
Für die ursprüngliche Anlage, hilfreiche Kritik und weitere Entwicklung meinen Dank an: Bernd, dmichelsen, uho, mhx, bushmills.
usenet:
* comp.arch.embedded
* comp.lang.forth