' Newton-Fraktal

TYPE Komplex
  r AS SINGLE
  i AS SINGLE
END TYPE

DECLARE SUB Aufaddieren (s1s AS Komplex, s2 AS Komplex)
DECLARE SUB Aufmultiplizieren (f1prod AS Komplex, f2 AS Komplex)
DECLARE SUB DividiereNegiere (z AS Komplex, n AS Komplex, q AS Komplex)
DECLARE SUB RGBPalette (f%, r%, g%, b%)
DECLARE FUNCTION Betrag2! (z1 AS Komplex, z2 AS Komplex)

DIM t AS Komplex, z AS Komplex, fx AS Komplex, fdiffx AS Komplex

PRINT "Eingabe Parameter: Nullstellen des Polynoms"
INPUT "Anzahl Nullstellen"; a%
DIM Nullst(a% - 1) AS Komplex, poly(a%) AS Komplex, poly1diff(a% - 1) AS Komplex
' Polynom initialisieren
poly(0).r = 1!
poly(0).i = 0!
FOR i% = 1 TO a%
  poly(i%).r = 0!
  poly(i%).i = 0!
NEXT i%

FOR i% = 1 TO a%
  PRINT i%;
  INPUT ". Nullstelle (Format Real,Imagin„r)"; Nullst(i% - 1).r, Nullst(i% - 1).i
  t.r = -Nullst(i% - 1).r
  t.i = -Nullst(i% - 1).i
  ' Polynom-Koeffizienten bilden (Faltungsprozess)
  FOR j% = i% TO 1 STEP -1
    Aufaddieren poly(j%), poly(j% - 1)
    Aufmultiplizieren poly(j% - 1), t
  NEXT j%
NEXT i%

INPUT "Grenze, wo die Nullstelle als erreicht gilt"; tol!
INPUT "Ausschrittbetrag (Grenze, falls das Newton-Verfahren divergiert)"; diverg!
INPUT "Maximale Iterationstiefe"; itMax%

' 1.Ableitung bilden
t.i = 0!
FOR i% = 1 TO a%
  poly1diff(i% - 1) = poly(i%)
  t.r = CSNG(i%)
  Aufmultiplizieren poly1diff(i% - 1), t
NEXT i%

PRINT "Eingabe Ausschnitt aus der Gauss'schen Ebene"
INPUT "Real min"; xmin!
INPUT "Imagin„r min"; ymin!
INPUT "Imagin„r max"; ymax!
' xmax kann aufgrund der Bildschirmproportionen berechnet werden
' denn wir wollen kein verzerrtes Fraktal
PRINT "Real max ="; xmin! + 4! * (ymax! - ymin!) / 3!

' Ab hier Berechnung
schr! = (ymax! - ymin!) / 480!
tol2! = tol! * tol!
diverg2! = diverg! * diverg!
SCREEN 12

RGBPalette 0, 63, 63, 63
FOR i% = 1 TO 3
  RGBPalette i%, 0, 63 - 5 * i%, 0
  RGBPalette i% + 3, 63 - 5 * i%, 0, 0
  RGBPalette i% + 6, 0, 0, 63 - 5 * i%
  RGBPalette i% + 9, 63 - 5 * i%, 63 - 5 * i%, 0
  RGBPalette i% + 12, 63 - 5 * i%, 0, 63 - 5 * i%
NEXT i%

FOR y% = 0 TO 479
  FOR x% = 0 TO 639
    ' Startwert vom Newton-Verfahren, wo wir beginnen wollen
    z.r = xmin! + CSNG(x%) * schr!
    z.i = ymax! - CSNG(y%) * schr!
    aIt% = 0
    ziel% = -1
    WHILE ziel% = -1 AND aIt% < itMax% AND z.r * z.r + z.i * z.i < diverg2!
      ' Ziel erreicht prfen
      FOR i% = 0 TO a% - 1
        IF Betrag2!(z, Nullst(i%)) < tol2! THEN
          ziel% = i%
        END IF
      NEXT i%
      IF ziel% = -1 THEN
        fx.r = 0!
        fx.i = 0!
        ' Wert f(z) berechnen (Polynom)
        FOR i% = a% TO 0 STEP -1
          Aufmultiplizieren fx, z
          Aufaddieren fx, poly(i%)
        NEXT i%
        ' Wert f'(z) berechnen (Polynom)
        fdiffx.r = 0!
        fdiffx.i = 0!
        FOR i% = a% - 1 TO 0 STEP -1
          Aufmultiplizieren fdiffx, z
          Aufaddieren fdiffx, poly1diff(i%)
        NEXT i%
        ' Wert z(n+)=z(n)-f(z(n))/f'(z(n)) berechnen
        DividiereNegiere fx, fdiffx, t
        Aufaddieren z, t
        aIt% = aIt% + 1
      END IF
    WEND
    IF ziel% = -1 THEN
      PSET (x%, y%), 0
    ELSE
      PSET (x%, y%), 3 * (ziel% MOD 5) + aIt% MOD 3 + 1
    END IF
  NEXT x%
NEXT y%

d$ = INPUT$(1)
SCREEN 0

SUB Aufaddieren (s1s AS Komplex, s2 AS Komplex)
  ' Entspricht s1s += s2; in C
  s1s.r = s1s.r + s2.r
  s1s.i = s1s.i + s2.i
END SUB

SUB Aufmultiplizieren (f1prod AS Komplex, f2 AS Komplex)
  ' Entspricht f1prod *= f2; in C
  DIM t AS Komplex
  t.r = f1prod.r * f2.r - f1prod.i * f2.i
  t.i = f1prod.r * f2.i + f1prod.i * f2.r
  f1prod = t
END SUB

FUNCTION Betrag2! (z1 AS Komplex, z2 AS Komplex)
  dr! = z2.r - z1.r
  di! = z2.i - z1.i
  Betrag2! = dr! * dr! + di! * di!
END FUNCTION

SUB DividiereNegiere (z AS Komplex, n AS Komplex, q AS Komplex)
  ' Diese Prozedur berechnet q = - z / n
  t! = n.r * n.r + n.i * n.i
  q.r = -(z.r * n.r + z.i * n.i) / t!
  q.i = (z.r * n.i - z.i * n.r) / t!
END SUB

SUB RGBPalette (f%, r%, g%, b%)
  PALETTE f%, CLNG(r% + 256 * g%) + 65536 * CLNG(b%)
END SUB