Das Quadratwurzel zu ziehen ist ein Standard-Problem. Es gibt verschiedene Verfahren dafür.
Hierbei wird mit dem Ausdruck
<latex>x' = {{n\over x} + x\over 2}</latex>
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
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 ( 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 ( n -- x ) dup 0< IF -24 throw THEN dup BEGIN 2dup / over + 2/ swap over - abs 1 <= UNTIL nip ;
: sqrt-3 ( u1 -- u2 ) dup 0< IF -24 throw THEN dup begin >r dup r@ / r@ + 2/ dup r> - abs 2 < until nip ;
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 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.
: 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 ;
\ 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/ ;
( 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 ;
Für die ursprüngliche Anlage, hilfreiche Kritik und weitere Entwicklung meinen Dank an: Bernd, dmichelsen, uho, mhx, bushmills.
usenet: