Benutzer-Werkzeuge

Webseiten-Werkzeuge


examples:square_root

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

<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

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

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
examples/square_root.txt · Zuletzt geändert: 2013-06-06 21:26 von 127.0.0.1