aboutsummaryrefslogtreecommitdiffstats
path: root/tests/t_parsys/float.i68
diff options
context:
space:
mode:
Diffstat (limited to 'tests/t_parsys/float.i68')
-rw-r--r--tests/t_parsys/float.i681226
1 files changed, 1226 insertions, 0 deletions
diff --git a/tests/t_parsys/float.i68 b/tests/t_parsys/float.i68
new file mode 100644
index 0000000..fbfc1e9
--- /dev/null
+++ b/tests/t_parsys/float.i68
@@ -0,0 +1,1226 @@
+; FLOAT.I68
+;-----------------------------------------------------------------------------
+; Fliesskommaroutinen fuer den PC-PAR 68000, Version ohne 68881
+; entnommen mc 11/88, c't...
+
+;-----------------------------------------------------------------------------
+; Definitionen
+
+vorz equ 0
+subflag equ 1
+maxexpo equ 255
+bias equ 127
+extend equ $10
+e equ $402df854 ; exp(1)
+ln2 equ $3f317218 ; ln(2)
+ln10 equ $40135d8e ; ln(10)
+eins equ $3f800000 ; 1.0
+zwei equ $40000000 ; 2.0
+pi2 equ $40c90fdb ; Pi*2
+pi equ $40490fdb ; Pi
+pihalf equ $3fc90fdb ; Pi/2
+
+;-----------------------------------------------------------------------------
+; Librarykopf:
+
+
+S_FloatLib: dc.l S_floatlibend-S_floatlibstart ; Laenge
+S_floatlibstart:
+ dc.l -1 ; Speicher fuer Zeiger
+ dc.b "FLOAT",0 ; Name
+ ds 0
+
+
+
+;-----------------------------------------------------------------------------
+; Sprungtabelle:
+
+ bra.l S_fadd
+ bra.l S_fsub
+ bra.l S_fmul
+ bra.l S_fdiv
+ bra.l S_fmul2
+ bra.l S_fsqrt
+ bra.l S_fabs
+ bra.l S_floatlibnop
+ bra.l S_fcmp
+ bra.l S_fitof
+ bra.l S_fftoi
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_fexp
+ bra.l S_fsinh
+ bra.l S_fcosh
+ bra.l S_ftanh
+ bra.l S_fcoth
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_fln
+ bra.l S_flog
+ bra.l S_fasinh
+ bra.l S_facosh
+ bra.l S_fatanh
+ bra.l S_facoth
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_fsin
+ bra.l S_fcos
+ bra.l S_ftan
+ bra.l S_fcot
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_floatlibnop
+ bra.l S_fasin
+ bra.l S_facos
+ bra.l S_fatan
+ bra.l S_facot
+
+
+;-----------------------------------------------------------------------------
+; Konstanten :
+
+S_Const1 dc.s 1.0
+
+;-----------------------------------------------------------------------------
+; Nullprozedur :
+
+S_floatlibnop: rts
+
+;-----------------------------------------------------------------------------
+; Addition : D0.S = D0.S + D1.S
+
+ ds 0
+S_fadd:
+ addq.l #1,_fadd_cnt.w
+ movem.l d1-d5,-(a7) ; benoetigte Register retten
+ rol.l #1,d0 ; Operanden rotieren und in Form
+ rol.l #1,d1 ; eeee eeee ffff ... fffs bringen
+ move.l d0,d2
+ sub.l d1,d2 ; Differenz beider Zahlen bilden
+ bcc.s fadd_1
+ exg d0,d1 ; ggf. vertauschen, so dass der
+fadd_1: move.b d0,d3 ; kleinere in Register D1 steht
+ and.b #1,d3 ; maskiere das Vorzeichenbit
+ btst #vorz,d2 ; haben beide gleiches Vorzeichen ?
+ beq.s fadd_2 ; bei verschiedenen Vorzeichen
+ bset #subflag,d3 ; Flag fuer Subtraktion setzen
+fadd_2: rol.l #8,d0 ; Form: ffff ... fffs eeee eeee
+ clr.w d4 ; Exponent der ersten Zahl
+ move.b d0,d4 ; wird im Register D4 aufgebaut
+ sne d0 ; falls ungleich Null, dann
+ ror.l #1,d0 ; implizite Eins, sonst implizite
+ clr.b d0 ; Null erzeugen, neu positionieren
+
+ rol.l #8,d1 ; jetzt das gleiche fuer den
+ clr.w d5 ; zweiten Operanden, der Exponent
+ move.b d1,d5 ; kommt ins Register D5
+ sne d1
+ ror.l #1,d1
+ clr.b d1
+
+; In den Registern D0 und D1 stehen jetzt nur noch die Mantissen
+; im Format ffff ... ffff 0000 0000, also linksbuendig, wobei die
+; ehemals implizite Null bzw. Eins nun explizit an erster Stelle steht.
+; In den Registern D4 und D5 stehen die Exponenten der beiden Zahlen.
+; Das Vorzeichen des Ergebnisses sowie die Subtraktionsflags sind im
+; Register D3 zwischengespeichert.
+
+ move.w d4,d2 ; Jetzt Differenz der Exponenten
+ sub.w d5,d2 ; berechnen
+ cmp.w #24,d2 ; groesser als 24 ?
+ bgt.s fadd_rnd ; ja, --> Ergebnis ist groessere Zahl
+ lsr.l d2,d1 ; Mantisse um (D2)-Bits verschieben
+ btst #subflag,d3 ; Subtraktion oder Addition ?
+ bne.s fadd_subtr ; ggf. zur Subtraktion springen
+ add.l d1,d0 ; die beiden Mantissen addieren
+ bcc.s fadd_rnd ; kein Ueberlauf --> zum Runden
+ roxr.l #1,d0 ; Ueberlauf einschieben
+ addq.w #1,d4 ; Exponent korrigieren
+ bra.s fadd_rnd ; und zum Runden
+
+fadd_subtr: sub.l d1,d0 ; die beiden Mantissen subtrahieren
+ beq.s fadd_zero ; bei Null ist das Gesamtergebnis Null
+ bmi.s fadd_rnd ; bei fuehrender Eins zum Runden
+fadd_nrm: tst.w d4 ; Exponent ist schon Null ?
+ beq.s fadd_rnd ; dann ist Ergebnis denormalisiert
+ subq.w #1,d4 ; Exponent erniedrigen
+ lsl.l #1,d0 ; Mantisse normalisieren bis
+ bpl.s fadd_nrm ; fuehrende Eins auftaucht
+
+fadd_rnd: add.l #$80,d0 ; jetzt Runden auf Bit hinter
+ bcc.s fadd_nov ; Mantisse
+ roxr.l #1,d0 ; bei Ueberlauf Mantisse normalisieren
+ addq.w #1,d4 ; und Exponent korrigieren
+fadd_nov: clr.b d0 ; Rest-Mantisse loeschen
+ tst.l d0 ; Ist die Mantisse komplett Null ?
+ beq.s fadd_zero ; ja, dann ist Ergebnis auch Null
+ cmp.w #maxexpo,d4 ; Exponent-Ueberlauf ?
+ blt.s fadd_nue
+ move.w #maxexpo,d4 ; Unendlich Exponent = maxexpo
+ clr.l d0 ; Mantisse = Null
+ bra.s fadd_den
+
+fadd_nue: tst.w d4 ; Exponent Null ( Zahl denormalisiert? )
+ beq.s fadd_den ; ja -->
+ lsl.l #1,d0 ; fuehrendes Bit wird nicht gespeichert
+fadd_den: move.b d4,d0 ; Exponent einsetzen
+ ror.l #8,d0 ; Form: eeee eeee ffff ... fffx
+ roxr.b #1,d3 ; Vorzeichen in Carry schieben
+ roxr.l #1,d0 ; Form: seee eeee efff ... ffff
+
+fadd_zero:
+ movem.l (a7)+,d1-d5 ; Register restaurieren
+ rts ; Ende, Ergebnis steht in D0.L
+
+;-----------------------------------------------------------------------------
+; Subtraktion : D0.S = D0.S - D1.S
+
+ ds 0
+S_fsub:
+ bchg #31,d1 ; Vorzeichen des zweiten Operanden
+ bra S_fadd ; invertieren und zur Addition springen
+
+
+;-----------------------------------------------------------------------------
+; Multiplikation : D0.S = D0.S * D1.S
+
+ ds 0
+S_fmul:
+ addq.l #1,_fmul_cnt.w
+ movem.l d1-d5,-(a7) ; benoetigte Register retten
+ move.l d0,d2 ; Operand 1 kopieren
+ eor.l d1,d2 ; EXOR um Vorzeichen zu bestimmen
+
+ swap d0 ; Registerhaelften Operand 1 vertauschen
+ move.l d0,d3 ; Operand 1 ab jetzt in Register D3
+ and.w #$7f,d3 ; Exponent und Vorzeichen loeschen
+ and.w #$7f80,d0 ; Exponent maskieren
+ beq.s fmul_dn1 ; gleich Null: Zahl ist denormalisiert
+ bset #7,d3 ; implizite Eins einsetzen
+ sub.w #$0080,d0 ; Bias kompensieren
+
+fmul_dn1: swap d1 ; jetzt Operand 2 behandeln
+ move.w d1,d4
+ and.w #$7f,d1
+ and.w #$7f80,d4
+ beq.s fmul_dn2
+ bset #7,d1
+ sub.w #$0080,d4 ; Bias kompensieren
+
+fmul_dn2: add.w d0,d4 ; Exponenten addieren
+ lsr.w #7,d4 ; richtig positionieren
+ sub.w #bias-3,d4 ; Bias-3 subtrahieren
+ cmp.w #-24,d4 ; totaler Unterlauf ?
+ blt.s fmul_zero ; ja, dann ist Ergebnis Null
+
+ move.w d3,d0 ; oberes Mantissenwort von Operand 1
+ mulu d1,d0 ; mal oberem Mantissenwort von Op2
+ swap d0 ; entspricht Verschiebung um 16 Bit
+
+; Das obere Wort von D0 ist nach der Multiplikation auf jeden Fall Null,
+; da die oberen Mantissenworte nur im Bereich 0 ... 255 liegen.
+; Das groete moegliche Ergebnis ist also 255 x 255 = 65025 = 0000FE01.
+; Nach der Vertauschung erhalten wir also eine Zahl der xxxx 0000.
+; Die untere Registerhaelfte von D0 koennen wir kurzzeitig als Zwischen-
+; speicher verwenden.
+
+ move.w d3,d0 ; oberes Wort von Operand 1 merken
+ swap d3 ; jetzt unteres Wort Op1 mal oberes Op2
+ move.w d1,d5
+ mulu d3,d5 ; Ergebnis steht im D5
+ swap d1 ; jetzt unteres Wort Op1 mal unteres Op2
+ mulu d1,d3 ; Ergebnis steht im D3
+ swap d3 ; entspricht Verschiebung um 16 Bit
+ mulu d0,d1 ; jetzt oberes Wort Op1 mal unteres Op2
+
+ move.w d3,d0 ; zum ersten Zwischenergebnis dazu
+ add.l d5,d0 ; jetzt alles aufaddieren
+ add.l d1,d0
+ beq.s fmul_res ; falls Mantisse Null auch Ergebnis Null
+ bmi.s fmul_rnd ; fuehrende Eins? dann zum Runden
+
+; Im Register D0.L befinden sich die oberen 32 Bit des Produktes,
+; im oberen Wort von D3 die restlichen 16 Bit.
+
+ tst.w d4 ; Exponent ist negativ ?
+ bmi.s fmul_unt ; ggf. Unterlauf behandeln
+
+fmul_nor: tst.w d4 ; Exponent = Null ?
+ beq.s fmul_rnd ; falls Null, dann zum Runden
+ roxl.l #1,d3 ; Im oberen Wort von D3 sind die
+ roxl.l #1,d0 ; niedrigsten Bits des Produktes
+ subq.w #1,d4 ; Exponent korrigieren
+ tst.l d0 ; Mantisse testen
+ bpl.s fmul_nor ; bis fuehrende Eins auftaucht
+
+fmul_rnd: add.l #$80,d0 ; Rundung
+ bcc.s fmul_nov
+ roxr.l #1,d0 ; Ueberlauf einschieben
+ addq.w #1,d4 ; Exponent korrigieren
+fmul_nov: cmp.w #maxexpo,d4 ; Exponent-Ueberlauf ?
+ blt.s fmul_nue
+fdiv_err: move.w #maxexpo,d4 ; Ueberlauf: Exponent = Maxexpo
+ clr.l d0 ; Mantisse = Null
+ bra.s fmul_den
+
+fmul_nue: tst.w d4 ; Exponent = Null ?
+ beq.s fmul_den ; falls Null, dann denormalisiert
+ lsl.l #1,d0 ; fuehrende Eins wird nicht abgespeichert
+
+fmul_den: move.b d4,d0 ; Exponent einsetzen
+ ror.l #8,d0 ; Form: eeee eeee ffff ... fffx
+ roxl.l #1,d2 ; Vorzeichen in Carry schieben
+ roxr.l #1,d0 ; und ins Ergebnis einsetzen
+
+fmul_res: movem.l (a7)+,d1-d5 ; Register restaurieren
+ rts
+
+fmul_zero: clr.l d0 ; Null erzeugen
+ bra.s fmul_res ; Ende, Ergebnis steht in D0.L
+
+fmul_unt: cmp.w #-24,d4 ; totaler Unterlauf ?
+ ble.s fmul_zero ; Dann ist das Ergebnis auf jeden Fall Null
+ neg.w d4 ; sonst Shift-Zaehler erzeugen
+ lsr.l d4,d0 ; und Zahl denormalisieren
+ clr.w d4 ; Exponent ist Null als Kennzeichen
+ bra.s fmul_rnd ; fuer eine denormalisierte Zahl
+
+;-----------------------------------------------------------------------------
+; Division : D0.S = D0.S / D1.S
+
+ ds 0
+S_fdiv:
+ addq.l #1,_fdiv_cnt.w
+ movem.l d1-d5,-(a7) ; benoetigte Register retten
+ move.l d0,d2 ; Operand 1 kopieren
+ eor.l d1,d2 ; EXOR um Vorzeichen zu bestimmen
+
+ swap d0 ; Registerhaelften Operand 1 vertauschen
+ move.l d0,d3 ; Operand 1 ab jetzt in Register D3
+ and.w #$7f,d3 ; Exponent und Vorzeichen loeschen
+ and.w #$7f80,d0 ; Exponent maskieren
+ beq.s fdiv_dn1 ; gleich Null: Zahl ist denormalisiert
+ bset #7,d3 ; implizite Eins einsetzen
+ sub.w #$0080,d0 ; Bias kompensieren
+
+fdiv_dn1: swap d1 ; jetzt Operand 2 behandeln
+ move.w d1,d4
+ and.w #$7f,d1
+ and.w #$7f80,d4
+ beq.s fdiv_dn2
+ bset #7,d1
+ sub.w #$0080,d4
+
+fdiv_dn2: sub.w d4,d0 ; Exponenten subtrahieren
+ move.w d0,d4 ; Exponent nach D4 kopieren
+ asr.w #7,d4 ; richtig positionieren
+ add.w #bias,d4 ; Bias addieren
+ cmp.w #-24,d4 ; totaler Ueberlauf ?
+ blt.s fmul_zero ; ja, dann ist Ergebnis Null
+
+ swap d1 ; Form: 0fff ... ffff 0000 0000
+ beq.s fdiv_err ; falls Divisor Null, dann wird
+ lsl.l #7,d1 ; als Ergebnis unendlich ausgegeben
+ swap d3
+ beq.s fmul_zero ; falls Divident Null --> Ergebnis Null
+ lsl.l #7,d3
+
+fdiv_nlp: btst #30,d1 ; ist der Divisor normalisiert ?
+ bne.s fdiv_nor ; ja, -->
+ addq.w #1,d4 ; nein, Exponent erhoehen
+ lsl.l #1,d1 ; Divisor verschieben bis Form 01ff ..
+ bra.s fdiv_nlp
+
+fdiv_nor: clr.l d0 ; Ergebnis vorbesetzen
+ add.w #25,d4 ; Exponent ist nicht groesser als Null
+
+fdiv_lop: move.l d3,d5 ; Divident zwischenspeichern
+ sub.l d1,d3 ; Divisor abziehen
+ eori #extend,ccr ; X-Bit invertieren
+ bcc.s fdiv_one ; kein Carry: Divisor passt
+ move.l d5,d3 ; zurueckkopieren (X-Bit unveraendert!)
+fdiv_one: roxl.l #1,d0 ; Ergebnis aufbauen
+ lsl.l #1,d3 ; Divident verschieben
+ subq.w #1,d4 ; Exponent erniedrigen
+ beq.s fdiv_den ; falls Null, dann denormalisiert
+ btst #24,d0 ; fuehrende Eins in Ergebnis-Mantisse?
+ beq.s fdiv_lop ; nein, weiter rechnen
+
+fdiv_den: lsl.l #7,d0 ; Mantisse positionieren
+ beq fmul_res ; Null ?
+ bra fmul_rnd ; zum Runden
+;-----------------------------------------------------------------------------
+; Multiplikation mit einer Zweierpotenz: D0.S=D0.S * 2^(D1.W)
+
+ ds 0
+S_fmul2:
+ addq.l #1,_fmul_cnt.w
+ movem.l d1-d2,-(a7) ; Register retten
+ move.l d0,d2 ; Vorzeichen in D2 Bit 31 merken
+ lsl.l #1,d0 ; Vorzeichen rausschieben
+ beq.s fmul2_zero ; falls Null, dann ist Ergebnis Null
+ rol.l #8,d0 ; Form: ffff ... fff0 eeee eeee
+ clr.w d2 ; auf Wort vorbereiten
+ move.b d0,d2 ; Exponent in D2
+ beq.s fmul2_den
+ tst.w d1 ; Multiplikation oder Division?
+ bmi.s fmul2_div ; (neg. Exponent entspr. Div.)
+
+ add.w d1,d2 ; Summe der Exponenten bilden
+ cmp.w #maxexpo,d2 ; Ueberlauf?
+ bge.s fmul2_over ; ja, Ergebnis ist unendlich
+fmul2_res: move.b d2,d0 ; Ergebnisexponent einsetzen
+ ror.l #8,d0 ; Form: eeee eeee ffff ... fffx
+ roxl.l #1,d2 ; Vorzeichen ins X-Bit
+ roxr.l #1,d0 ; und ins Ergebnis einschieben
+fmul2_zero: movem.l (a7)+,d1-d2 ; Register restaurieren
+ rts
+
+fmul2_over: move.w #maxexpo,d2 ; Unendlich: Exponent = maxexpo
+ clr.l d0 ; Mantisse = Null
+ bra.s fmul2_res
+
+fmul2_div: add.w d1,d2 ; Summe der Exponenten bilden
+ bgt.s fmul2_res ; Unterlauf? nein --> Ergebnis
+ ori #Extend,ccr ; implizite Eins real machen
+ roxr.l #1,d0 ; Form: 1fff ... ffff xxxx xxxx
+fmul2_dnr: tst.w d2 ; Exponent = Null ?
+ beq.s fmul2_res ; ja, Ergebnis ist denormalisiert
+ lsr.l #1,d0 ; Mantisse denormalisieren
+ beq.s fmul2_zero ; totaler Unterlauf: Ergebnis ist Null
+ addq.w #1,d2 ; Exponent korrigieren
+ bra.s fmul2_dnr
+fmul2_ddd: add.w d1,d2 ; Summe der Exponenten bilden
+ bra.s fmul2_dnr ; mit denormalisiereter Eingabe bearbeiten
+
+fmul2_den: tst.w d1 ; Multiplikation oder Division
+ bmi.s fmul2_ddd
+ clr.b d0 ; Form: ffff ... fff0 0000 0000
+fmul2_nor: lsl.l #1,d0 ; Mantisse nach links schieben
+ bcs.s fmul2_stp ; bis fuehrende Eins auftaucht
+ subq.w #1,d1 ; oder zweiter Exponent Null wird
+ bne.s fmul2_nor
+ bra.s fmul2_res ; Ergebnis abliefern
+fmul2_stp: add.w d1,d2 ; Rest zum Exponenten addieren
+ bra.s fmul2_res ; Bias stimmt auch ( jetzt 127 statt 126)
+
+;-----------------------------------------------------------------------------
+; Vergleich zweier Zahlen: cmp d0,d1
+
+S_fcmp:
+ bclr #31,d0 ; Zahl 1 >=0 ?
+ bne.s fcmp_2
+fcmp_1:
+ bclr #31,d1 ; Zahl 2 >=0 ?
+ bne.s fcmp_12
+fcmp_11:
+ cmp.l d1,d0 ; beide Zahlen >=0
+ rts ; dann Betraege vergleichen
+fcmp_12:
+ moveq.l #1,d0 ; Zahl 1 >=0 und Zahl 2 <0
+ cmp.l #-1,d0
+ rts
+fcmp_2:
+ bclr #31,d1 ; Zahl 2 >=0 ?
+ bne.s fcmp_22
+fcmp_21:
+ moveq.l #-1,d0 ; Zahl 1 <0 und Zahl 2 >=0
+ cmp.w #1,d0 ; dann kleiner
+ rts
+fcmp_22:
+ neg.l d0
+ neg.l d1
+ cmp.l d1,d0 ; beide Zahlen <0, dann ver-
+ rts ; kehrtherum vergleichen
+
+
+;-----------------------------------------------------------------------------
+; Longint-->Gleitkomma
+; D0.L --> D0.S
+
+S_fitof:
+ movem.l d1-d2,-(a7) ; Register retten
+ tst.l d0 ; Integer ist Null ?
+ beq.s fitof_res; Ergebnis ist auch Null
+ smi d1 ; Vorzeichen in D1 merken
+ bpl.s fitof_pos
+ neg.l d0 ; ggf. Integer negieren
+fitof_pos: move.w #bias+32,d2 ; Exponent vorbesetzen
+fitof_shift: subq.w #1,d2 ; Mantisse verschieben
+ lsl.l #1,d0 ; bis fuehrende Eins rausfliegt
+ bcc.s fitof_shift
+ move.b d2,d0 ; Exponent einsetzen
+ ror.l #8,d0 ; Zahl positionieren
+ roxr.b #1,d1 ; Vorzeichen in X-Bit
+ roxr.l #1,d0 ; und ins Ergebnis
+fitof_res: movem.l (a7)+,d1-d2 ; fertig
+ rts
+
+;-----------------------------------------------------------------------------
+; Gleitkomma --> Longint:
+; D0.S --> D0.L
+
+S_fftoi:
+ movem.l d1-d2,-(a7) ; Register retten
+ roxl.l #1,d0 ; Vorzeichen in Carry
+ scs d1 ; in D1 merken
+ rol.l #8,d0 ; Form: ffff ... fffx eeee eeee
+ move.b d0,d2 ; Exponent extrahieren
+ sub.b #bias,d2 ; Bias subtrahieren
+ bmi.s fftoi_zero ; kleiner Null -> Ergebnis = Null
+ cmp.b #31,d2 ; Ueberlauf?
+ bge.s fftoi_over
+ ori #extend,ccr ; Implizite Eins explizit machen
+ roxr.l #1,d0
+ clr.b d0 ; Form: 1fff ... ffff 0000 0000
+fftoi_shft:
+ lsr.l #1,d0 ; jetzt Verschiebung bis
+ addq.b #1,d2 ; Exponent stimmt
+ cmp.b #31,d2
+ bne.s fftoi_shft
+ tst.b d1 ; Zahl negativ ?
+ bpl.s fftoi_pos
+ neg.l d0 ; ggf. Ergebnis negieren
+fftoi_pos:
+ movem.l (a7)+,d1-d2 ; Register wieder holen
+ rts
+fftoi_zero:
+ clr.l d0 ; Unterlauf; Ergebnis ist Null
+ bra.s fftoi_pos
+fftoi_over:
+ move.l #$7fffffff,d0 ; Ueberlauf: Maxint zurueckgeben
+ tst.b d1 ; positiv oder negativ ?
+ bpl.s fftoi_pos
+ not.l d0 ; Einser-Komplement erzeugt Minint
+ bra.s fftoi_pos
+
+;-----------------------------------------------------------------------------
+; Quadratwurzel : D0.S-->D0.S
+
+ ds 0
+fsqrt_domainerror:
+ move.l #$ffc00000,d0 ; -NAN zurueckgeben
+ movem.l (a7)+,d1-d4
+ rts
+fsqrt_sq0:
+ clr.l d0
+ movem.l (a7)+,d1-d4
+ rts
+S_fsqrt:
+ addq.l #1,_fsqrt_cnt.w
+ movem.l d1-d4,-(a7) ; D1-D4 werden sonst zerstoert
+ move.l d0,d4
+ bmi.s fsqrt_domainerror ; Fehler bei negativem Argument
+ swap d4 ; MSW des Arguments
+ and.l #$7f80,d4 ; Exponent isolieren
+ beq.s fsqrt_sq0 ; Zahl ist 0, wenn Exponent 0
+ and.l #$007fffff,d0 ; Mantisse isolieren
+ sub.w #$7f*$80,d4 ; Exponent im Zweierkomplement
+ bclr #7,d4 ; Exponent ungerade? (und LSB auf 0)
+ beq.s fsqrt_evenexp
+ add.l d0,d0 ; ja: Mantisse * 2
+ add.l #$01000000-$00800000,d0 ; Hidden Bit setzen, 1.Iteration
+
+fsqrt_evenexp:
+ ; 1. Iteration fuer geraden Exponenten: Hidden Bit nicht setzen
+ asr.w #1,d4 ; Exponent/2 mit Vorzeichen
+ add.w #$7f*$80,d4 ; Exponent wieder in Offset-Darst.
+ swap d4 ; neuen Exponenten im MSW aufheben
+ lsl.l #7,d0 ; x ausrichten
+ move.l #$40000000,d2 ; xroot nach erster Iteration
+ move.l #$10000000,d3 ; m2=2 << (MaxBit-1);
+fsqrt_loop10:
+ move.l d0,d1 ; xx2 = x
+fsqrt_loop11:
+ sub.l d2,d1 ; xx2 -= root
+ lsr.l #1,d2 ; xroot >>= 1
+ sub.l d3,d1 ; x2 -= m2
+ bmi.s fsqrt_dontset1
+ move.l d1,d0 ; x = xx2
+ or.l d3,d2 ; xroot += m2
+ lsr.l #2,d3 ; m2 >>= 2
+ bne.s fsqrt_loop11
+ bra.s fsqrt_d0d1same
+fsqrt_dontset1:
+ lsr.l #2,d3 ; m2 >>= 2
+ bne.s fsqrt_loop10 ; Schleife 15* abarbeiten
+ ; Bit 22..8
+ ; 17. Iteration (Bit 7) mit separatem Code durchfuehren:
+ move.l d0,d1 ; xx2 = x
+fsqrt_d0d1same:
+ sub.l d2,d1 ; xx2 -= root
+ ror.l #1,d2 ; xroot >>= 1 mitsamt Carry...
+ swap d2 ; auf neues Alignment umstellen
+ subq.l #1,d1 ; Carry von 0-0x4000: x2 -= m2
+ ; Teil 1
+ bmi.s fsqrt_dontset7
+ or.l #-$40000000,d1 ; 0 - 0x4000: x2 -= m2, Teil 2
+ move.l d1,d0 ; x = xx2
+ or.w #$4000,d2 ; xroot += m2
+fsqrt_dontset7:
+ swap d0 ; x auf neues Alignment umstellen
+
+ move.w #$1000,d3 ; m2 - Bit 16..31 bereits 0
+fsqrt_loop20:
+ move.l d0,d1 ; xx2 = x
+fsqrt_loop21:
+ sub.l d2,d1 ; xx2 -= xroot
+ lsr.l #1,d2 ; xroot >>= 1
+ sub.l d3,d1 ; x2 -= m2
+ bmi.s fsqrt_dontset2
+ move.l d1,d0 ; x = xx2
+ or.w d3,d2 ; xroot += m2
+ lsr.w #2,d3 ; m2 >>= 2
+ bne.s fsqrt_loop21
+ bra.s fsqrt_finish
+fsqrt_dontset2:
+ lsr.w #2,d3 ; m2 >>= 2
+ bne.s fsqrt_loop20 ; Schleife 7 * abarbeiten (n=6..0)
+fsqrt_finish:
+ sub.l d2,d0 ; Aufrunden notwendig ?
+ bls.s fsqrt_noinc
+ addq.l #1,d2 ; wenn ja, durchfuehren
+fsqrt_noinc:
+ bclr #23,d2 ; Hidden Bit loeschen
+ or.l d4,d2 ; Exponent und Mantisse kombinieren
+ move.l d2,d0 ; Ergebnis
+ movem.l (a7)+,d1-d4
+ rts ; Z-,S-, und V-Flag o.k.
+
+;-----------------------------------------------------------------------------
+; Absolutbetrag: D0.S--> D0.S
+
+ ds 0
+
+S_fabs: bclr #31,d0 ; ganz einfach...
+ rts
+
+;-----------------------------------------------------------------------------
+; Exponentialfunktion: D0.S--> D0.S
+
+; Die "krummen" Konstanten legen wir als hex ab, damit es keine Vergleichs-
+; fehler durch Rundungsvarianzen gibt.
+
+S_fexp_Const0: dc.l $3FB8AA3B ; ld(exp(1.0)) = ld(e) = 1/ln(2)
+S_fexp_ConstA: dc.l $3D0DF4E0 ; 0.034657359038 Polynomkonstanten
+S_fexp_ConstB: dc.l $411F4606 ; 9.9545957821
+S_fexp_ConstC: dc.l $441A7E3A ; 617.97226953
+S_fexp_ConstD: dc.l $42AED5C2 ; 87.417498202
+
+ ds 0
+S_fexp: movem.l d1-d5,-(sp)
+
+ bclr #31,d0 ; Vorzeichen loeschen und nach D2 retten
+ sne d2
+
+ move.l S_fexp_Const0(pc),d1 ; auf 2erpotenz umrechnen
+ bsr S_fmul
+
+ move.l d0,d3 ; in Ganzzahlanteil und Nach-
+ bsr S_fftoi ; kommastellen (z) zerlegen
+ move.l d0,d4 ; Ganzzahlanteil nach D4
+ bsr S_fitof
+ move.l d0,d1
+ move.l d3,d0
+ bsr S_fsub
+ move.l d0,d3
+
+ move.l d0,d1 ; z^2 berechnen
+ bsr S_fmul
+ move.l d0,d5 ; noch zu gebrauchen
+
+ move.l S_fexp_ConstD(pc),d1 ; --> D+z^2
+ bsr S_fadd
+ move.l d0,d1 ; --> C/(..)
+ move.l S_fexp_ConstC(pc),d0
+ bsr S_fdiv
+ move.l d0,d1 ; --> B-(..)
+ move.l S_fexp_ConstB(pc),d0
+ bsr S_fsub
+ move.l d3,d1 ; --> (..)-z
+ bsr S_fsub
+ exg d0,d5 ; Ergebnis retten
+ move.l S_fexp_ConstA(pc),d1 ; A*z^2 berechnen
+ bsr S_fmul
+ move.l d5,d1 ; ergibt Nenner
+ bsr S_fadd
+ move.l d0,d1 ; Quotient bilden
+ move.l d3,d0
+ bsr S_fdiv
+ moveq #1,d1 ; verdoppeln
+ bsr S_fmul2
+ move.l S_Const1(pc),d1 ; 1 addieren
+ bsr S_fadd
+ move.l d4,d1 ; Potenzieren
+ bsr S_fmul2
+
+ tst.b d2 ; war Argument negativ ?
+ beq.s S_fexp_ArgPos
+ move.l d0,d1 ; dann Kehrwert bilden
+ move.l S_Const1(pc),d0
+ bsr S_fdiv
+
+Terminate:
+S_fexp_ArgPos: movem.l (sp)+,d1-d5
+
+ rts
+
+;------------------------------------------------------------------------------
+; Sinus hyperbolicus: D0.S-->D0.S
+
+S_fsinh:
+ movem.l d1-d2,-(a7) ; Register retten
+ bsr S_fexp ; exp(x) berechnen
+ move.l d0,d2 ; in D2 merken
+ move.l d0,d1 ; exp(-x)=1/exp(x) berechnen
+ move.l #eins,d0
+ bsr S_fdiv
+ move.l d0,d1 ; Teilergebnisse subtrahieren
+ move.l d2,d0
+ bsr S_fsub
+ move.w #-1,d1 ; halbieren
+ bsr S_fmul2
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+
+;------------------------------------------------------------------------------
+; Cosinus hyperbolicus: D0.S-->D0.S
+
+S_fcosh:
+ movem.l d1-d2,-(a7) ; Register retten
+ bsr S_fexp ; exp(x) berechnen
+ move.l d0,d2 ; in D2 merken
+ move.l d0,d1 ; exp(-x)=1/exp(x) berechnen
+ move.l #eins,d0
+ bsr S_fdiv
+ move.l d2,d1 ; Teilergebnisse addieren
+ bsr S_fadd
+ move.w #-1,d1 ; halbieren
+ bsr S_fmul2
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Tangens hyperbolicus: D0.S-->D0.S
+
+S_ftanh:
+ movem.l d1-d3,-(a7) ; Register sichern
+ bsr S_fexp ; exp(x) berechnen
+ move.l d0,d2 ; in D2 merken
+ move.l d0,d1 ; exp(-x)=1/exp(x) berechnen
+ move.l #eins,d0
+ bsr S_fdiv
+ move.l d0,d3 ; in D3 merken
+ move.l d2,d1 ; Summe=Nenner berechnen
+ bsr S_fadd
+ exg d0,d2 ; jetzt exp(x) in D0, Nenner
+ move.l d3,d1 ; in D2
+ bsr S_fsub ; Zaehler berechnen
+ move.l d2,d1 ; Quotient berechnen
+ bsr S_fdiv
+ movem.l (a7)+,d1-d3 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Cotangens hyperbolicus: D0.S-->D0.S
+
+S_fcoth:
+ tst.l d0 ; Argument Null ?
+ beq.s S_fcoth_valerr ; dann zur Fehlerroutine
+ movem.l d1-d3,-(a7) ; Register sichern
+ bsr S_fexp ; exp(x) berechnen
+ move.l d0,d2 ; in D2 merken
+ move.l d0,d1 ; exp(-x)=1/exp(x) berechnen
+ move.l #eins,d0
+ bsr S_fdiv
+ move.l d0,d3 ; in D3 merken
+ move.l d0,d1 ; Differenz=Nenner berechnen
+ move.l d2,d0
+ bsr S_fsub
+ exg d0,d2 ; jetzt exp(x) in D0, Nenner
+ move.l d3,d1 ; in D2
+ bsr S_fadd ; Zaehler berechnen
+ move.l d2,d1 ; Quotient berechnen
+ bsr S_fdiv
+ movem.l (a7)+,d1-d3 ; Register zurueck
+ rts
+S_fcoth_valerr:
+ move.l #$7f800000,d0 ; +INF zurueckgeben
+ rts
+
+;-----------------------------------------------------------------------------
+; nat. Logarithmus: D0.S-->D0.S
+
+ ds 0
+S_fln:
+ tst.l d0 ; Argument <=0 ?
+ ble S_fln_errval
+ movem.l d1-d7,-(a7) ; Register retten
+ move.l d0,d3 ; Argument sichern
+
+ move.l #eins,d1 ; Zahl>1?
+ bsr S_fsub ; ( dies ist sinnvoll bei
+ tst.l d0 ; Zahlen <<1 );
+ smi d7 ; und die Vorzeichenumkehr merken
+ bpl.s S_fln_gr1 ; ja-->o.k.
+ move.l d3,d1 ; ansonsten Kehrwert bilden
+ move.l #eins,d0
+ bsr S_fdiv
+ move.l d0,d3
+
+S_fln_gr1:
+ clr.l d2 ; Merker = Null
+S_fln_nrm:
+ move.l d3,d0 ; Zahl > 1 ?
+ move.l #eins,d1
+ bsr S_fsub
+ bmi.s S_fln_isok
+ beq.s S_fln_isok
+ sub.l #$00800000,d3 ; ja-->Zahl durch 2 teilen...
+ addq.w #1,d2 ; ...und Merker erhoehen
+ bra.s S_fln_nrm ; nochmal probieren
+S_fln_isok:
+ move.l d0,d3 ; Zahl um Eins erniedrigt abspeichern
+ move.l d0,d4 ; yz:=y
+ moveq.l #1,d6 ; zaehler:=1
+ clr.l d5 ; Summe:=0
+ bchg #31,d3 ; Multiplikator negativ
+S_fln_loop:
+ move.l d6,d0 ; Zaehler in Gleitkomma wandeln
+ bsr S_fitof
+ move.l d0,d1 ; s:=s+yz/zaehler*vz
+ move.l d4,d0
+ bsr S_fdiv
+ move.l d5,d1
+ bsr S_fadd
+ cmp.l d5,d0 ; noch signifikant ?
+ beq.s S_fln_loopend
+ move.l d0,d5
+ addq.w #1,d6 ; zaehler:=zaehler+1
+ cmp.w #10,d6 ; Schleife fertig ?
+ beq.s S_fln_loopend
+ move.l d4,d0 ; yz:=yz*y
+ move.l d3,d1
+ bsr S_fmul
+ move.l d0,d4
+ bra.s S_fln_loop
+S_fln_loopend:
+ move.l d2,d0 ; Merker in Gleitkomma
+ bsr S_fitof
+ move.l #ln2,d1 ; * ln(2)
+ bsr S_fmul
+ move.l d5,d1 ; s:=s+merker
+ bsr S_fadd
+
+ tst.b d7 ; noch Vorzeichen tauschen ?
+ beq.s S_fln_end
+ bchg #31,d0
+S_fln_end:
+ movem.l (a7)+,d1-d7 ; Register zurueck
+ rts
+S_fln_errval:
+ move.l #$ffc00000,d0 ; -NAN zurueckgeben
+ rts
+
+;-----------------------------------------------------------------------------
+; 10er-Logarithmus : D0.S --> D0.S
+
+S_flog:
+ tst.l d0 ; Argument <=0 ?
+ ble.s S_flog_errval
+ bsr S_fln ; nat. Logarithmus bilden
+ move.l #ln10,d1 ; umrechnen
+ bsr S_fdiv
+ rts
+S_flog_errval:
+ move.l #$ffc00000,d0 ; -NAN zurueckgeben
+ rts
+
+;-----------------------------------------------------------------------------
+; Areasinus hyperbolicus: D0.S-->D0.S == ln[x+sqrt(x*x+1)]
+
+S_fasinh:
+ movem.l d1-d2,-(a7)
+ move.l d0,d2 ; Argument sichern
+ move.l d0,d1 ; quadrieren
+ bsr S_fmul
+ move.l #eins,d1 ; 1 addieren
+ bsr S_fadd
+ bsr S_fsqrt ; Wurzel ziehen
+ move.l d2,d1 ; Argument addieren
+ bsr S_fadd
+ bsr S_fln ; Logarithmus des ganzen
+ movem.l (a7)+,d1-d2
+ rts
+
+;-----------------------------------------------------------------------------
+; Areacosinus hyperbolicus: D0.S-->D0.S == ln[x+sqrt(x*x-1)]
+
+S_facosh:
+ movem.l d1-d2,-(a7) ; Register sichern
+ move.l d0,d2 ; Argument sichern
+ move.l #eins,d1 ; Argument <1 ?
+ bsr S_fcmp
+ bmi.s S_facosh_errval
+ move.l d2,d0 ; Argument zurueck
+ move.l d0,d1 ; quadrieren
+ bsr S_fmul
+ move.l #eins,d1 ; 1 abziehen
+ bsr S_fsub
+ bsr S_fsqrt ; Wurzel ziehen
+ move.l d2,d1 ; Argument addieren
+ bsr S_fadd
+ bsr S_fln ; Logarithmus des ganzen
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_facosh_errval:
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ move.l #$ffc00000,d0 ; NAN zurueckgeben
+ rts
+
+;-----------------------------------------------------------------------------
+; Areatangens hyperbolicus: D0.S-->D0.S == 0.5*ln((1+x)/(1-x))
+
+S_fatanh:
+ movem.l d1-d2,-(a7) ; Register sichern
+ move.l d0,d2 ; Argument sichern
+ bclr #31,d0 ; Vorzeichen weg
+ cmp.l #eins,d0
+ beq.s S_fatanh_inf ; =1-->INF
+ bhi.s S_fatanh_errval ; >1-->NAN
+ move.l d2,d1 ; Nenner berechnen
+ move.l #eins,d0
+ bsr S_fsub
+ exg d0,d2 ; Zaehler berechnen
+ move.l #eins,d1
+ bsr S_fadd
+ move.l d2,d1 ; Quotient daraus
+ bsr S_fdiv
+ bsr S_fln ; logarithmieren
+ move.w #-1,d1 ; halbieren
+ bsr S_fmul2
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_fatanh_inf:
+ move.l #$ff000000,d0 ; vorzeichenbehaftete Unend-
+ roxr.l #1,d0 ; lichkeit
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_fatanh_errval:
+ move.l #$7fc00000,d0 ; NAN geben
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Areakotangens hyperbolicus: D0.S--> D0.S == 0.5*ln((1+x)/(x-1))
+
+S_facoth:
+ movem.l d1-d2,-(a7) ; Register sichern
+ move.l d0,d2 ; Argument sichern
+ roxl.l #1,d0 ; Vorzeichen in X-Flag
+ cmp.l #eins*2,d0
+ beq.s S_facoth_inf ; =1-->INF
+ bmi.s S_facoth_errval ; <1-->NAN
+ move.l d2,d0 ; Nenner berechnen
+ move.l #eins,d1
+ bsr S_fsub
+ exg d0,d2 ; Zaehler berechnen
+ move.l #eins,d1
+ bsr S_fadd
+ move.l d2,d1 ; Quotient daraus
+ bsr S_fdiv
+ bsr S_fln ; logarithmieren
+ move.w #-1,d1 ; halbieren
+ bsr S_fmul2
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_facoth_inf:
+ move.l #$ff000000,d0 ; vorzeichenbehaftete Unend-
+ roxr.l #1,d0 ; lichkeit
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_facoth_errval:
+ move.l #$7fc00000,d0 ; NAN geben
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Kosinusfunktion: D0.S--> D0.S
+
+ ds 0
+S_fcos:
+ movem.l d1-d6,-(a7) ; Register retten
+ bclr #31,d0 ; cos(-x)=cos(x)
+
+ move.l #pi2,d1 ; auf Bereich 0..2*Pi reduzieren
+S_fcos_subtr:
+ cmp.l d1,d0 ; x>=2*Pi ?
+ blo.s S_fcos_subend ; ja-->Ende
+ bchg #31,d1 ; fuer Subtraktion
+ bsr S_fadd ; reduzieren
+ bchg #31,d1 ; Subtrahend wieder richtig
+ bra.s S_fcos_subtr
+S_fcos_subend:
+ cmp.l #pi,d0 ; x>Pi ?
+ blo.s S_fcos_nosub
+ exg d0,d1 ; ja-->cos(x)=cos(2*Pi-x)
+ bsr S_fsub
+S_fcos_nosub:
+ move.l d0,d1 ; wir brauchen nur x^2
+ bsr S_fmul
+ bset #31,d0
+ move.l d0,d3 ; -x^2 in D3
+ move.l d0,d4 ; D4 enthaelt laufende Potenz von x^2
+ ; inkl. Vorzeichen
+ move.l #zwei,d5 ; D5 enthaelt laufende Fakultaet
+ move.l #eins,d2 ; D2 enthaelt Summe
+ moveq.l #2,d6 ; D6 enthaelt Zaehler
+S_fcos_loop:
+ move.l d5,d1 ; s:=s+yz/zaehler
+ move.l d4,d0
+ bsr S_fdiv
+ move.l d2,d1
+ bsr S_fadd
+ cmp.l d2,d0 ; Veraendert sich Summe noch ?
+ beq.s S_fcos_end
+ move.l d0,d2
+ addq.b #2,d6 ; i:=i+1
+ cmp.b #22,d6 ; i=11 ?
+ beq.s S_fcos_end
+ move.w d6,d0 ; Fakultaet erhhen: *(2n-1)*(2n)
+ mulu.w d6,d0 ; =4*n^2-2*n
+ sub.w d6,d0
+ bsr S_fitof ; dazumultiplizieren
+ move.l d5,d1
+ bsr S_fmul
+ move.l d0,d5
+ move.l d4,d0 ; yz:=yz*y
+ move.l d3,d1
+ bsr S_fmul
+ move.l d0,d4
+ bra.s S_fcos_loop
+S_fcos_end:
+ ; Ergebnis bereits in D0
+ movem.l (a7)+,d1-d6 ; Register zurueck
+ rts
+
+;----------------------------------------------------------------------------
+; Sinus : D0.S-->D0.S
+
+S_fsin:
+ move.l d1,-(a7) ; Register retten
+ move.l #pihalf,d1 ; sin(x)=cos(x-pi/2)
+ bsr S_fsub
+ bsr S_fcos
+ move.l (a7)+,d1 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Tangens: D0.S-->D0.S
+
+S_ftan:
+ movem.l d1-d4,-(a7) ; Register retten
+ tst.l d0 ; Vorzeichen merken
+ smi d4
+ bclr #31,d0
+ move.l #pi,d1 ; auf Bereich 0..Pi reduzieren
+S_ftan_subtr:
+ cmp.l d1,d0 ; x>=Pi ?
+ blo.s S_ftan_subend ; ja-->Ende
+ bchg #31,d1 ; fuer Subtraktion
+ bsr S_fadd ; reduzieren
+ bchg #31,d1 ; Subtrahend wieder richtig
+ bra.s S_ftan_subtr
+S_ftan_subend:
+ move.l d0,d2 ; Argument merken
+ bsr S_fcos ; Nenner rechnen
+ move.l d0,d3 ; Nenner merken
+ move.l d0,d1 ; sqr(1-x^2) rechnen
+ bsr S_fmul
+ move.l d0,d1
+ move.l #eins,d0
+ bsr S_fsub
+ bsr S_fsqrt
+ move.l d3,d1
+ bsr S_fdiv ; Quotient
+ tst.b d4 ; Vorzeichen dazu
+ beq.s S_ftan_noneg
+ bchg #31,d0
+S_ftan_noneg:
+ movem.l (a7)+,d1-d4 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Kotangens: D0.S-->D0.S
+
+S_fcot:
+ movem.l d1-d4,-(a7) ; Register retten
+ tst.l d0 ; Vorzeichen merken
+ smi d4
+ bclr #31,d0
+ move.l #pi,d1 ; auf Bereich 0..Pi reduzieren
+S_fcot_subtr:
+ cmp.l d1,d0 ; x>=Pi ?
+ blo.s S_fcot_subend ; ja-->Ende
+ bchg #31,d1 ; fuer Subtraktion
+ bsr S_fadd ; reduzieren
+ bchg #31,d1 ; Subtrahend wieder richtig
+ bra.s S_fcot_subtr
+S_fcot_subend:
+ move.l d0,d2 ; Argument merken
+ bsr S_fcos ; Zaehler rechnen
+ move.l d0,d3 ; Zaehler merken
+ move.l d0,d1 ; sqr(1-x^2) rechnen
+ bsr S_fmul
+ move.l d0,d1
+ move.l #eins,d0
+ bsr S_fsub
+ bsr S_fsqrt
+ move.l d0,d1
+ move.l d3,d0
+ bsr S_fdiv ; Quotient
+ tst.b d4 ; Vorzeichen dazu
+ beq.s S_fcot_noneg
+ bchg #31,d0
+S_fcot_noneg:
+ movem.l (a7)+,d1-d4 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Arcustangens: D0.S-->D0.S
+
+S_fatan:
+ movem.l d1-d6,-(a7) ; Register sichern
+ subq.l #2,a7 ; Platz fuer Hilfsvariablen
+ tst.l d0 ; Vorzeichen merken...
+ smi (a7)
+ bclr #31,d0 ; ...und loeschen
+ cmp.l #eins,d0 ; Argument>1 ?
+ shi 1(a7) ; ja :
+ bls.s S_fatan_no1 ; nein :
+ move.l d0,d1 ; ja : Kehrwert bilden
+ move.l #eins,d0
+ bsr S_fdiv
+S_fatan_no1:
+ move.l #3,d2 ; Zaehler initialisieren
+ move.l d0,d5 ; Summe initialisieren
+ move.l d0,d4 ; laufende Potenz = x^1
+ move.l d0,d1 ; x^2 berechnen
+ bsr S_fmul
+ move.l d0,d3
+ bset #31,d3 ; immer mit -x^2 multiplizieren
+S_fatan_loop:
+ move.l d4,d0 ; naechste Potenz ausrechnen
+ move.l d3,d1
+ bsr S_fmul
+ move.l d0,d4
+ move.l d2,d0 ; Nenner in Gleitkomma
+ bsr S_fitof
+ move.l d0,d1 ; Division ausfuehren
+ move.l d4,d0
+ bsr S_fdiv
+ move.l d5,d1 ; zur Summe addieren
+ bsr S_fadd
+ cmp.l d0,d5 ; noch signifikant ?
+ beq.s S_fatan_endloop ; nein-->Ende
+ move.l d0,d5
+ addq.l #2,d2 ; Zaehler erhoehen
+ cmp.l #61,d2 ; fertig ?
+ bne.s S_fatan_loop
+S_fatan_endloop:
+ move.l d5,d0 ; Ergebnis in D0 bringen
+ tst.b 1(a7) ; war Argument < 1 ?
+ beq.s S_fatan_not2
+ bchg #31,d0 ; ja : Erg.=Pi/2-Erg
+ move.l #pihalf,d1
+ bsr S_fadd
+S_fatan_not2:
+ tst.b (a7) ; Vorzeichen dazu
+ beq.s S_fatan_not1
+ bset #31,d0
+S_fatan_not1:
+ addq.l #2,a7 ; Hilfsvar. abraeumen
+ movem.l (a7)+,d1-d6 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Arcuskotangens: D0.S-->D0.S
+
+S_facot:
+ move.l d1,-(a7) ; Register sichern
+ bsr S_fatan ; acot(x)=pi/2-atan(x)
+ bchg #31,d0
+ move.l #pihalf,d1
+ bsr S_fadd
+ move.l (a7)+,d1 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Arcussinus: D0.S --> D0.S
+
+S_fasin:
+ movem.l d1-d2,-(a7) ; Register retten
+ move.l d0,d2 ; Argument merken
+ move.l d0,d1 ; Quadrat berechnen
+ bsr S_fmul
+ bset #31,d0 ; 1-x^2 bilden
+ move.l #eins,d1
+ bsr S_fadd
+ tst.l d0 ; Sonderfaelle abfangen
+ bmi.s S_fasin_errval ; <0 geht nicht
+ beq.s S_fasin_inf ; Endwerte
+ bsr S_fsqrt ; Wurzel ...
+ move.l d0,d1 ; ... und Quotient
+ move.l d2,d0
+ bsr S_fdiv
+ bsr S_fatan ; zuletzt das wichtigste
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_fasin_errval:
+ move.l #$7fc00000,d0 ; NAN zurueckgeben
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_fasin_inf:
+ move.l #pihalf,d0 ; +- pi/2 zurueckgeben
+ and.l #$80000000,d2 ; Vorzeichen dazu
+ or.l d2,d0
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+
+;-----------------------------------------------------------------------------
+; Arcuskosinus: D0.S --> D0.S
+
+S_facos:
+ tst.l d0 ; Argument=0 ?
+ beq.s S_facos_inf
+ move.l d0,d2 ; Argument merken
+ move.l d0,d1 ; Quadrat berechnen
+ bsr S_fmul
+ bset #31,d0 ; 1-x^2 bilden
+ move.l #eins,d1
+ bsr S_fadd
+ tst.l d0 ; Sonderfaelle abfangen
+ bmi.s S_facos_errval ; <0 geht nicht
+ bsr S_fsqrt ; Wurzel ...
+ move.l d2,d1 ; ... und Quotient
+ bsr S_fdiv
+ bsr S_fatan ; zuletzt das wichtigste
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_facos_errval:
+ move.l #$7fc00000,d0 ; NAN zurueckgeben
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+S_facos_inf:
+ move.l #pihalf,d0 ; +- pi/2 zurueckgeben
+ and.l #$80000000,d2 ; Vorzeichen dazu
+ or.l d2,d0
+ movem.l (a7)+,d1-d2 ; Register zurueck
+ rts
+
+S_floatlibend:
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-20 10:41:51 +0000