; 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: