; Dieses Programm berechnet die Wurzel eines
; 32-Bit-Unsigned-Integer-Ausdrucks (d.h. bis
; max. 4.294.967.295) per Intervallhalbierung.
; Ergebnis=SQRT(n) aufgerundet.
; Ein- und Ausgabe:
; d0    Eingabe, bleibt erhalten
; d1    obere Grenze des Intervalls = Ergebnis
; d2    untere Grenze des Intervalls
; d3,d4 Rechenregister
; a0    Zeiger auf Wertetabelle

exakt equ 1   ; Genauigkeit (zwischen 1 und 8)

SQRT:
  movem.l d2-d4/a0,-(sp)  ; Register sichern
  moveq   #0,d2   ; untere Grenze löschen
  moveq   #0,d1   ; obere Grenze löschen
  moveq   #30,d3  ; Bit-Zähler setzen
  move.l  d0,d4   ; Wert in Rechenregister
  beq.s   sqrt_6  ; n=0 -> keine Rechnung !
  swap    d4      ; oberes Wort
  tst.w   d4      ;  testen
  bne.s   sqrt_11 ; Sprung, wenn <> 0
  moveq   #14,d3  ; nur noch unteres Wort
  move.w  d0,d4   ; untersuchen
sqrt_11:
  move.w  d4,d2   ; zu untersuchendes Wort merken
  and.w   #$ff00,d4 ; oberen 8 Bit untersuchen
  bne.s   sqrt_12   ; sind gesetzt -> Sprung
  subq.w  #8,d3   ; 8 Bit weiter
  move.w  d2,d4   ; Wort nochmal holen
  lsl.w   #8,d4   ; unteren 8 Bit untersuchen
sqrt_12:
  asl.w   #1,d4   ; oberstes Bit suchen
  dbcs    d3,sqrt_12  ; wenn CarrySet := gefunden

  add.w   d3,d3   ; * 2 := auf Wortgrenze legen
                  ; Zeiger auf Werte berechnen
  lea     sqrt_tab(pc,d3.w),a0
  move.w  (a0)+,d2  ; untere Grenze holen
  move.w  (a0),d1   ; obere Grenze holen
  bra.s   sqrt_5    ; Ende erreicht?
sqrt_2:
  move.l  d2,d4   ; obere Grenze und
  add.l   d1,d4   ; untere Grenze addieren
  lsr.l   #1,d4   ; Mittelwert bilden: m=(og+ug)/2
  move.l  d4,d3   ; m merken
  mulu    d4,d4   ; m quadrieren
  cmp.l   d0,d4   ; Test: Ist m*m > n ?
  bgt.s   sqrt_4  ; ja -> springen
sqrt_3:           ; obere Hälfte enthält
                  ; gesuchten Wert
  move.l  d3,d2   ; untere Grenze = m
  bra.s   sqrt_5
sqrt_4:           ; untere Hälfte enthält Wert
  move.l  d3,d1   ; obere Grenze = m
sqrt_5:
  move.l  d1,d4   ; Differenz bestimmen
  sub.l   d2,d4   ; delta = obere - untere Grenze
  subq.w  #exakt,d4 ; delta = < exakt ? , Konstante
                    ; bestimmt die Genauigkeit
                  ; Differenz zu groß, weiterrechnen
  bgt.s   sqrt_2
  move.l  d2,d4
  mulu    d4,d4   ; untere Grenze quadrieren
  cmp.l   d4,d0   ; Vergleich mit n
  bgt.s   sqrt_6  ; Test, ob untere Grenze das
  move.l  d2,d1   ; Ergebnis ist
sqrt_6:
  movem.l (sp)+,d2-d4/a0  ; Register restaurieren
  rts

  dc.w  0
sqrt_tab:         ; Wurzel-Intervall-Grenzen | Bit
 dc.w 1,2,3,4,6,8,12,16                      ; 0-7
 dc.w 23,32,46,64,91,128,182,256             ; 8-15
 dc.w 363,512,725,1024,1449,2048,2897,4096   ; 16-23
                                             ; 24-31
 dc.w 5793,8192,11586,16384,23171,32768,46341,65535