Quadratwurzel ziehen mittels Assembler

[Dr. Irata]

Im Rahmen des Ray Casting bin ich auf eine weitere Hürde gestossen: Zur Berechnung der Distanz zwischen Player und Wand/Hürde brauche ich die Quadratwurzel.
Mehrere Möglichkeiten fallen mir nach langer Überlegung dazu ein: Tabellen, bei einem großen Playfield wird das aber recht speicherintensiv... oder berechnen! Die mir bekannte schriftliche Berechnung ist meines Erachtens schwierig umzusetzen, daher habe ich weiter recherchiert und bin auf eine andere einfachere Lösung gestossen.... es geht um die Heron-Methode!

Als Beispiel nehme ich mal die Wurzel aus 12 - laut Taschenrechner = 3,464
Man definiert ein Rechteck mit 1x12 und nähert sich dann dem Quadrat an.
Zunächst wird ein x0 gebildet: x0 = (12+1)/2 = 6,5
Dann x1: x1 = (6,5+12/6,5) /2 = 4,173
Dann x2: x2 = (4,173 + 12/4,173) / 2 = 3,52
Dann x3: x3 = (3,52 + 12/3,52) / 2 = 3,465

Mit jedem Iterationsschritt wird das Ergebnis genauer.

Wenn mir das Ergebnis etwas gröber reicht und ich die Kommastellen weglasse kommt folgendes raus:

x0 = 6,5 = 6
x1 = 4
x2 = 3,5 = 3
x3 = 3,5 = 3

Edit: Bei hohen Zahlen sollte man eine Iteration mehr nehmen!!
So liefert die Methode bei Wurzel aus 134 (=11,5758) nach 3 Schritten mit Abrundung das Ergebnis 13 und nach einem Schritt mehr kommt man auf 11,5 also abgerundet auf 11 und aufgerundet auf 12!


Als nächstes gilt es diesen Heronalgorithmus in Assembler umzusetzen....

Gibt es noch einfachere Möglichkeiten oder Umsetzungen????

[pps]

Hier gibt es eine sehr gute Übersicht mitsamt Performance Test für 6502 Assembler und 16 bit Integer SQRT: https://github.com/TobyLobster/sqrt_test

Eine andere Idee wäre, die Wurzel durch "Quadratisches Ausgleichen" der Gleichung zu entfernen. Vielleicht ist das am Ende schneller.

[slx]

CORDIC?

[Dr. Irata]

CORDIC kannte ich noch nicht - sehr interessant / vielen Dank - scheint aber eher langsam zu sein und setzt doch ein wenig mathematisches Verständnis voraus (Vektorrotation etc.) - das dann noch in Assembler umzusetzen ist sicherlich sportlich.

[Mathy]

.
Hallo Peter

Vielleicht kann Dietrich dir helfen, er ist Mathematiker. Und kommt auch zur Fujiama.

Tschüß

Mathy

[Dr. Irata]

Gute Idee - da werde ich ihn unbedingt mal drauf ansprechen und ne kleine Mathestunde von ihm erbeten .. bin ja selber so eine Art Hobbymathematiker

[Irgendwer]

Meine Routine aus "Disc o' Pop":

Code: Alles auswählen

; UWORD __fastcall__
; SM_Sqrt(ULONG lValue_);
; 
.define result sreg
.define input  ptr4
.define T      ptr1

.export _SM_Sqrt
.proc 	_SM_Sqrt

		sta input
		stx input+1
		lda sreg
		sta input+2
		lda sreg+1
		sta input+3

		lda #0
        sta result        ; R=0
        sta result+1
        sta input+4
        ;sta T+15        ; (T+1) is zero until last iteration; (T+0) is always 0
 
        clc
        ldy #14          ; 15 iterations (14-->0) + last iteration
loopsq:  
        lda result       ; (2*R+D) LSR 1; actually: R+(D LSR 1)
        ora stablo,y     ; using ORA instead of ADC is ok because the bit to be set
        sta T+2          ; will have not been affected yet
        lda result+1
        ora stabhi,y
        sta T+3
        bcs skip0      	 ; takes care of large numbers; if set, input>T
 
        lda input+3
        cmp T+3
        bcc skip1        ; T <= input    (branch if T>input)
        bne skip0
        lda input+2
        cmp T+2
        bcc skip1
skip0:  ;sec
        lda input+2      ; input=input-T
        sbc T+2
        sta input+2
        lda input+3
        sbc T+3
        sta input+3
        lda result       ; R=R+D
        ora stablo+1,y
        sta result
        lda result+1
        ora stabhi+1,y
        sta result+1
skip1:
        asl input       ; input=input*2
        rol input+1 
        rol input+2
        rol input+3
        dey             ; implicit: D=D/2, by the decrement of .Y
        bpl loopsq
lastiter:               ; code for last iteration
        bcs skp0        ; takes care of large numbers; if set, input>T
        				; during last iteration D=1, so [(2*R+D) LSR 1] makes D the MSB of T+1
        lda input+3
        cmp result+1    ; (T+3) = resultHI
        bcc skp1        ; T <= input    branch if T>input
        bne skp0
        lda input+2
        cmp result      ; (T+2) = resultLO
        bcc skp1
        bne skp0
        lda input+1
        cmp #$80        ; value of (T+1) during last iteration
        bcc skp1
skp0:
	    ;sec
        lda input+1
        sbc #$80        ; (T+1) during last iteration
        sta input+1
        lda input+2     ; input=input-T
        sbc result      ; (T+2)
        sta input+2
        lda input+3
        sbc result+1    ; (T+3)
        sta input+3
        inc result      ; R=R+D with D=1
skp1:   
        lda result
        ldx result+1
        rts

.endproc

.data 
stabhi: .byte 0,0,0,0,0,0,0,0
stablo: .byte $01,$02,$04,$08,$10,$20,$40,$80
        .byte 0,0,0,0,0,0,0,0

[Dr. Irata]

sehr cool ... welchem Algorithmus liegt dem zugrunde??

[LarsImNetz]

Hallo Peter,
vielleicht mal zur Definition: Welchen Zahlentyp verwendest Du für den Abstand?

Der Abstand lässt sich auch prima berechnen mit:

Abstand := abs(x1-x2) + abs(y1-y2)

Das habe ich bei Pacmen so gemacht, reicht für den groben Abstand völlig aus.

Ansonsten nimm eine der Versionen von Toby.

LG
Lars

[Dr. Irata]

Hallo Lars,
das geht natürlich im Prinzip auch - und ich denke bei Pacman ist das super schnell und sinnvoll - die Näherung (x1-x2) + (y1-y2)/2 wenn y1-y2 die kurze Seite ist, ist deutlich genauer. Ich denke beim Raytracing sollte es etwas genauer sein, da vielleicht sonst eine Verzerrung auftritt.... ich weiß es aber noch nicht.
Ich werde beides testen und dann posten!
LG
Peter

[Dr. Irata]

Meine Routine aus "Disc o' Pop":

Code: Alles auswählen

; UWORD __fastcall__
; SM_Sqrt(ULONG lValue_);
; 
.define result sreg
.define input  ptr4
.define T      ptr1

.export _SM_Sqrt
.proc 	_SM_Sqrt

		sta input
		stx input+1
		lda sreg
		sta input+2
		lda sreg+1
		sta input+3

		lda #0
        sta result        ; R=0
        sta result+1
        sta input+4
        ;sta T+15        ; (T+1) is zero until last iteration; (T+0) is always 0
 
        clc
        ldy #14          ; 15 iterations (14-->0) + last iteration
loopsq:  
        lda result       ; (2*R+D) LSR 1; actually: R+(D LSR 1)
        ora stablo,y     ; using ORA instead of ADC is ok because the bit to be set
        sta T+2          ; will have not been affected yet
        lda result+1
        ora stabhi,y
        sta T+3
        bcs skip0      	 ; takes care of large numbers; if set, input>T
 
        lda input+3
        cmp T+3
        bcc skip1        ; T <= input    (branch if T>input)
        bne skip0
        lda input+2
        cmp T+2
        bcc skip1
skip0:  ;sec
        lda input+2      ; input=input-T
        sbc T+2
        sta input+2
        lda input+3
        sbc T+3
        sta input+3
        lda result       ; R=R+D
        ora stablo+1,y
        sta result
        lda result+1
        ora stabhi+1,y
        sta result+1
skip1:
        asl input       ; input=input*2
        rol input+1 
        rol input+2
        rol input+3
        dey             ; implicit: D=D/2, by the decrement of .Y
        bpl loopsq
lastiter:               ; code for last iteration
        bcs skp0        ; takes care of large numbers; if set, input>T
        				; during last iteration D=1, so [(2*R+D) LSR 1] makes D the MSB of T+1
        lda input+3
        cmp result+1    ; (T+3) = resultHI
        bcc skp1        ; T <= input    branch if T>input
        bne skp0
        lda input+2
        cmp result      ; (T+2) = resultLO
        bcc skp1
        bne skp0
        lda input+1
        cmp #$80        ; value of (T+1) during last iteration
        bcc skp1
skp0:
	    ;sec
        lda input+1
        sbc #$80        ; (T+1) during last iteration
        sta input+1
        lda input+2     ; input=input-T
        sbc result      ; (T+2)
        sta input+2
        lda input+3
        sbc result+1    ; (T+3)
        sta input+3
        inc result      ; R=R+D with D=1
skp1:   
        lda result
        ldx result+1
        rts

.endproc

.data 
stabhi: .byte 0,0,0,0,0,0,0,0
stablo: .byte $01,$02,$04,$08,$10,$20,$40,$80
        .byte 0,0,0,0,0,0,0,0

welcher Idee bzw. Routine/Algorithmus steckt denn dahinter???

[LarsImNetz]

Ich glaube:
https://de.wikipedia.org/wiki/Heron-Verfahren

In Java könnte man es wie folgt machen:

Code: Alles auswählen

public int sqrt(int WurzelAus) {
  int Wert;
  int nextWert = (WurzelAus + 1) / 2;
  do {
    Wert = nextWert;
    nextWert = (Wert + (WurzelAus / Wert)) / 2;
  } while (nextWert != Wert);
  return Wert;
}

[Dr. Irata]

... ja Heron ... leider englisch und nicht super kommentiert...

Ich werde hier mal für MADS eine Heron-Routine machen, die wir dann auch ins Wiki stellen können...

[Dr. Irata]

.... ich habe das Wurzel ziehen mittels Heron-Iteration zunächst mal in Basic umgesetzt und geschaut, wie genau das Ergebnis nach x-Iterationsschritten aussieht.
Wenn man im 8-Bit Bereich bleibt, dann reichen 6-7 Schritte völlig aus für ein sehr genaues Ergebnis. Geht man in den 16-Bit Bereich und will z.B. von 65000 die Quadratwurzel ziehen, dann braucht man 10 Durchläufe.

Hier zunächst der Basic-Code:

Code: Alles auswählen

5	Graphics 0
10	Input A
20	B=1
30	C=A
40	For F=1 To 6
50	C=C+B
60	C=C/2
70	B=A/C
80	Next F
90	Print C

Mit B=1 und C=A definiere ich das erste Rechteck C x B (9x1) - nach Heron wird nun Schritt für Schritt die Seite C verkleinert und B vergrößert, bis beide Seiten sich annähern. Wenn beide Seiten annähernd gleich sind, dann steht das Ergebnis fest - die eine Seite (hier C) hat dann das Ergebnis Wurzel von A!
Schon schlau die alten Griechen....

Bevor ich nun den Code in Assembler umsetzen, brauche ich noch eine Routine fürs Teilen!

[Dr. Irata]

... und hier jetzt der Assemblercode für die Quadratwurzel im Wertebereich 0-255 (8 Bit):

Wurzel aus var1 -> Ergebnis in var3

Code: Alles auswählen

			 .proc qwurzel	;von var1
			 ldx #0
			 ldy #6
			 lda var1
			 sta var3
			 lda #1
			 sta var2
	lp1
			clc
			lda var3
			adc var2
			sta var3
			lsr var3
			lda var1
			;------------ hier Divisionsroutine
			ldx #255
			sec
	lp2		
			sbc var3
			inx
			bcs lp2
			stx var2
			;----------- Division Ende
			dey
			bpl lp1
			rts
			.endp

[pps]

Y sollte aber anfangs definiert werden... Oder habe ich was übersehen?

Und dann sinnvoller rückwärts... Von 6 bis negativ...

[Dr. Irata]

y muss auf Null sein - stimmt
y rückwärts macht Sinn und ist schneller!
Habe es gerade oben geändert!! Tak!!

[Dr. Irata]

als nächster Schritt wäre jetzt noch interessant, die Quadratwurzel auf eine 16 Bit- Zahl zu erweitern, dann bebommt man ruck zuck auf noch Pythagoras... Freiwillige vor

[pps]

als nächster Schritt wäre jetzt noch interessant, die Quadratwurzel auf eine 16 Bit- Zahl zu erweitern, dann bebommt man ruck zuck auf noch Pythagoras... Freiwillige vor

Aus meinem verlinkten Test mehrerer Routinen finde ich am Ende das sehr interessant (sqrt7). Hat nur 38 Bytes und ist dennoch eine der schnellsten Routinen insgesamt:

http://6502org.wikidot.com/software-math-sqrt

Code: Alles auswählen

   LDA #0
   STA ROOT
   STA REM
   LDX #8
L1 SEC
   LDA NUMH
   SBC #$40
   TAY
   LDA REM
   SBC ROOT
   BCC L2
   STY NUMH
   STA REM
L2 ROL ROOT
   ASL NUML
   ROL NUMH
   ROL REM
   ASL NUML
   ROL NUMH
   ROL REM
   DEX
   BNE L1
 

[Dr. Irata]

tolle Routine!
Ich finde es in der Mathematik und in der Informatik immer wieder faszinierend, wie viele verschiedene Ansätze zur gleichen Lösung führen!!

[Mathy]

.
Hallo Peter

... dann bebommt man ruck zuck auf noch Pythagoras...

Pita Gyros! hmmmmm....

Tschüß

Mathy

[Dr. Irata]

Ich habe jetzt die Routine für die Quadratwurzel einer 16-Bit Zahl nochmal analysiert und getestet und entsprechend der Variablen angepasst. Diese Routine kommt auch ins Atari-Wiki und ist hinsichtlich Schnelligkeit und Kreativität ein echter Schatz!!!
In Gedenken an Lee Davison, der leider schon früh gestorben ist....

In der Routine wird die Quadratwurzel von var1/var1+1 gebildet. Das Ergebnis ist dann in var2, der Rest in var3.
Ich habe versucht es zu durchdringen - es ist mir nicht vollständig gelungen! Insgesamt wurde hier ein schriftliches Wurzelziehen auf 2-Bit Ebene umgesetzt (daher die vielen rol und asl Anweisungen). Die Routine erscheint wirklich stark optimiert und ist entsprechend auch schön kompakt und schnell! Ich denke Lee hat da lange dran getüftelt und hat meinen uneingeschränkten Respekt dafür!!

Code: Alles auswählen

		.proc qwurzel16bit
		; Quadratwurzel 16 Bit
		; nach Lee Davison
		; Quadratwurzel var1/var1+1
		
		lda #0
		sta var2
		sta var3
		ldx #8
lp1
		sec
		lda var1+1
		sbc #64
		tay
		lda var3
		sbc var2
		bcc lp2
		sty var1+1
		sta var3
lp2
		rol var2
		asl var1
		rol var1+1
		rol var3
		asl var1
		rol var1+1
		rol var3
		dex
		bne lp1
		rts
		.endp

[Dr. Irata]

.... vielleicht ist hier auch der CORDIC Algorithmus verbaut .... mit diesem werde ich mich auch noch demnächst etwas näher beschäftigen!