p-q-Formel für Assembler

[Dr. Irata]

Als nächstes wichtiges Mathe-Modul möchte ich hier an dieser Stelle noch eine schöne Routine für die p-q-Formel erarbeiten und einstellen.

Was steckt dahinter?
Die Lösung einer quadratischen Gleichung in der Normalform x²+px+q=0 mit quadratischer Ergänzung nach x.
Das Ergebnis ist die allgemeine pq Formel:

Bildschirmfoto 2024-01-18 um 12.49.57.png (13.36 KiB) 3946 mal betrachtet

Und dafür bauen wir jetzt eine Routine, alle dafür benötigten Mittel haben wir bereits im Vorfeld erarbeitet.

[Dr. Irata]

... hat jetzt etwas länger gedauert, aber diese Assemblerroutine steht ja noch aus.

Die pq-Formel hat grundsätzlich 0,1 oder 2 Lösungen.
Das muss man relativ am Anfang der Routine checken, sonst gibt es Chaos.
Der erste Teil der Routine wird also das Ergebnis unter der Wurzel betrachten, nachdem es ausgerechnet worden ist:
Wenn das Ergebnis negativ ist, gibt es insgesamt keine Lösung für x, wenn das Ergebnis 0 ist, dann gibt es genau eine Lösung für x und wenn das Ergebnis positiv ist, dann gibt es zwei Lösungen (x1 und x2).
Also 1. Schritt das Ergebnis unter der Wurzel ausrechen und falls das Ergebnis größer 0 ist, die Wurzel davon ziehen.

[Dr. Irata]

Insgesamt nutzen wir für dieses Beispiel keine Fließkommazahlen, daher wird das Ergebnis immer nur eine Näherung sein. Allerdings brauchen wir von Anfang an auch negative Zahlen, daher beschäftigen wir uns zunächst damit:
Da ich mit 8 Bit Zahlen rechnen möchte, nehme ich für die Kennung des Vorzeichens eine eine weitere Variable...

[cas]

"Dumme" (oder nicht so "dumme") Frage: wofür wird die Formel benötigt? Was kann ich damit machen? Warum sollte ich mich damit beschäftigen?

[Dr. Irata]

Die pq-Formel ist wohl eine der essentiellen Formeln in der Mathematik. Damit löst man quadratische Gleichungen!
Wenn man sich in der Mathematik bewegt, wird man wohl um die pq-Formel kaum herumkommen. Ob man die jetzt zwingend beim programmieren braucht, weiß ich auch nicht. Also zumindest beim Raycaster braucht es doch einiges an Mathematik (Pythagoras z.B.) - die pq-Formel ergänzt quasi diese kleine mathematische "Reise".....

[Dr. Irata]

... puh ... war doch aufwendiger, als gedacht! Aber noch in diesem Jahr fertig gestellt!!
Hier nun ist der Code für die pq-Formel (wegen fehlender Fließkommazahlen natürlich nur eine Näherung):

Bitte gerne testen... hat mich gerade mein Hirn zum drehen gebracht

LG und guten Rutsch
Peter

Code: Alles auswählen

var1	.by 08		;p
var2	.by 01		;q
var3	.by 00		;x1
var4	.by 00		;x2

minus1	.by 00		;Minus für var1
minus2	.by 00		;Minus für var2
minus3	.by 00		;Minus für var3
minus4	.by 00		;Minus für var4

temp1	.by 00,00	;Zwischenspeicher1
temp2	.by 00,00	;Zwischenspeicher2

anzahl	.by 02		;Anzahl der Lösungen


		
		.proc main
		
		;....................... pq-Formel .......
		; x*x+px+q=0
		; x1,2 = -p/2 ±√(p*p/4 - q)
		
		
		; ----------------- p*p ausrechnen ----------------------------------------
		
			  
		lda var1       	   ; Lade p in den Akku			
           	sta temp1          ; Speichere p im niedrigeren Byte von temp1
           	lda #0             ; Setze Akku auf 0 (wird als Akkumulator verwendet)
           	ldx #8             ; Initialisiere Schleifen-Zähler (8 Bits zu verarbeiten)
 
           	lsr temp1          ; Shift temp1 (niedriges Byte) um 1 nach rechts
lp1:
           	bcc lp2            	; Überspringe Addition, wenn kein Carry-Flag gesetzt
           	clc               	; Clear Carry vor Addition
           	adc var1          	; Addiere var1 (aktuelles p) zum Akkumulator
lp2:
            	ror             	; Rotate Right Akku (höheres Bit nach Carry)
            	ror temp1         	; Rotate Right temp1 (Carry ins höchstwertige Bit)
            	dex			; Dekrementiere Schleifen-Zähler
            	bne lp1			; Wiederhole, solange Zähler nicht null
 
           	sta temp1+1       	; Speichere höheres Byte des Ergebnisses in temp1+1
					; jetzt durch 4 teilen
		clc
		ror temp1+1
		ror temp1
		clc
		ror temp1+1
		ror temp1
			
		;----------------- Subtraktion q -------------------------------------
		#if .byte minus2=#0		; für den Fall q ist positiv
		sec 
           	lda temp1 
            	sbc var2 
            	bcs weiter 
            	#if .byte temp1+1=#0
            	lda #0
            	sta anzahl	;keine Lösung für x
		jmp fertig
		#end
            	dec temp1+1 
            	jmp weiter
           	#end
           	clc					; für den Fall q ist negativ
            	lda temp1
            	adc var2
            	sta temp1
            	lda temp1+1
            	adc #0
           	sta temp1+1
            	lda temp1
            
weiter     	sta temp1 
		
		;-----------------------------------------
            	; Quadratwurzelberechnung (16-Bit)
           	; Eingabe: temp1/temp1+1 (16-Bit-Wert)
           	; Ausgabe: Quadratwurzel in temp2/temp2+1
           	;-----------------------------------------
 

           	lda #0             ; Akku auf 0 setzen
           	sta temp2           ; temp2 (niedriges Byte des Ergebnisses) initialisieren
           	sta temp2+1         ; temp2+1 (hohes Byte des Ergebnisses) initialisieren
      		ldx #8             ; Anzahl der Bits (8 Iterationen)
 
lp3
            	sec                ; Carry-Flag setzen (für Subtraktion)
           	lda temp1+1         ; Höheres Byte von temp1 laden
            	sbc #64            ; Subtrahiere 64 (Schrittgröße)
            	tay                ; Zwischenergebnis in Y speichern
           	lda temp2+1         ; Lade aktuelles Ergebnis (hohes Byte von temp2)
            	sbc temp2           ; Subtrahiere Zwischenergebnis
           	bcc lp4            ; Wenn Carry = 0, Rücknahme der Subtraktion
 
            	sty temp1+1         ; Speichere angepasstes Höheres Byte von temp1
           	sta temp2+1         ; Aktualisiere hohes Byte des Ergebnisses
 
lp4
           	rol temp2           ; Rotate temp2 (niedriges Byte) nach links
           	asl temp1           ; Shift temp1 (niedriges Byte) nach links
           	rol temp1+1         ; Shift temp1+1 (höheres Byte) nach links
           	rol temp2+1         ; Shift temp2+1 (hohes Byte des Ergebnisses) nach links
 
           	asl temp1           ; Wiederholen für Präzision
           	rol temp1+1
           	rol temp2+1
 
           	dex                 ; Dekrementiere Schleifen-Zähler
           	bne lp3             ; Wiederhole, bis alle Bits verarbeitet sind
		;------------------- Wurzelberechnung fertig -----------------------------
		
		lsr var1	    ; p/2 bilden
			
		#if .byte temp2=#0
		lda #1 
		sta anzahl
		sta minus3
		mva var1 var3	    ; p/2 als Ergebnis in x1 speichern
		#if .byte minus1=#1
		mva #0 minus3	    ; in diesem Fall wird das Ergebnis positiv
		#end
		jmp fertig
		#end
		
				    ; ab hier jetzt x1 und x2 berechnen
				    ; x1 = -p/2 + Wurzel  /  x2 = -p/2 - Wurzel

		clc
		lda var1	    ; p/2
		adc temp2	    ; + Wurzelwert 
		sta var3	    ; x1
		
		lda var1	    ; p/2
		sec
		sbc temp2	    ; - Wurzelwert
		bpl positiv	    ; den Wert als Betrag bilden
negativ
		eor #$ff
		clc
		adc #1
positiv	
		sta var4	    ; x2
		
				    ; jetzt noch die Vorzeichen korrigieren
		#if .byte temp2>var1
		mva #1 minus3
		mva #0 minus4
		#if .byte minus1=#1 ; wenn p = minus
		mva #0 minus3
		mva #1 minus4
		#end
		jmp fertig
		#end
		
		lda #0
		sta minus3
		sta minus4
		#if .byte minus1=#0
		lda #1
		sta minus3
		sta minus4
		#end
fertig		
			

[LarsImNetz]

Erstmal danke für den Code, aber ich hätte da noch ein paar Anmerkungen:
Sieh es bitte nicht als Kritik gegen Dich Peter, ich finde super gut, was Du hier etwas sendest.

Also:

Code: Alles auswählen

           	sta temp1+1       	; Speichere höheres Byte des Ergebnisses in temp1+1
					; jetzt durch 4 teilen
		clc
		ror temp1+1
		ror temp1
		clc
		ror temp1+1
		ror temp1

ist zu ersetzen durch:

Code: Alles auswählen

  lsr
  ror temp1
  lsr
  ror temp1
  sta temp1+1

weil kürzer und schneller

2. IMHO sieht der Code scheußlich aus. Ich frage mich manchmal, wozu man Variablen überhaupt benennen kann. Warum benennst Du 'p' var1 und 'q' var2?

3. Zu was wird '#if .byte minus2=#0'. Das ist kein Assembler, auch weil der Mac65/Atasm das nicht übersetzt.

4. Welche Variablen darf/muss ich wie setzen?
- p ist scheinbar nur 8bit breit
- q ist scheinbar auch nur 8bit breit.

5. Das Ergebnis landet wohl in x1,2. Wann genau gilt: Es gibt nichts, oder nur x1 oder x1 und x2. Das erschließt sich mir nicht so ganz.
6. Was mit '.proc Name' beginnt, sollte IMHO mit '.end' enden. Oder war das '.endproc' ...

Das scheint mir noch nicht ganz fertig zu sein.

JM2C

[LarsImNetz]

Einen hab ich noch:
Dein Code enthält die Wurzelberechnung nach Lee Davison (http://6502org.wikidot.com/software-math-sqrt).
Die Routinen von (https://www.codebase64.org/doku.php?id=base:fast_sqrt) sind etwas flotter. (so ca. 40%)

IMHO hatten wir das schon mal...

LG
Lars

[Dr. Irata]

Hallo Lars,

ich wollte den Code unbedingt gestern noch reinstellen und vielleicht hätte es noch mehr Codetest und Bereinigung gebraucht.
Um 01:15 konnte ich mich aber einfach nicht mehr so gut konzentrieren.
Ich werde im neuen Jahr den Code noch überarbeiten, die #if / then ersetzen und die Variablen besser lesbar anpassen und noch mehr Erklärungen reinbringen.
Ich schaue mir auch nochmal die andere Wurzelberechnung an...
Danke für die Kritik. Der Code ist auch komplett frei und jeder darf ihn optimieren (wie du es schon beim teilen gemacht hast).

Liebe Grüße und jetzt aber erstmal mit gutem Rutsch ins neue Jahr kommen...

Peter

[Dr. Irata]

Hallo Ihr Lieben und ein frohes neues Jahr wünsche ich allen!
Ich hoffe, alle sind gut ins neue Jahr hineingekommen.

Ich habe jetzt mal die Routine überarbeitet! Hoffe, daß es jetzt besser gefällt.
Die langsamere Wurzeliteration habe ich noch so belassen, das will ich erst noch in meinen Code umsetzen und kann es dann hier einbauen.

Code: Alles auswählen

p	.by 00		;p
q	.by 00		;q
x1	.by 00		;x1
x2	.by 00		;x2

minusp	.by 00		;Minus für p
minusq	.by 00		;Minus für q
minusx1	.by 00		;Minus für x1
minusx2	.by 00		;Minus für x2

temp1	.by 00,00	;Zwischenspeicher1
temp2	.by 00,00	;Zwischenspeicher2

anzahl	.by 00		;Anzahl der Lösungen


		
		;....................... pq-Formel .......
		; x*x+px+q=0
		; x1,2 = -p/2 ±√(p*p/4 - q)
		; p und q werden in diesem Beispiel bewußt auf 8 Bit
		; beschränkt, da man ansonsten bei 16 Bit für p und q
		; Zwischenspeicher und Wurzelfunktion auf 32 Bit hätte
		; ausrechnen müssen - das würde hier den Rahmen sprengen

		; ----------------- p*p ausrechnen ----------------------------------------
		
		lda #8		; 8-Bit-Wert in p laden
		sta p
		lda #0
		sta minusp
		
		lda #10		; 8-Bit-Wert in q laden
		sta q
		lda #1
		sta minusq
		
		;--------------- Start pq-Routine -----------------------------------------
			  
		lda p       	; Lade p in den Akku			
           	sta temp1       ; Speichere p im niedrigeren Byte von temp1
           	lda #0          ; Setze Akku auf 0 (wird als Akkumulator verwendet)
           	ldx #8          ; Initialisiere Schleifen-Zähler (8 Bits zu verarbeiten)
 
           	lsr temp1       ; Shift temp1 (niedriges Byte) um 1 nach rechts
lp1
           	bcc lp2         ; Überspringe Addition, wenn kein Carry-Flag gesetzt
           	clc             ; Clear Carry vor Addition
           	adc p          	; Addiere p (aktuelles p) zum Akkumulator
lp2
            	ror		; Rotate Right Akku (höheres Bit nach Carry)
            	ror temp1	; Rotate Right temp1 (Carry ins höchstwertige Bit)
            	dex		; Dekrementiere Schleifen-Zähler
            	bne lp1		; Wiederhole, solange Zähler nicht null
 
           	       	
				; p*p durch 4 teilen
		lsr
  		ror temp1
  		lsr
  		ror temp1
  		sta temp1+1
			
		;----------------- Subtraktion q -------------------------------------
		; für den Fall q ist positiv
		lda minusq
		cmp #0
		bne sp2
		sec 
           	lda temp1 
            	sbc q 
            	bcs weiter 
            	
            	lda temp1+1
            	cmp #0
            	bne sp1
            	lda #0
            	sta anzahl	;keine Lösung für x
		; hier Textausgabe: keine Lösung für x
		jmp fertig
		
sp1		
            	dec temp1+1 
            	jmp weiter
           	
sp2           	
           	; für den Fall q ist negativ
		clc
            	lda temp1
            	adc q
            	sta temp1
            	lda temp1+1
            	adc #0
           	sta temp1+1
            	lda temp1
            
weiter     	sta temp1 
		
		;-----------------------------------------
            	; Quadratwurzelberechnung (16-Bit)
           	; Eingabe: temp1/temp1+1 (16-Bit-Wert)
           	; Ausgabe: Quadratwurzel in temp2/temp2+1
           	;-----------------------------------------
 

           	lda #0		; Akku auf 0 setzen
           	sta temp2       ; temp2 (niedriges Byte des Ergebnisses) initialisieren
           	sta temp2+1     ; temp2+1 (hohes Byte des Ergebnisses) initialisieren
      		ldx #8          ; Anzahl der Bits (8 Iterationen)
 
lp3
            	sec             ; Carry-Flag setzen (für Subtraktion)
           	lda temp1+1     ; Höheres Byte von temp1 laden
            	sbc #64         ; Subtrahiere 64 (Schrittgröße)
            	tay             ; Zwischenergebnis in Y speichern
           	lda temp2+1     ; Lade aktuelles Ergebnis (hohes Byte von temp2)
            	sbc temp2       ; Subtrahiere Zwischenergebnis
           	bcc lp4         ; Wenn Carry = 0, Rücknahme der Subtraktion
 
            	sty temp1+1     ; Speichere angepasstes Höheres Byte von temp1
           	sta temp2+1     ; Aktualisiere hohes Byte des Ergebnisses
 
lp4
           	rol temp2       ; Rotate temp2 (niedriges Byte) nach links
           	asl temp1       ; Shift temp1 (niedriges Byte) nach links
           	rol temp1+1     ; Shift temp1+1 (höheres Byte) nach links
           	rol temp2+1     ; Shift temp2+1 (hohes Byte des Ergebnisses) nach links
 
           	asl temp1       ; Wiederholen für Präzision
           	rol temp1+1
           	rol temp2+1
 
           	dex             ; Dekrementiere Schleifen-Zähler
           	bne lp3         ; Wiederhole, bis alle Bits verarbeitet sind
		;------------------- Wurzelberechnung fertig -----------------------------
		
		lsr p	    	; p/2 bilden
		
		lda temp2
		cmp #0
		bne sp4
		lda #1 
		sta anzahl
		sta minusx1
		lda p 
		sta x1	; p/2 als Ergebnis in x1 speichern
		
		
		lda minusp
		cmp #1
		bne sp3
		lda #0 
		sta minusx1	; in diesem Fall wird das Ergebnis positiv
		
sp3		
		; hier Textausgabe: "Lösung für x="x1
		jmp fertig
		
sp4		
				; ab hier jetzt x1 und x2 berechnen
				; x1 = -p/2 + Wurzel  /  x2 = -p/2 - Wurzel

		clc
		lda p		; p/2
		adc temp2	; + Wurzelwert 
		sta x1		; x1
		
		lda p		; p/2
		sec
		sbc temp2	; - Wurzelwert
		bpl positiv	; den Wert als Betrag bilden
negativ
		eor #$ff
		clc
		adc #1
positiv	
		sta x2	    	; x2
		
				; jetzt noch die Vorzeichen korrigieren
		lda p
		cmp temp2
		bcs sp6
		lda #1 
		sta minusx1
		lda #0 
		sta minusx2
		
		lda minusp
		cmp #1
		bne sp5
		lda #0 
		sta minusx1
		lda #1 
		sta minusx2
		
sp5		
		jmp fertig
		
sp6		
		lda #0
		sta minusx1
		sta minusx2
		lda minusp
		cmp #0
		bne sp7
		lda #1
		sta minusx1
		sta minusx2
sp7
		; hier Textausgabe: "Lösung für x1,2=" x1,x2
			
fertig		

[Dr. Irata]

... anbei final die pq-Routine mit der schnellen Wurzelroutine:

Code: Alles auswählen

p	.by 00			;p
q	.by 00			;q
x1	.by 00			;x1
x2	.by 00			;x2

minusp		.by 00		;Minus für p
minusq		.by 00		;Minus für q
minusx1		.by 00		;Minus für x1
minusx2		.by 00		;Minus für x2
temp		.by 00,00
temp1		.by 00,00	;Zwischenspeicher1
temp2		.by 00		;Zwischenspeicher2

anzahl		.by 00		;Anzahl der Lösungen


			;Iterationstabelle für die Wurzelberechnung:

iter:  	 .byte $01,$02,$04,$08,$10,$20,$40,$80	
		
		
		;....................... pq-Formel .......
		; x*x+px+q=0
		; x1,2 = -p/2 ±√(p*p/4 - q)
		; p und q werden in diesem Beispiel bewußt auf 8 Bit
		; beschränkt, da man ansonsten bei 16 Bit für p und q
		; Zwischenspeicher und Wurzelfunktion auf 32 Bit hätte
		; ausrechnen müssen - das würde hier den Rahmen sprengen

		; ----------------- p*p ausrechnen ----------------------------------------
		
		lda #8		; 8-Bit-Wert in p laden
		sta p
		lda #0
		sta minusp
		
		lda #10		; 8-Bit-Wert in q laden
		sta q
		lda #1
		sta minusq
		
		;--------------- Start pq-Routine -----------------------------------------
			  
		lda p       	; Lade p in den Akku			
           	sta temp1       ; Speichere p im niedrigeren Byte von temp1
           	lda #0          ; Setze Akku auf 0 (wird als Akkumulator verwendet)
           	ldx #8          ; Initialisiere Schleifen-Zähler (8 Bits zu verarbeiten)
 
           	lsr temp1       ; Shift temp1 (niedriges Byte) um 1 nach rechts
lp1
           	bcc lp2         ; Überspringe Addition, wenn kein Carry-Flag gesetzt
           	clc             ; Clear Carry vor Addition
           	adc p          	; Addiere p (aktuelles p) zum Akkumulator
lp2
            	ror		; Rotate Right Akku (höheres Bit nach Carry)
            	ror temp1	; Rotate Right temp1 (Carry ins höchstwertige Bit)
            	dex		; Dekrementiere Schleifen-Zähler
            	bne lp1		; Wiederhole, solange Zähler nicht null
 
           	       	
				; p*p durch 4 teilen
		lsr
  		ror temp1
  		lsr
  		ror temp1
  		sta temp1+1
			
		;----------------- Subtraktion q -------------------------------------
		; für den Fall q ist positiv
		lda minusq
		cmp #0
		bne sp2
		sec 
           	lda temp1 
            	sbc q 
            	bcs weiter 
            	
            	lda temp1+1
            	cmp #0
            	bne sp1
            	lda #0
            	sta anzahl	;keine Lösung für x
		; hier Textausgabe: keine Lösung für x
		jmp fertig
		
sp1		
            	dec temp1+1 
            	jmp weiter
           	
sp2           	
           	; für den Fall q ist negativ
		clc
            	lda temp1
            	adc q
            	sta temp1
            	lda temp1+1
            	adc #0
           	sta temp1+1
            	lda temp1
            
weiter     	sta temp1 
		
		;-----------------------------------------
            	; Quadratwurzelberechnung (16-Bit)
           	; Eingabe: temp1/temp1+1 (16-Bit-Wert)
           	; Ausgabe: Quadratwurzel in temp2/temp2+1
           	;-----------------------------------------
 

		ldy #$00   
        	ldx #$07
        	clc         
loopwu
       	 	tya
       		ora iter-1,X
       		sta temp+1   
        	lda temp1+1
        	bcs wu0  
       		cmp temp+1
        	bcc wu1  ; temp <= a
wu0
        	sbc temp+1
        	sta temp1+1   
        	tya
        	ora iter,x
        	tay        
wu1
        	asl temp1
        	rol temp1+1  
        	dex
        	bne loopwu
        	; letzte Iteration
        	bcs wu2
        	sty temp+1
        	lda temp1
        	cmp #$80
        	lda temp1+1
        	sbc temp+1
        	bcc wu3
wu2
        	iny        
wu3
		sty temp2
		;------------------- Wurzelberechnung fertig -----------------------------
		
		lsr p	    	; p/2 bilden
		
		lda temp2
		cmp #0
		bne sp4
		lda #1 
		sta anzahl
		sta minusx1
		lda p 
		sta x1	; p/2 als Ergebnis in x1 speichern
		
		
		lda minusp
		cmp #1
		bne sp3
		lda #0 
		sta minusx1	; in diesem Fall wird das Ergebnis positiv
		
sp3		
		; hier Textausgabe: "Lösung für x="x1
		jmp fertig
		
sp4		
				; ab hier jetzt x1 und x2 berechnen
				; x1 = -p/2 + Wurzel  /  x2 = -p/2 - Wurzel

		clc
		lda p		; p/2
		adc temp2	; + Wurzelwert 
		sta x1		; x1
		
		lda p		; p/2
		sec
		sbc temp2	; - Wurzelwert
		bpl positiv	; den Wert als Betrag bilden
negativ
		eor #$ff
		clc
		adc #1
positiv	
		sta x2	    	; x2
		
				; jetzt noch die Vorzeichen korrigieren
		lda p
		cmp temp2
		bcs sp6
		lda #1 
		sta minusx1
		lda #0 
		sta minusx2
		
		lda minusp
		cmp #1
		bne sp5
		lda #0 
		sta minusx1
		lda #1 
		sta minusx2
		
sp5		
		jmp fertig
		
sp6		
		lda #0
		sta minusx1
		sta minusx2
		lda minusp
		cmp #0
		bne sp7
		lda #1
		sta minusx1
		sta minusx2
sp7
		; hier Textausgabe: "Lösung für x1,2=" x1,x2
			
fertig