Pythagoras goes Assembler

[Dr. Irata]

Wer kennt ihn nicht, den Satz des Pythagoras:
"In allen rechtwinkligen Dreiecken ist die Summe der Flächen der Katheten- Quadrate gleich der Fläche des Quadrates der Hypotenuse"

A^2 + B^2 = C^2 ----> C = (A^2 + B^2)^0.5 also: C = Wurzel aus A Quadrat plus B Quadrat

Wir haben jetzt alles, um diese wichtige Formel in Assembler umzusetzen.
Im Anschluß brauchen wie noch die p/q-Formel, sin/cos, e/ln

Nun zu Pythagoras (den brauche ich übrigens dringend für den Raycaster):

Das Problem ist in Assembler mit 8 Bit, daß man sehr schnell bei A^2 + B^2 auf eine 16-Bit Zahl kommt und dann auch noch die Wurzel dieser 16-Bit Zahl ziehen muss. Das Ergebnis kann wieder eine 8 Bit Zahl sein, aber zwischendurch muss man sehr wahrscheinlich durch die 16 Bit marschieren. Das ist jetzt aber alles gar kein Problem mehr. Wir bilden erst das Produkt A*A, dann B*B, macht dann eine 16-Bit Addiotion von beiden Produkten und davon dann ein 16-Bit Wurzel ziehen - und schwubs hat man das Ergebnis. Daraus bauen wir eine schöne kompakte Subroutine...

Liebe Grüße
Peter

[LarsImNetz]

Hi Peter,
eigentlich suchst Du die Funktion hypot(a,b)
Berechnen einer Hypotenuse.
In meiner Compiler-Lib habe ich 2 Implementierungen gefunden, die Dir beide vielleicht nicht gefallen.
Aber hier mal ein Beispiel:

; @see: https://arxiv.org/pdf/1904.09481.pdf

Code: Alles auswählen

; hypot(a,b) returns sqrt(a^2+b^2) big problem is overrun
; if a == 0 then
;    h = b
;  else
;    r = b / a
;    h = a * sqrt(1 + r * r)
;  end if
; BUT: This works only with float values

Ich könnte meinen Code hier einfach abladen, aber:
ich habe ein Problem mit den Parametern, wenn mein Compiler einen Funktionsaufruf macht:

Code: Alles auswählen

 a := function_call(a,b,c)

Werden die Parameter a,b,c vom Compiler immer erst einmal auf 16 Bit "aufgeblasen" und dann auf einem eigenen Heap abgelegt:
Dieser @heap_ptr liegt irgendwo im Speicher (Cassettenbuffer, oder am Ende des Sources...)

Code: Alles auswählen

 lda a
 ldy #1
 sta (@heap_ptr),y
 lda #0
 iny
 sta (@heap_ptr),y
 lda b
 ldy #3
 sta (@heap_ptr),y
 lda b+1
 iny
 sta (@heap_ptr),y
 lda c
 ldy #5
 sta (@heap_ptr),y
 lda #0
 iny
 sta (@heap_ptr),y
 jsr function_call_iii
 ...

In allen selbst geschriebenen Funktionen fängt mein Source deshalb meist erstmal an die Parameter auszulesen:

Code: Alles auswählen

function_call_iii
 ldy #1
 lda (@heap_ptr),y
 sta parameter_a
 ldy #3
 lda (@heap_ptr),y
 sta parameter_b
 iny
 lda (@heap_ptr),y
 sta parameter_b+1
 iny
 lda (@heap_ptr),y
 sta parameter_c 
 ...

Das ist straightforward für meinen Compiler, diesen ganzen .proc krams (von Mads/cc65) kann der atasm Assembler nicht, den ich ja vorzugsweise verwende. Des Weiteren kommen Rückgabewerte immer in den X-Y-Registern zurück, immer der Form: (y+x*256)

Desweiteren gibt es 'lokale' Variablen, diese beginnen immer mit '?'. Und ein '.local' ändert intern im Atasm nur einen Zähler, so das Variablen, die mit '?' beginnen solange gleich sind, bis ein erneutes '.local' kommt.

Soviel zu meinen Routinen.

Hier eine Implementierung:

Code: Alles auswählen

 .local
 ?POW2 = @REG+1

; Funktion berechnet die Hypotenuse von 2 Parametern (jeweils 8 Bit gross)
@HYPOT
@HYPOT_II
  LDY #1
  LDA (@HEAP_PTR),Y ; wert a holen
  TAX

  LDY #3
  LDA (@HEAP_PTR),Y ; wert b holen
  TAY

  clc
  lda @pow2_low,x
  adc @pow2_low,y
  sta ?pow2
  lda @pow2_high,x
  adc @pow2_high,y
  sta ?pow2+1
; TODO: wenn carry gesetzt ist, ist der Wert > 65535
;  lda #0
;  adc #0
;  sta ?pow2+2

  jmp @SQRT_INNER

@pow2_low
 .byte <0, <1, <4, <9, <16, <25, <36, <49, <64, <81, <100, <121, <144, <169, <196, <225
 .byte <256, <289, <324, <361, <400, <441, <484, <529, <576, <625, <676, <729, <784, <841, <900, <961
 .byte <1024, <1089, <1156, <1225, <1296, <1369, <1444, <1521, <1600, <1681, <1764, <1849, <1936, <2025, <2116, <2209
 .byte <2304, <2401, <2500, <2601, <2704, <2809, <2916, <3025, <3136, <3249, <3364, <3481, <3600, <3721, <3844, <3969
 .byte <4096, <4225, <4356, <4489, <4624, <4761, <4900, <5041, <5184, <5329, <5476, <5625, <5776, <5929, <6084, <6241
 .byte <6400, <6561, <6724, <6889, <7056, <7225, <7396, <7569, <7744, <7921, <8100, <8281, <8464, <8649, <8836, <9025
 .byte <9216, <9409, <9604, <9801, <10000, <10201, <10404, <10609, <10816, <11025, <11236, <11449, <11664, <11881, <12100, <12321
 .byte <12544, <12769, <12996, <13225, <13456, <13689, <13924, <14161, <14400, <14641, <14884, <15129, <15376, <15625, <15876, <16129
 .byte <16384, <16641, <16900, <17161, <17424, <17689, <17956, <18225, <18496, <18769, <19044, <19321, <19600, <19881, <20164, <20449
 .byte <20736, <21025, <21316, <21609, <21904, <22201, <22500, <22801, <23104, <23409, <23716, <24025, <24336, <24649, <24964, <25281
 .byte <25600, <25921, <26244, <26569, <26896, <27225, <27556, <27889, <28224, <28561, <28900, <29241, <29584, <29929, <30276, <30625
 .byte <30976, <31329, <31684, <32041, <32400, <32761, <33124, <33489, <33856, <34225, <34596, <34969, <35344, <35721, <36100, <36481
 .byte <36864, <37249, <37636, <38025, <38416, <38809, <39204, <39601, <40000, <40401, <40804, <41209, <41616, <42025, <42436, <42849
 .byte <43264, <43681, <44100, <44521, <44944, <45369, <45796, <46225, <46656, <47089, <47524, <47961, <48400, <48841, <49284, <49729
 .byte <50176, <50625, <51076, <51529, <51984, <52441, <52900, <53361, <53824, <54289, <54756, <55225, <55696, <56169, <56644, <57121
 .byte <57600, <58081, <58564, <59049, <59536, <60025, <60516, <61009, <61504, <62001, <62500, <63001, <63504, <64009, <64516, <65025

@pow2_high
 .byte >0, >1, >4, >9, >16, >25, >36, >49, >64, >81, >100, >121, >144, >169, >196, >225
 .byte >256, >289, >324, >361, >400, >441, >484, >529, >576, >625, >676, >729, >784, >841, >900, >961
 .byte >1024, >1089, >1156, >1225, >1296, >1369, >1444, >1521, >1600, >1681, >1764, >1849, >1936, >2025, >2116, >2209
 .byte >2304, >2401, >2500, >2601, >2704, >2809, >2916, >3025, >3136, >3249, >3364, >3481, >3600, >3721, >3844, >3969
 .byte >4096, >4225, >4356, >4489, >4624, >4761, >4900, >5041, >5184, >5329, >5476, >5625, >5776, >5929, >6084, >6241
 .byte >6400, >6561, >6724, >6889, >7056, >7225, >7396, >7569, >7744, >7921, >8100, >8281, >8464, >8649, >8836, >9025
 .byte >9216, >9409, >9604, >9801, >10000, >10201, >10404, >10609, >10816, >11025, >11236, >11449, >11664, >11881, >12100, >12321
 .byte >12544, >12769, >12996, >13225, >13456, >13689, >13924, >14161, >14400, >14641, >14884, >15129, >15376, >15625, >15876, >16129
 .byte >16384, >16641, >16900, >17161, >17424, >17689, >17956, >18225, >18496, >18769, >19044, >19321, >19600, >19881, >20164, >20449
 .byte >20736, >21025, >21316, >21609, >21904, >22201, >22500, >22801, >23104, >23409, >23716, >24025, >24336, >24649, >24964, >25281
 .byte >25600, >25921, >26244, >26569, >26896, >27225, >27556, >27889, >28224, >28561, >28900, >29241, >29584, >29929, >30276, >30625
 .byte >30976, >31329, >31684, >32041, >32400, >32761, >33124, >33489, >33856, >34225, >34596, >34969, >35344, >35721, >36100, >36481
 .byte >36864, >37249, >37636, >38025, >38416, >38809, >39204, >39601, >40000, >40401, >40804, >41209, >41616, >42025, >42436, >42849
 .byte >43264, >43681, >44100, >44521, >44944, >45369, >45796, >46225, >46656, >47089, >47524, >47961, >48400, >48841, >49284, >49729
 .byte >50176, >50625, >51076, >51529, >51984, >52441, >52900, >53361, >53824, >54289, >54756, >55225, >55696, >56169, >56644, >57121
 .byte >57600, >58081, >58564, >59049, >59536, >60025, >60516, >61009, >61504, >62001, >62500, >63001, >63504, >64009, >64516, >65025

Aufmerksame Leser sehen sofort, das SQRT() fehlt, das kommt hier:

Code: Alles auswählen

 .local
; fast integer SQRT

?THI = @REG
?MLO = @REG+1
?MHI = @REG+2

; find more here https://www.codebase64.org/doku.php?id=base:fast_sqrt
@SQRT
@SQRT_I
      LDY #1
      LDA (@HEAP_PTR),Y
      STA ?MLO
      INY
      LDA (@HEAP_PTR),Y
      STA ?MHI

        CLC         ; clear bit 16 of M
@SQRT_INNER
        LDY #$00    ; R = 0
        LDX #$07
?loop
        TYA
        ORA ?stab-1,X
        STA ?THI     ; (R ASL 8) | (D ASL 7)
        LDA ?MHI
        BCS ?skip0  ; M >= 65536? then T <= M is always true
        CMP ?THI
        BCC ?skip1  ; T <= M
?skip0
        SBC ?THI
        STA ?MHI     ; M = M - T
        TYA
        ORA ?stab,x
        TAY         ; R = R OR D
?skip1
        ASL ?MLO
        ROL ?MHI     ; M = M ASL 1
        DEX
        BNE ?loop

        ; last iteration

        BCS ?skip2
        STY ?THI
        LDA ?MLO
        CMP #$80
        LDA ?MHI
        SBC ?THI
        BCC ?skip3
?skip2
        INY         ; R = R OR D (D is 1 here)
?skip3
        RTS

?stab   .BYTE $01,$02,$04,$08,$10,$20,$40,$80

LG
Lars

[LarsImNetz]

BTW: Ein Test meinerseits hat ergeben, das die hier angegebene Routine für die Wurzel ca. 45% schneller ist als die von Lee Davison. Aber probiere es selbst aus. Ich habe es auch nicht geglaubt.

[Dr. Irata]

Hallo Lars,
ja so etwas suchte ich - interessanter Ansatz muss ich sagen, sieht aber recht kompliziert aus.
Allerdings haben wir jetzt alles, was man braucht, um c stabil auszurechnen - eine schnelle Multiplikationsroutine mit 16-Bit Ergebnis, die Addition zweier 16-Bit-Zahlen habe ich vorhin gemacht (das ist recht einfach und kurz) und dann diese schnelle Wurzelfunktion einer 16 Bit Zahl. Damit kann man nun recht einfach c ausrechnen!! Ich stelle das in Kürze hier rein, dann können wir zusammen den Code noch optimieren und haben es für gut!
Liebe Grüße
Peter

p.s. An gute und schnelle Algorithmen bin ich immer sehr interessiert!!

[Dr. Irata]

... 45% schneller ist natürlich mega interessant!!

Die Berechnung der Hypothenuse mit r=a/b und dann Wurzel aus (1+r^2) ist ja nur umgeformter Pythagoras. Der Vorteil dabei ist, daß man mit deutlich niedrigeren Zahlen rechnen kann, da man nun nur noch das Quadrat eines Quotienten bilden muss. a=20 und b=10 würde normal Wurzel aus (400+100) also Wurzel aus (500) = 22 Rest 16 = 22,36
Bildet man vorher r=a/b = 20/10 = 2 und dann Wurzel (1+ r*r) also Wurzel 5 = 2,236 * 10 = 22,36 hat man deutlich kleinere Zahlen, allerdings wird man wenn man kein Fließkomma nutzt einen höheren Fehler produzieren!

[Dr. Irata]

... kannst du die Routine noch in reinen Assembler übersetzten? Ich verstehe die Syntax mit dem ? und dem @ nicht...

[LarsImNetz]

Das ist reiner MAC65 Code, das @ hat keinerlei extra Bedeutung. Es ist wie ein Buchstabe. Auch MADS erlaubt es als ersten Wert in einer Variablen, aber hier nochmal der komplette Code, etwas mehr kommentiert.

Code: Alles auswählen

.local

REG=214        ; irgendwo in der Zeropage
_HEAP_PTR=1024 ; wo etwas Platz ist im Speicher

POW2=REG+1    ; + 1 weil SQRT dort den Wert erwartet zum Ziehen der Wurzel

; Funktion berechnet die Hypotenuse von 2 Parametern (jeweils 8 Bit gross)
HYPOT
HYPOT_II
  LDY #1
  LDA (_HEAP_PTR),Y ; wert a holen
  TAX

  LDY #3
  LDA (_HEAP_PTR),Y ; wert b holen
  TAY

; Berechnung von n*n mittels Tabelle, schneller geht es nicht, kostet aber 512 Bytes Speicher!
  CLC
  LDA pow2_low,x
  ADC pow2_low,y
  STA pow2
  LDA pow2_high,x
  ADC pow2_high,y
  STA pow2+1     ; pow2 ist x*x + y*y
                                        

; TODO: wenn carry gesetzt ist, ist der Wert > 65535
;  lda #0
;  adc #0
;  sta pow2+2

  JMP SQRT_INNER

; Enthaelt für den Index n das Low-Byte von n*n 
pow2_low
 .byte <0, <1, <4, <9, <16, <25, <36, <49, <64, <81, <100, <121, <144, <169, <196, <225
 .byte <256, <289, <324, <361, <400, <441, <484, <529, <576, <625, <676, <729, <784, <841, <900, <961
 .byte <1024, <1089, <1156, <1225, <1296, <1369, <1444, <1521, <1600, <1681, <1764, <1849, <1936, <2025, <2116, <2209
 .byte <2304, <2401, <2500, <2601, <2704, <2809, <2916, <3025, <3136, <3249, <3364, <3481, <3600, <3721, <3844, <3969
 .byte <4096, <4225, <4356, <4489, <4624, <4761, <4900, <5041, <5184, <5329, <5476, <5625, <5776, <5929, <6084, <6241
 .byte <6400, <6561, <6724, <6889, <7056, <7225, <7396, <7569, <7744, <7921, <8100, <8281, <8464, <8649, <8836, <9025
 .byte <9216, <9409, <9604, <9801, <10000, <10201, <10404, <10609, <10816, <11025, <11236, <11449, <11664, <11881, <12100, <12321
 .byte <12544, <12769, <12996, <13225, <13456, <13689, <13924, <14161, <14400, <14641, <14884, <15129, <15376, <15625, <15876, <16129
 .byte <16384, <16641, <16900, <17161, <17424, <17689, <17956, <18225, <18496, <18769, <19044, <19321, <19600, <19881, <20164, <20449
 .byte <20736, <21025, <21316, <21609, <21904, <22201, <22500, <22801, <23104, <23409, <23716, <24025, <24336, <24649, <24964, <25281
 .byte <25600, <25921, <26244, <26569, <26896, <27225, <27556, <27889, <28224, <28561, <28900, <29241, <29584, <29929, <30276, <30625
 .byte <30976, <31329, <31684, <32041, <32400, <32761, <33124, <33489, <33856, <34225, <34596, <34969, <35344, <35721, <36100, <36481
 .byte <36864, <37249, <37636, <38025, <38416, <38809, <39204, <39601, <40000, <40401, <40804, <41209, <41616, <42025, <42436, <42849
 .byte <43264, <43681, <44100, <44521, <44944, <45369, <45796, <46225, <46656, <47089, <47524, <47961, <48400, <48841, <49284, <49729
 .byte <50176, <50625, <51076, <51529, <51984, <52441, <52900, <53361, <53824, <54289, <54756, <55225, <55696, <56169, <56644, <57121
 .byte <57600, <58081, <58564, <59049, <59536, <60025, <60516, <61009, <61504, <62001, <62500, <63001, <63504, <64009, <64516, <65025

; Enthaelt für den Index n das High-Byte von n*n 
pow2_high
 .byte >0, >1, >4, >9, >16, >25, >36, >49, >64, >81, >100, >121, >144, >169, >196, >225
 .byte >256, >289, >324, >361, >400, >441, >484, >529, >576, >625, >676, >729, >784, >841, >900, >961
 .byte >1024, >1089, >1156, >1225, >1296, >1369, >1444, >1521, >1600, >1681, >1764, >1849, >1936, >2025, >2116, >2209
 .byte >2304, >2401, >2500, >2601, >2704, >2809, >2916, >3025, >3136, >3249, >3364, >3481, >3600, >3721, >3844, >3969
 .byte >4096, >4225, >4356, >4489, >4624, >4761, >4900, >5041, >5184, >5329, >5476, >5625, >5776, >5929, >6084, >6241
 .byte >6400, >6561, >6724, >6889, >7056, >7225, >7396, >7569, >7744, >7921, >8100, >8281, >8464, >8649, >8836, >9025
 .byte >9216, >9409, >9604, >9801, >10000, >10201, >10404, >10609, >10816, >11025, >11236, >11449, >11664, >11881, >12100, >12321
 .byte >12544, >12769, >12996, >13225, >13456, >13689, >13924, >14161, >14400, >14641, >14884, >15129, >15376, >15625, >15876, >16129
 .byte >16384, >16641, >16900, >17161, >17424, >17689, >17956, >18225, >18496, >18769, >19044, >19321, >19600, >19881, >20164, >20449
 .byte >20736, >21025, >21316, >21609, >21904, >22201, >22500, >22801, >23104, >23409, >23716, >24025, >24336, >24649, >24964, >25281
 .byte >25600, >25921, >26244, >26569, >26896, >27225, >27556, >27889, >28224, >28561, >28900, >29241, >29584, >29929, >30276, >30625
 .byte >30976, >31329, >31684, >32041, >32400, >32761, >33124, >33489, >33856, >34225, >34596, >34969, >35344, >35721, >36100, >36481
 .byte >36864, >37249, >37636, >38025, >38416, >38809, >39204, >39601, >40000, >40401, >40804, >41209, >41616, >42025, >42436, >42849
 .byte >43264, >43681, >44100, >44521, >44944, >45369, >45796, >46225, >46656, >47089, >47524, >47961, >48400, >48841, >49284, >49729
 .byte >50176, >50625, >51076, >51529, >51984, >52441, >52900, >53361, >53824, >54289, >54756, >55225, >55696, >56169, >56644, >57121
 .byte >57600, >58081, >58564, >59049, >59536, >60025, >60516, >61009, >61504, >62001, >62500, >63001, >63504, >64009, >64516, >65025

 .local

; fast integer SQRT

THI = REG
MLO = REG+1

; find more here https://www.codebase64.org/doku.php?id=base:fast_sqrt
SQRT
SQRT_I
      LDY #1
      LDA (_HEAP_PTR),Y
      STA MLO
      INY
      LDA (_HEAP_PTR),Y
      STA MLO+1

      CLC         ; clear bit 16 of M
SQRT_INNER
      LDY #$00    ; R = 0
      LDX #$07
loop
      TYA
      ORA stab-1,X
      STA THI     ; (R ASL 8) | (D ASL 7)
      LDA MLO+1
      BCS skip0  ; M >= 65536? then T <= M is always true
      CMP THI
      BCC skip1  ; T <= M
skip0
      SBC THI
      STA MLO+1     ; M = M - T
      TYA
      ORA stab,x
      TAY         ; R = R OR D
skip1
      ASL MLO
      ROL MLO+1     ; M = M ASL 1
      DEX
      BNE loop

      ; last iteration

      BCS skip2
      STY THI
      LDA MLO
      CMP #$80
      LDA MLO+1
      SBC THI
      BCC skip3
skip2
      INY         ; R = R OR D (D is 1 here)
skip3
; Das Ergebnis ist im Register (y+x*256)
      RTS

stab
      .BYTE $01,$02,$04,$08,$10,$20,$40,$80

[Dr. Irata]

... und hier nun meine Routine für die Berechnung von c nach Pythagoras --> Wurzel aus (a*a+b*b) :

Code: Alles auswählen

		; Pythagoras c = (a*a + b*b)^0,5
		; a=var1
		; b=var2
		; c=var3
		
		.proc pythagoras
		; a^2 ausrechnen
		;............... var1 * var1 = temp1/temp1+1 .......
		;.... Multiplikation 8 Bit * 8 Bit auf 16 Bit ....
		lda var1
		sta temp1
		lda #0
		ldx #8
		lsr temp1
lp1
		bcc lp2
		clc
		adc var1
lp2
		ror
		ror temp1
		dex
		bne lp1
		sta temp1+1
		
		; b^2 ausrechnen
		;............... var2 * var2 = temp2/temp2+1 .......
		;.... Multiplikation 8 Bit * 8 Bit auf 16 Bit ....
		lda var2
		sta temp2
		lda #0
		ldx #8
		lsr temp2
lp3
		bcc lp4
		clc
		adc var2
lp4
		ror
		ror temp2
		dex
		bne lp3
		sta temp2+1
		
		; 16 Bit Addition
		; temp1/temp1+1  +  temp2/temp2+1  =  temp3/temp3+1
		clc
		lda temp1
		adc temp2
		sta temp3
		lda temp1+1
		adc temp2+1
		adc #0
		sta temp3+1
		
		; Quadratwurzel 16 Bit
		; nach Lee Davison
		; Quadratwurzel temp3/temp3+1 ---> Ergebnis in var3
		ldx #8
lp5
		sec
		lda temp3+1
		sbc #64
		tay
		lda temp4
		sbc var3
		bcc lp6
		sty temp3+1
		sta temp4
lp6
		rol var3
		asl temp3
		rol temp3+1
		rol temp4
		asl temp3
		rol temp3+1
		rol temp4
		dex
		bne lp5
		rts
		.endp
		;das Ergebnis steht in var3

... viel Spass beim Testen und Optimieren!!

Liebe Grüße
Peter

[Dr. Irata]

Die Routine ist ja als Subroutine "gekapselt" - das kann man zum testen natürlich auch anders machen und das .proc / .endp weglassen.
Außerdem müssen natürlich alle Variablen vorher initialisiert werden im Assembler z.B.
var1 .byte 00