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/ ;

• 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 ;

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