erweiterte Näherung zu Pythagoras

[Dr. Irata]

... der Pythagoras läßt mich einfach nicht los!

In meinen langen Erläuterungen und Versuchen einen schönen Raycaster auf dem A8 zu bauen (ich denke das Projekt werde ich wieder im Herbst aufnehmen), bin ich zwingend abhängig auf die Berechnungen von Entfernungen. In Assembler ist das gar nicht so einfach, denn quadrieren und Wurzel ziehen bedarf besonderer Routinen (die ich ja später hier alle hergeleitet habe). Diese Routinen existierten damals so nicht und ich war auf eine ganz gute Näherung angewiesen, die ich dann auch im Netz finden konnte - in den 50 igern und 60 igern gerade bei den einfachen Rechnern waren diese Näherungen notwendig und einige Mathematiker echt erfinderisch.
Die Näherung für c^2 = a^2+b^2 also c = sqrt (a^2+b^2) ist ja:
c = a+b/2 für a>b
Die Näherung ist echt gut und super schnell. Leider ist die an bestimmten Stellen fehlerhaft. Diesen Fehler kann man aber ganz gut ausgleichen:
Tendenziell ist die Näherung immer zu groß. Wenn man jetzt aber vom Ergebnis c die erste Stelle von a+1 subtrahiert solange b>10 ist, dann ist das Ergebnis wirklich nahezu perfekt.
Beispiel: a=40 , b=30 --> c=40-5+30/2 = 50 (mathematisch 50)
Beispiel: a=35 , b=25 --> c=35-4+25/2 = 43 (mathematisch 43)
Beispiel: a=72, b=62 --> c=72-8+62/2 = 95 (mathematisch 95)

[Dr. Irata]

Also: c = sqrt (a^2+b^2)
wenn a>=b
wenn b>=10
dann c=a-(a/10)+b/2
wenn b<10
dann c=a+b/2

wenn b>a
wenn a>=10
dann c=b-(b/10)+a/2
wenn a<10
dann c=b+a/2

und jetzt eine Aufgabe in die "Gemeinde":

Wer kann das in möglichst schnellen Assembler-Code umsetzten?
Bin gespannt auf eure innovativen Lösungen...

[RhoSigma]

Äh, hab zwar keinen Code für Dich, ABER warum Basis 10? Die division durch 10 ist lahm. Nimm doch 16. Oder 8?

[Dr. Irata]

... alle guten Ideen sind hier willkommen. Ich brauche ja die erste Ziffer der Zahl a - also geteilt durch 10. Das ist blöd in Assembler. Durch 8 teilen passt bestimmt auch noch für die Näherung - oder im Dezimalmodus den vorderen Nibble extrahieren...
Es muss halt schnell sein!
Wenn ich durch 8 teile, dann fällt sogar das +1 weg ... also ist noch besser.
Das funktioniert bis 70 sehr gut. Viel größere Entfernungen werde wir am A8 bei den möglichen Spielen nicht brauchen.
Sehr gut!!

[Kveldulfur]

Hallo!

Code: Alles auswählen

		ORG $2000
		
DIV10TAB
		.BYTE 0,0,0,0,0,0,0,0,0,0
		.BYTE 1,1,1,1,1,1,1,1,1,1
		.BYTE 2,2,2,2,2,2,2,2,2,2
		.BYTE 3,3,3,3,3,3,3,3,3,3
		.BYTE 4,4,4,4,4,4,4,4,4,4
		.BYTE 5,5,5,5,5,5,5,5,5,5
		.BYTE 6,6,6,6,6,6,6,6,6,6
		.BYTE 7,7,7,7,7,7,7,7,7,7
		.BYTE 8,8,8,8,8,8,8,8,8,8
		.BYTE 9,9,9,9,9,9,9,9,9,9
		.BYTE 10,10,10,10,10,10,10,10,10,10
		.BYTE 11,11,11,11,11,11,11,11,11,11
		.BYTE 12,12,12,12,12,12,12,12,12,12
		.BYTE 13,13,13,13,13,13,13,13,13,13
		.BYTE 14,14,14,14,14,14,14,14,14,14
		.BYTE 15,15,15,15,15,15,15,15,15,15
		.BYTE 16,16,16,16,16,16,16,16,16,16
		.BYTE 17,17,17,17,17,17,17,17,17,17
		.BYTE 18,18,18,18,18,18,18,18,18,18
		.BYTE 19,19,19,19,19,19,19,19,19,19
		.BYTE 20,20,20,20,20,20,20,20,20,20
		.BYTE 21,21,21,21,21,21,21,21,21,21
		.BYTE 22,22,22,22,22,22,22,22,22,22
		.BYTE 23,23,23,23,23,23,23,23,23,23
		.BYTE 24,24,24,24,24,24,24,24,24,24
		.BYTE 25,25,25,25,25,25		
		
		
VarA		.BYTE 50
VarB		.BYTE 30
VarC		.BYTE 00

Calc:		LDA VarA		; A>=B?
		CMP VarB
		BCS AgtB
		JMP BgtA
		
AgtB:		LDA VarB		; B>=10?
		CMP #10
		BCS Bgt10
		
		LDA VarB
		LSR			; b/2
		CLC
		ADC VarA
		STA VarC
		
		RTS
		
Bgt10:		SEC
		LDX VarA	
		TXA
		SBC DIV10TAB,X		; a/10
		STA VarC		; VarC = a-(a/10)
		
		LDA VarB		; b/2
		LSR
		
		CLC
		ADC VarC		; a-(a/10)
		STA VarC		; VarC = b/2 + a-(a/10) 
		
		RTS
		
BgtA:		LDA VarA		; a>=10?
		CMP #10
		BCS Agt10
		
		LDA VarA
		LSR			; a/2
		CLC
		ADC VarB
		STA VarC
		
		RTS
		
Agt10:		SEC
		LDX VarB	
		TXA
		SBC DIV10TAB,X		; b/10
		STA VarC		; VarC = b-(b/10)
		
		LDA VarA		; a/2
		LSR
		
		CLC
		ADC VarC		; b-(b/10)
		STA VarC		; VarC = a/2 + b-(b/10) 
		
		RTS
	
Start:		JSR Calc
		BRK
				
		RUN Start

Wenn nur 8 Bit benötigt werden, kann man geteilt 10 am Besten in einer Tabelle realisieren.
Ich hoffe der Code hat keine Flüchtigkeitsfehler

Grüße
Janko

[Dr. Irata]

Hi Janko,
sehr super und danke.
Ich habe das aber nochmal geprüft. Nimm mal nicht durch 10 sondern durch 8. Dann kann man das +1 für den zu subtrahierenden Wert weglassen. Das +1 hatte ich ganz oben erwähnt, in der Beschreibung für den Code aber vergessen.
Wenn man aber statt durch 10 nur durch 8 teilt, dann fällt praktischerweise das +1 weg.

Wieviele Zyklen braucht das dann?

[Dr. Irata]

Janko hat mir gerade geschrieben, daß das so etwa 34-51 Zyklen braucht.
Das ist echt super schnell - und läßt sich sicher noch ein wenig optimieren.
Wer schafft es schneller?

Diese Näherung ist wirklich gut, schwächelt nur etwas bei b<10 ... da werde ich mal schauen, ob man das noch für den Bereich besser schafft.

[RhoSigma]

Code: Alles auswählen

		ORG $2000
		
VarA		.BYTE 50
VarB		.BYTE 30
VarC		.BYTE 00

Calc:
		LDA VarA		; A>=B?
		CMP VarB
		BCS AgtB

BgtA:
		CMP #10		; a>=10?
		BCS Agt10
		
		LSR			; a/2
		CLC
		ADC VarB
		STA VarC
		
		RTS
		
AgtB:
		LDA VarB		; B>=10?
		CMP #10
		BCS Bgt10
		
		LSR			; b/2
		CLC
		ADC VarA
		STA VarC
		
		RTS
		
Bgt10:
		LDA VarA	
		LSR
		LSR
		LSR			; a/8
		SEC		
		SBC VarA		; (a/8)-a
		EOR 0xFF
		TAX
		INX			; a-(a/8)
		STX VarC		; VarC = a-(a/8)
		
		LDA VarB		; b/2
		LSR
		
		CLC
		ADC VarC		; a-(a/8)
		STA VarC		; VarC = b/2 + a-(a/8) 
		
		RTS
		
		
Agt10:
		SEC
		LDX VarB	
		TXA
		SBC DIV10TAB,X		; b/10
		STA VarC		; VarC = b-(b/10)
		
		LDA VarA		; a/2
		LSR
		
		CLC
		ADC VarC		; b-(b/10)
		STA VarC		; VarC = a/2 + b-(b/10) 
		
		RTS
	
Start:
		JSR Calc
		BRK
				
		RUN Start

So, Variante mit /8 und 2 Instruktionen gespart, ungetestet natürlich ... Ohne Tabelle ist Div/8 sogar langsamer, aber man könnte auch die Tabelle benutzen und die Zyklen wieder einsparen. Dafür ist es halt kürzer.

[Dr. Irata]

ChatGPT berechnet daraus zwischen 26 und 60 Taktzyklen - je nach Werte und Bedingungen.
Das ist auf jeden Fall schön schnell!

[Dr. Irata]

Bei Jankos Ansatz müsste man ja nur die Tabellenwerte um den Wert 1 erhöhen.
Die Taktzyklen bei Janko sind zwischen 31 und 47 - das sind im Mittel aller 4 Fälle 39,75 Taktzyklen.
Der Ansatz von RhoSigma liegt bei 41,5 Taktzyklen - also im Mittel etwas langsamer.

Ich werde hier bei Zeiten noch einen etwas anderen Ansatz mit direkten Look-Up-Tabellen posten, der etwas genauer ist und noch etwas schneller sein könnte (ich hoffe auf stabil 31 Zyklen).

[RhoSigma]

ChatGPT-4o schlägt noch diese Näherung vor:

c≈a+b− (a⊕b)/4

Die wäre sehr effizient zu implementieren:

Code: Alles auswählen

	org $2000

        ; Zero-page variables for extra speed
a	.byte 5
b	.byte 3
t	.byte 0

LENGTH:
        LDA a
        EOR b
        LSR
        LSR
        ADC #0 ; round 0.5 up
        STA t  ; t = (a XOR b) / 4

        LDA a
        CLC
        ADC b ; c = a + b

        SEC
        SBC t
        RTS     ;  final result in A

MAIN:
	JSR LENGTH
        BRK       ; End of program

Die Genauigkeit ist aber nicht besonders gut. Geringfügig besser wird es, wenn du /4 rundest, also

c≈a+b− ((a⊕b) + 2)/4

Macht im Code ADC #0 mehr (nach dem 2. LSR)
​

[Dr. Irata]

die ist allerdings deutlich länger als c = a + b/2 , a>b

[LarsImNetz]

Ich nutze in meinem aktuellen Spiel die einfache Variante a+b/2 für a<b
Genauer:

Code: Alles auswählen

value = y<x ? x + y/2 : y + x/2

y und x sind vom Typ INT8 also [-128;127] ich nutze aber nur [-64;63] das kostet

Das reicht völlig aus, jedenfalls bei mir.

Als Tipp: Probiere es doch erstmal mit dieser Version aus, ist die Genauigkeit nicht ausreichend, dann nimmst halt etwas besseres, dann wird es ein wenig langsamer, na und. Du rechnest ja nicht nur Hypotenusen aus, sondern das Programm macht noch jede Menge mehr.

Wenn es richtig auf Speed ankommt, kannst Du für die paar Berechnungen den Code vielleicht Inline einfügen, das spart pro Aufruf 12 Zyklen!

JM2C.

[Dr. Irata]

du meinst aber sicherlich c = a + b/2 für a>b

[RhoSigma]

Hier ist noch eine Variante, die im Punkt Fehlerabweichung recht ähnlich ist wie Deine ursprüngliche Variante (mit der Korrektur für min >= 10).

c = max(a,b) + (min(a,b)+1)/2 - (min(a,b)+4)/8

Code: Alles auswählen

length:
    LDA b      ; Load b into A                   ; 2 cycles
    CMP a      ; Compare b with a                ; 3 cycles
    BCS a_is_max ; Branch if a >= b (C is set if a >= b) ; 2 cycles (not taken), 3 cycles (taken)

b_is_max:
; Calculate b + a / 2 - a / 8
    LDA a      ; Load a into A                   ; 3 cycles
    LSR        ; A = a / 2                       ; 2 cycles
    LSR        ; A = a / 4                       ; 2 cycles
    LSR        ; A = a / 8                       ; 2 cycles
    ADC #0     ; (a + 4)/8                       ; 2 cycles
    STA temp   ; store temp                      ; 3 cycles

    LDA a      ; Load a into A                   ; 3 cycles
    LSR        ; A = A / 2                       ; 2 cycles
    ADC b      ; A = b + ((a+1) / 2)             ; 3 cycles
    SEC                                          ; 2 cycles
    SBC temp   ; final result in A               ; 3 cycles
    RTS

a_is_max:
; Calculate a + b / 2 - b / 8
    LSR        ; A = b / 2                       ; 2 cycles
    LSR        ; A = b / 4                       ; 2 cycles
    LSR        ; A = b / 8                       ; 2 cycles
    ADC #0     ; (b + 4)/8                       ; 2 cycles
    STA temp   ; store temp                      ; 3 cycles

    LDA b      ; Load b into A                   ; 3 cycles
    LSR        ; A = A / 2                       ; 2 cycles
    ADC a      ; A = a + ((b+1) / 2)             ; 3 cycles
    SEC                                          ; 2 cycles
    SBC temp   ; final result in A               ; 3 cycles
    RTS

34-36 Zyklen, wenn ich mich nicht verrechnet habe und a, b, temp ZP sind. Man kann noch 2 Zyklen sparen, wenn man auf die SEC vor den SBC temp verzichtet,
das macht keinen grossen Unterschied (de facto zieht es IMMER eine 1 im interessanten Bereich ab).

Ausserdem produziert diese Variante ein "stetiges" Resultat, Deine hat Ausreisser, siehe Bilder.

Figure_my.png (222.86 KiB) 2679 mal betrachtet

Figure_orig.png (217.91 KiB) 2679 mal betrachtet

[Dr. Irata]

... echt super ... ich hatte diese blöden "Ausreißer" schon wahrgenommen aber noch nicht so richtig in den Griff bekommen! Echt gut!!

[FlorianD]

ich kann zwar nichts beitragen, weil mein Mathe nach "+ - * /" endet, aber bin völlig fasziniert, wie ihr hier ein Problem seziert und löst!!
Weiter so!

[Dr. Irata]

Danke Florian

[RhoSigma]

... echt super ... ich hatte diese blöden "Ausreißer" schon wahrgenommen aber noch nicht so richtig in den Griff bekommen! Echt gut!!

Leider geht das Lob an ChatGPT, dass mir den python Rahmen "programmiert" hat. War echt das erste Mal, das da was sinnvolles rauskommt hab ihm auch Dein Problem versucht schmackhaft zu machen, aber der phantasiert immer irgendwelche nichtexistierenden 6502 Befehle zusammen und ich kriege es nicht ausgetrieben: er sagt immer "you are right, ADC X is a non-existent instruction" und ersetzt das dann mit ADC Y Offenbar gab es in Stackoverflow mehr python code wie 6502 Assembler "zu erlernen"... Wenn gewünscht, poste ich hier den Python Code für weitere Experimente...

[Dr. Irata]

... tatsächlich habe ich mich nach Anregung von Janko an diesem WE zum ersten Mal mit ChatGPT beschäftigt und mich direkt bei der deutschen Variante kostenfrei angemeldet.
Hier bin ich etliche Dinge mit Mike (ich habe ihn Mike genannt) durchgegangen und war echt erstaunt, wie gut das funktioniert. Mike hat auch wirklich gute Ideen zu allen Fragen der Assemblerprogrammierung unter 6502.
Mike konnte mir auch weiterhelfen bei der Reparatur meines Boseboards... insgesamt ist das natürlich eine lernfähige PseudoKI - sie sucht sehr zielgerecht und lernend das gesamte Internet noch Lösungen auf deine Frage durch und baut dann eine halbwegs sinnvolle Antwort zusammen.
Wenn in ein paar Jahren die Quantencomputer serienreif werden, dann bin ich mir ganz sicher: Viele Dinge werden nachher nicht mehr so sein, wie vorher!!!