/* libgcc1 routines for 68000 w/o floating-point hardware. Copyright (C) 1994, 1996, 1997, 1998 Free Software Foundation, Inc. This file is part of GNU CC. GNU CC is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. In addition to the permissions in the GNU General Public License, the Free Software Foundation gives you unlimited permission to link the compiled version of this file with other programs, and to distribute those programs without any restriction coming from the use of this file. (The General Public License restrictions do apply in other respects; for example, they cover modification of the file, and distribution when not linked into another program.) This file is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; see the file COPYING. If not, write to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /* As a special exception, if you link this library with files compiled with GCC to produce an executable, this does not cause the resulting executable to be covered by the GNU General Public License. This exception does not however invalidate any other reasons why the executable file might be covered by the GNU General Public License. */ /* Use this one for any 680x0; assumes no floating point hardware. The trailing " '" appearing on some lines is for ANSI preprocessors. Yuk. Some of this code comes from MINIX, via the folks at ericsson. D. V. Henkel-Wallace (gumby@cygnus.com) Fete Bastille, 1992 */ /* Floating point support largely rewritten to improve compliance with IEEE 754-1985, and improve speed on Coldfire family processors. Wayne Deeter (wrd@deetour.net) November-December 1999 */ /* These are predefined by new versions of GNU cpp. */ #ifndef __USER_LABEL_PREFIX__ #define __USER_LABEL_PREFIX__ _ #endif #ifndef __REGISTER_PREFIX__ #define __REGISTER_PREFIX__ #endif #ifndef __IMMEDIATE_PREFIX__ #define __IMMEDIATE_PREFIX__ # #endif /* ANSI concatenation macros. */ #define CONCAT1(a, b) CONCAT2(a, b) #define CONCAT2(a, b) a ## b /* Use the right prefix for global labels. */ #define SYM(x) CONCAT1 (__USER_LABEL_PREFIX__, x) /* Use the right prefix for registers. */ #define REG(x) CONCAT1 (__REGISTER_PREFIX__, x) /* Use the right prefix for immediate values. */ #define IMM(x) CONCAT1 (__IMMEDIATE_PREFIX__, x) #define d0 REG (d0) #define d1 REG (d1) #define d2 REG (d2) #define d3 REG (d3) #define d4 REG (d4) #define d5 REG (d5) #define d6 REG (d6) #define d7 REG (d7) #define a0 REG (a0) #define a1 REG (a1) #define a2 REG (a2) #define a3 REG (a3) #define a4 REG (a4) #define a5 REG (a5) #define a6 REG (a6) #define fp REG (fp) #define sp REG (sp) #if defined(L_floatex) || defined(L_double) || defined(L_float) || defined(L_extended) || defined(L_cvt) | This is an attempt at a decent floating point (single, double and | extended double) code for the GNU C compiler. It should be easy to | adapt to other compilers (but beware of the local labels!). | Starting date: 21 October, 1990 | It is convenient to introduce the notation (s,e,f) for a floating point | number, where s=sign, e=exponent, f=fraction. We will call a floating | point number fpn to abbreviate, independently of the precision. | Let MAX_EXP be in each case the maximum exponent (255 for floats, 1023 | for doubles and 16383 for long doubles). We then have the following | different cases: | 1. Normalized fpns have 0 < e < MAX_EXP. They correspond to | (-1)^s x 1.f x 2^(e-bias-1). | 2. Denormalized fpns have e=0. They correspond to numbers of the form | (-1)^s x 0.f x 2^(-bias). | 3. +/-INFINITY have e=MAX_EXP, f=0. | 4. Quiet NaN (Not a Number) have all bits set. | 5. Signaling NaN (Not a Number) have s=0, e=MAX_EXP, f=1. |============================================================================= | exceptions |============================================================================= | This is the floating point condition code register (_fpCCR): | | struct { | int _exception_bits; | int _trap_enable_bits; | int _sticky_bits; | int _rounding_mode; | unsigned int _unused :16; | unsigned int _op2 :3; | unsigned int _dst :3; | unsigned int _op1 :3; | unsigned int _last_operation : 7; | union { | float sf; | double df; | long double xf; | } _operand1; | union { | float sf; | double df; | long double xf; | } _operand2; | union { | int si; | float sf; | double df; | long double xf; | } _result; | } _fpCCR; .globl SYM (_fpCCR) .globl $_exception_handler Q_NaN = 0xffffffff SF_MAX_EXP = 0xff SF_BIAS = 127 SF_SIGNIF_DIG = 24 SF_INFINITY = 0x7f800000 DF_MAX_EXP = 0x7ff DF_BIAS = 1023 DF_SIGNIF_DIG = 53 DF_INFINITY = 0x7ff00000 XF_MAX_EXP = 0x7fff XF_BIAS = 16383 XF_SIGNIF_DIG = 64 XF_INFINITY = 0x7fff0000 | | Many routines work with reduced exponent biases as the sf & df | formats treat a zero exponent as a special case for denormalized | numbers. When a number is broken into its parts, and has its | possibly hidden bit restored, it makes sense to do it this way | as the denorms have the same exponent as the smallest normalized | numbers. They just need a special exponent value when they're ' | packed. SF_RBIAS = SF_BIAS-1 DF_RBIAS = DF_BIAS-1 | Offsets: EBITS = 0 TRAPE = 4 STICK = 8 ROUND = 12 LASTO = 16 OPER1 = 20 OPER2 = 32 RESLT = 44 | The following exception types are supported: INEXACT_RESULT = 0x0001 UNDERFLOW = 0x0002 OVERFLOW = 0x0004 DIVIDE_BY_ZERO = 0x0008 INVALID_OPERATION = 0x0010 | The allowed rounding modes are: UNKNOWN = -1 ROUND_TO_NEAREST = 0 | round result to nearest representable value ROUND_TO_ZERO = 1 | round result towards zero ROUND_TO_PLUS = 2 | round result towards plus infinity ROUND_TO_MINUS = 3 | round result towards minus infinity | Formats: NIL = 0 SINGLE_FLOAT = 1 DOUBLE_FLOAT = 2 LONG_FLOAT = 3 INT = 5 DST_FM = 1024 OP1_FM = 128 OP2_FM = 8192 DST_NIL = DST_FM*NIL DST_FLOAT = DST_FM*SINGLE_FLOAT DST_DOUBLE = DST_FM*DOUBLE_FLOAT DST_LDOUBLE = DST_FM*LONG_FLOAT DST_INT = DST_FM*INT OP1_NIL = OP1_FM*NIL OP1_FLOAT = OP1_FM*SINGLE_FLOAT OP1_DOUBLE = OP1_FM*DOUBLE_FLOAT OP1_LDOUBLE = OP1_FM*LONG_FLOAT OP1_INT = OP1_FM*INT OP2_NIL = OP2_FM*NIL OP2_FLOAT = OP2_FM*SINGLE_FLOAT OP2_DOUBLE = OP2_FM*DOUBLE_FLOAT OP2_LDOUBLE = OP2_FM*LONG_FLOAT OP2_INT = OP2_FM*INT | The allowed values for the last operation are: NOOP = 0 ADD = 1 SUB = 2 MULTIPLY = 3 DIVIDE = 4 NEGATE = 5 COMPARE = 6 EXTENDSFDF = 7 TRUNCDFSF = 8 FIXDFSI = 9 FIXSFSI = 10 EXTENDDFXF = 11 TRUNCXFDF = 12 EXTENDSFXF = 13 TRUNCXFSF = 14 FIXXFSI = 15 FLOATSISF = 16 #endif /* L_floatex || L_double || L_float || L_extended || L_cvt */ #ifdef L_floatex .data .even SYM (_fpCCR): __exception_bits: .long 0 __trap_enable_bits: .long 0 __sticky_bits: .long 0 __rounding_mode: .long ROUND_TO_NEAREST __last_operation: .long NOOP __operand1: .long 0 .long 0 .long 0 __operand2: .long 0 .long 0 .long 0 __result: .long 0 .long 0 .long 0 |============================================================================= | __clear_sticky_bits |============================================================================= | The sticky bits are normally not cleared (thus the name), whereas the | exception type and exception value reflect the last computation. | This routine is provided to clear them (you can also write to _fpCCR, | since it is globally visible). .globl SYM (__clear_sticky_bit) .text .even | void __clear_sticky_bits(void); SYM (__clear_sticky_bit): lea SYM (_fpCCR),a0 clrl a0@(STICK) rts |============================================================================= | $_exception_handler |============================================================================= .text .even | This is the common exit point if an exception occurs. | NOTE: it is NOT callable from C! | It expects the exception type in d7, the format (SINGLE_FLOAT, | DOUBLE_FLOAT or LONG_FLOAT) in d6, and the last operation code in d5. | It sets the corresponding exception and sticky bits, and the format. | Depending on the format if fills the corresponding slots for the | operands which produced the exception (all this information is provided | so if you write your own exception handlers you have enough information | to deal with the problem). | Then checks to see if the corresponding exception is trap-enabled, | in which case it pushes the address of _fpCCR and traps through | trap FPTRAP (15 for the moment). #ifndef FPTRAP FPTRAP = 15 #endif $_exception_handler: lea SYM (_fpCCR),a1 orl d7,a1@(STICK) | Set sticky bits movel a1@(TRAPE),d4 | Any traps taken? andl d7,d4 jne 3f | Yep, go take it 9: movel IMM (7*DST_FM),d4 | Find destination format andl sp@,d4 cmpl IMM (DST_LDOUBLE),d4 | Result by pointer? jeq 2f | Yep, go store 1: moveml sp@(4),d2-d7 | Restore data registers unlk a6 | and return rts 2: movel d0,a0@ | Save away the result movel d1,a0@(4) movel d2,a0@(8) movel a0,d0 | & return the pointer jra 1b | We've got trap. Save away all the info & go do it ' 3: movel d7,a1@(EBITS) | Set __exception_bits movel d0,a1@(RESLT) | Save the result movel d1,a1@(RESLT+4) movel d2,a1@(RESLT+8) movel a0,d7 | Save a0 (may be result ptr) lea a6@(8),a0 | Point to operands movel sp@,d3 | Get operation info movel d3,a1@(LASTO) | and save movel IMM (3),d4 lsrl IMM (7),d3 andl d3,d4 | Get OP1 format jeq 2f | None? lea a1@(OPER1),a1 1: movel a0@+,a1@+ | Copy words ... subl IMM (1),d4 jne 1b 2: movel IMM (3),d4 lsrl IMM (6),d3 andl d3,d4 | Get OP2 format jeq 2f lea SYM (_fpCCR)+OPER2,a1 1: movel a0@+,a1@+ subl IMM (1),d4 jne 1b 2: pea SYM (_fpCCR) | Push address of _fpCCR trap IMM (FPTRAP) | and trap lea sp@(4),sp lea SYM (_fpCCR)+RESLT,a1 movel a1@+,d0 | Restore result movel a1@+,d1 movel a1@+,d2 movel d7,a0 | and pointer jra 9b | and return #endif /* L_floatex */ #ifdef L_mulsi3 .text .proc .globl SYM (__mulsi3) SYM (__mulsi3): movew sp@(4), d0 /* x0 -> d0 */ muluw sp@(10), d0 /* x0*y1 */ movew sp@(6), d1 /* x1 -> d1 */ muluw sp@(8), d1 /* x1*y0 */ #ifndef __mcf5200__ addw d1, d0 #else addl d1, d0 #endif swap d0 clrw d0 movew sp@(6), d1 /* x1 -> d1 */ muluw sp@(10), d1 /* x1*y1 */ addl d1, d0 rts #endif /* L_mulsi3 */ #ifdef L_udivsi3 .text .proc .globl SYM (__udivsi3) SYM (__udivsi3): #ifndef __mcf5200__ movel d2, sp@- movel sp@(12), d1 /* d1 = divisor */ movel sp@(8), d0 /* d0 = dividend */ cmpl IMM (0x10000), d1 /* divisor >= 2 ^ 16 ? */ jcc L3 /* then try next algorithm */ movel d0, d2 clrw d2 swap d2 divu d1, d2 /* high quotient in lower word */ movew d2, d0 /* save high quotient */ swap d0 movew sp@(10), d2 /* get low dividend + high rest */ divu d1, d2 /* low quotient */ movew d2, d0 jra L6 L3: movel d1, d2 /* use d2 as divisor backup */ L4: lsrl IMM (1), d1 /* shift divisor */ lsrl IMM (1), d0 /* shift dividend */ cmpl IMM (0x10000), d1 /* still divisor >= 2 ^ 16 ? */ jcc L4 divu d1, d0 /* now we have 16 bit divisor */ andl IMM (0xffff), d0 /* mask out divisor, ignore remainder */ /* Multiply the 16 bit tentative quotient with the 32 bit divisor. Because of the operand ranges, this might give a 33 bit product. If this product is greater than the dividend, the tentative quotient was too large. */ movel d2, d1 mulu d0, d1 /* low part, 32 bits */ swap d2 mulu d0, d2 /* high part, at most 17 bits */ swap d2 /* align high part with low part */ tstw d2 /* high part 17 bits? */ jne L5 /* if 17 bits, quotient was too large */ addl d2, d1 /* add parts */ jcs L5 /* if sum is 33 bits, quotient was too large */ cmpl sp@(8), d1 /* compare the sum with the dividend */ jls L6 /* if sum > dividend, quotient was too large */ L5: subql IMM (1), d0 /* adjust quotient */ L6: movel sp@+, d2 rts #else /* __mcf5200__ */ /* Coldfire implementation of non-restoring division algorithm from Hennessy & Patterson, Appendix A. */ link a6,IMM (-12) moveml d2-d4,sp@ movel a6@(8),d0 movel a6@(12),d1 clrl d2 | clear p moveq IMM (31),d4 L1: addl d0,d0 | shift reg pair (p,a) one bit left addxl d2,d2 movl d2,d3 | subtract b from p, store in tmp. subl d1,d3 jcs L2 | if no carry, bset IMM (0),d0 | set the low order bit of a to 1, movl d3,d2 | and store tmp in p. L2: subql IMM (1),d4 jcc L1 moveml sp@,d2-d4 | restore data registers unlk a6 | and return rts #endif /* __mcf5200__ */ #endif /* L_udivsi3 */ #ifdef L_divsi3 .text .proc .globl SYM (__divsi3) SYM (__divsi3): movel d2, sp@- moveq IMM (1), d2 /* sign of result stored in d2 (=1 or =-1) */ movel sp@(12), d1 /* d1 = divisor */ jpl L1 negl d1 #ifndef __mcf5200__ negb d2 /* change sign because divisor <0 */ #else negl d2 /* change sign because divisor <0 */ #endif L1: movel sp@(8), d0 /* d0 = dividend */ jpl L2 negl d0 #ifndef __mcf5200__ negb d2 #else negl d2 #endif L2: movel d1, sp@- movel d0, sp@- jbsr SYM (__udivsi3) /* divide abs(dividend) by abs(divisor) */ addql IMM (8), sp tstb d2 jpl L3 negl d0 L3: movel sp@+, d2 rts #endif /* L_divsi3 */ #ifdef L_umodsi3 .text .proc .globl SYM (__umodsi3) SYM (__umodsi3): movel sp@(8), d1 /* d1 = divisor */ movel sp@(4), d0 /* d0 = dividend */ movel d1, sp@- movel d0, sp@- jbsr SYM (__udivsi3) addql IMM (8), sp movel sp@(8), d1 /* d1 = divisor */ #ifndef __mcf5200__ movel d1, sp@- movel d0, sp@- jbsr SYM (__mulsi3) /* d0 = (a/b)*b */ addql IMM (8), sp #else mulsl d1,d0 #endif movel sp@(4), d1 /* d1 = dividend */ subl d0, d1 /* d1 = a - (a/b)*b */ movel d1, d0 rts #endif /* L_umodsi3 */ #ifdef L_modsi3 .text .proc .globl SYM (__modsi3) SYM (__modsi3): movel sp@(8), d1 /* d1 = divisor */ movel sp@(4), d0 /* d0 = dividend */ movel d1, sp@- movel d0, sp@- jbsr SYM (__divsi3) addql IMM (8), sp movel sp@(8), d1 /* d1 = divisor */ #ifndef __mcf5200__ movel d1, sp@- movel d0, sp@- jbsr SYM (__mulsi3) /* d0 = (a/b)*b */ addql IMM (8), sp #else mulsl d1,d0 #endif movel sp@(4), d1 /* d1 = dividend */ subl d0, d1 /* d1 = a - (a/b)*b */ movel d1, d0 rts #endif /* L_modsi3 */ #ifdef L_double | Entry points: .globl SYM (__adddf3) .globl SYM (__subdf3) .globl SYM (__muldf3) .globl SYM (__divdf3) .globl SYM (__negdf2) .globl SYM (__cmpdf2) .text .even |============================================================================= |============================================================================= | double precision routines |============================================================================= |============================================================================= | A double precision floating point number (double) has the format: | | struct _double { | unsigned int sign : 1; /* sign bit */ | unsigned int exponent : 11; /* exponent, biased by 1023 */ | unsigned int fraction : 52; /* fraction */ | } double; | | Thus sizeof(double) = 8 (64 bits). | | All the routines are callable from C programs, and return the result | in the register pair d0-d1. They also preserve all registers except | d0-d1 and a0-a1. |============================================================================= | __subdf3 |============================================================================= | double __subdf3(double, double); SYM (__subdf3): link a6,IMM (-24) moveml d2-d7,sp@ pea SUB+OP1_DOUBLE+OP2_DOUBLE+DST_DOUBLE movel a6@(16),d5 | get second operand bchg IMM (31),d5 | change sign of second operand jra 1f | and fall through, so we always add |============================================================================= | __adddf3 |============================================================================= | double __adddf3(double, double); SYM (__adddf3): link a6,IMM (-24) moveml d2-d7,sp@ pea ADD+OP1_DOUBLE+OP2_DOUBLE+DST_DOUBLE movel a6@(16),d5 | get second operand 1: movel a6@(20),d6 | movel a6@(8),d0 | get first operand movel a6@(12),d1 | movel d0,a1 | Get a's sign bit in a1 ' movel d5,a0 | Get b's sign in a0 ' bclr IMM (31),d0 | Clear signs from working regs bclr IMM (31),d5 subl d0,a1 subl d5,a0 movel IMM (0xfff00000),d2 | Exponent mask & -hidden bit movel d0,d4 | Get a's exponent in d4 ' andl d2,d4 addl d2,d4 | & adjust bias subl d4,d0 | Clear exponent & set hidden bit movel d5,d3 | Get b's exponent in d3 ' andl d2,d3 addl d2,d3 | & adjust bias subl d3,d5 | Clear exponent & set hidden bit swap d4 | Right justify exponents swap d3 lsrl IMM (4),d4 lsrl IMM (4),d3 movel IMM (0x7fe),d2 | Special case literal cmpl d2,d4 | Zero, denorm, INFINITY, NaN? jcc Ladddf$a$spec | Yep, go handle Ladddf$1: cmpl d2,d3 | Zero, denorm, INFINITY, NaN? jcc Ladddf$b$spec | Yep, go handle Ladddf$2: | Both operands are finite and non-zero. Here we shift the numbers until | their exponents are the same. subl d4,d3 | compare the exponents jeq Ladddf$3 | If equal we don't need to shift ' jmi 9f | Branch if second exponent is lower addl d3,d4 | Put largest exponent into d4 addl a0,a1 | Swap a & b addl a1,a0 addl a0,a1 movel d0,d2 | Swap high words movel d5,d0 movel d2,d5 movel d1,d2 | Swap low words movel d6,d1 movel d2,d6 jra 1f | Here we have a's exponent larger than b's, so we have to shift b. 9: negl d3 | Get positive right-shift count | If difference is too large we don't shift, just summarize b's | non-zeroness in d3 1: cmpl IMM (DF_SIGNIF_DIG+2),d3 jge Ladddf$b$small addl d4,a0 | Save exponent movel IMM (32), d4 cmpl d4,d3 | If difference >= 32, shift by longs jge 5f subl d3,d4 | Get left shift count movel d6,d2 lsll d4,d2 lsrl d3,d6 movel d5,d7 lsll d4,d7 orl d7,d6 lsrl d3,d5 movel IMM (0),d3 | Show zero lsb's ' movew a0,d4 | Restore exponent subl d4,a0 jra Ladddf$4 5: subl d4,d3 | Get right-shift count subl d3,d4 | and left-shift count movel d3,d7 movel d6,d3 lsll d4,d3 movel d6,d2 lsrl d7,d2 movel d5,d6 lsll d4,d5 orl d5,d2 lsrl d7,d6 movel IMM (0), d5 movew a0,d4 | Restore exponent subl d4,a0 jra Ladddf$4 | If one of the numbers is pretty small (difference of exponents >= | DF_SIGNIF_DIG+1) we just summarize it to show an Inexact operation. Ladddf$b$small: movel IMM (2),d2 | Show b is piddling movel IMM (0),d3 movel d3,d5 movel d3,d6 jra Ladddf$4 | Now go do operation & rounding Ladddf$3: movel IMM (0),d2 | Equal exponents; zap the extension | words; d3 is already zero Ladddf$4: | Now we have the numbers in d0-d1 and d5-d6/d2-d3, the exponent in d4, and | the signs in a0,a1. | Here we have to decide whether to add or subtract the numbers: cmpl a0,a1 | Compare signs jne Lsubdf$0 | If the signs are different we get | to subtract addl d6,d1 | else we add the numbers addxl d5,d0 | | Before rounding normalize so bit #DF_SIGNIF_DIG-1 is set (we consider | the case of denormalized numbers in the rounding routine itself). | As in the addition (not in the subtraction!) we could have set | one more bit we check this: btst IMM (DF_SIGNIF_DIG),d0 jeq Ladddf$unnormal movel IMM (31),d6 movel d2,d5 | We needn't bother about the exact ' lsll d6,d5 | value in d3 as we're already normalized; ' orl d5,d3 | as it may affect rounding, only lsrl IMM (1),d2 | its zero/nonzero nature matters. movel d1,d5 lsll d6,d5 orl d5,d2 lsrl IMM (1),d1 movel d0,d5 lsll d6,d5 orl d5,d1 lsrl IMM (1),d0 addl IMM (1),d4 jra Lround$common | Go finish up! Lsubdf$0: | Here we do the subtraction. negl d3 negxl d2 subxl d6,d1 | subxl d5,d0 | jeq Ladddf$ezer | if zero just exit jpl Ladddf$unnormal | if positive skip the following movel a0,a1 | Use b's sign ' negl d1 | and negate result negxl d0 | We don't need to touch d2,d3 as ' | the exponents are equal. d2,d3 | are therfore both zero | Normalization for add/subtract | When the exponents of the operands differ by more than one, there will | at most be one left shift. This means that whenever there is more than | one left shift, the result is exact, and rounding may be skipped. Ladddf$unnormal: movel IMM (0x100000),d5 | MSB position cmpl d5,d0 | Already normalized? jcc Lround$common lea SYM (_fpCCR),a0 btst IMM (1),a0@(TRAPE+3) | Check for Underflow trap enabled jne Ladddf$normal | Yep, don't check for denorm ' tstl d4 | Already denormalized? jle Ladddf$exact | Yep, don't shift ' | First we find the MSB of the value. Find the most significant non-zero | word. tstl d0 | In the first word? jeq 3f | Nope, go search second word | We needn't touch d3 here. If it is non-zero, only one left shift will ' | be required. If the exponents differ by at most one, this will ensure | that both d2 and d3 are zero, and the result is exact. addl d2,d2 | Perhaps shift last bit out addxl d1,d1 addxl d0,d0 subl IMM (1),d4 | Denormal yet? jeq Ladddf$exact | Yep, but at least it's exact! ' cmpl d5,d0 jcc Lround$common | or perhaps not? | Left shift until it's either normalized, or we've reached the denormal | exponent value 1: addl d1,d1 | Try moving over one place addxl d0,d0 subl IMM (1),d4 | Denormal yet? jeq Ladddf$exact | Yep, go adjust other words cmpl d5,d0 | Normalized yet? jcs 1b | Nope, try another place Ladddf$exact: lsll IMM (4),d4 | All we need to do now is swap d4 | assemble the double float! addl d4,d0 | Put in the exponent Ladddf$ret: addl a1,d0 | and the sign Ladddf$pzer: moveml sp@(4),d2-d7 | We're done! ' | XXX if frame pointer is ever removed, stack pointer must | be adjusted here. unlk a6 rts | d0 is zero, the exponents of the operands are at most one different. | d1 is zero only in one special case: a is 1.0*2^n and b is the next | smallest representable number, and the difference is taken. In this | case all result bits are zero except for bit 31 of d2. In any case, | when we get here, d3 is zero! 3: movel IMM (20),d3 | Left shift count - 1 subl IMM (21),d4 | Getting near denormal? jle 9f | Yep, speed forward to it | Well, the bit we're looking for is somewhere in d1, (or MSB of d2), ' | and we know we can shift left at least 21 places and still not be | a denormal. movel IMM (11),d6 | Right shift count movel d1,d0 lsrl d6,d0 addl d2,d2 | Get bit from d2 addxl d1,d1 lsll d3,d1 | and shift over the rest cmpl d5,d0 | Normalized yet? jcc Ladddf$exact jra 1b 9: addl d2,d2 | Squeeze last bit from d2 addxl d1,d1 addxl d0,d0 addl d3,d4 | To denormal yet? jle Ladddf$exact 1: addl d1,d1 | Shift over another place addxl d0,d0 subl IMM (1),d4 | Done? jgt 1b jra Ladddf$exact | As the Underflow trap is enabled, we just normalize the result | without checking for the denormal exponent value. Ladddf$normal: tstl d0 | In the first word? jeq 3f | Nope, go search second word | We needn't touch d3 here. If it is non-zero, only one left shift will ' | be required. If the exponents differ by at most one, both d2 and d3 | are zero, and the result is exact. addl d2,d2 | Perhaps shift last bit out addxl d1,d1 addxl d0,d0 subl IMM (1),d4 | Adjust exponent cmpl d5,d0 jcc Lround$common | Normalized yet? | Left shift until it's normalized ' 1: addl d1,d1 | Try moving over one place addxl d0,d0 subl IMM (1),d4 | Adjust exponent cmpl d5,d0 | Normalized yet? jcs 1b | Nope, try another place Ladddf$uflock: tstl d4 | Underflow? jgt Ladddf$exact | Nope, assemble & return jlt 2f | Yep, go report cmpl d5,d0 | 2^Emin jne Ladddf$exact tstl d1 jne Ladddf$exact 2: addl IMM (1536),d4 | Adjust exponent lsll IMM (4),d4 | All we need to do now is swap d4 | assemble the double float! addl d4,d0 | Put in the exponent addl a1,d0 | and the sign movel IMM (UNDERFLOW),d7 | Show error jra $_exception_handler | and go trap | d0 is zero, the exponents of the operands are at most one different. | d1 is zero only in one special case: a is 1.0*2^n and b is the next | smallest representable number, and the difference is taken. In this | case all result bits are zero except for bit 31 of d2. In any case, | when we get here, d3 is zero! 3: subl IMM (21),d4 | Getting near denormal? | Well, the bit we're looking for is somewhere in d1, (or MSB of d2), ' | and we know we can shift left at least 21 places and still not be | a denormal. movel IMM (11),d6 | Right shift count movel d1,d0 lsrl d6,d0 movel IMM (20),d6 | Left shift count - 1 addl d2,d2 | Get bit from d2 addxl d1,d1 lsll d6,d1 | and shift over the rest cmpl d5,d0 | Normalized yet? jcc Ladddf$uflock jra 1b | a is either zero, denormalized, +/-INFINITY, or NaN Ladddf$a$spec: jeq Ladddf$a$nf | Zero or denorm? movel IMM (0),d4 | Fixup to proper exponent andl IMM (0xfffff),d0 | Strip off the false hidden bit jne Ladddf$1 | Denormalized? tstl d1 | or zero? jne Ladddf$1 cmpl d2,d3 | Is b zero, denorm, INFINITY, NaN? jcs Ladddf$b | Nope, then return it jeq Ladddf$b$nf | b zero or denorm? movel IMM (0),d3 | Fixup to proper exponent andl IMM (0xfffff),d5 | Strip off the false hidden bit jne Ladddf$b | Denormalized? tstl d6 | or zero? jne Ladddf$b cmpl a0,a1 | Compare signs jeq Ladddf$ret | and return zero if same Ladddf$ezer: lea SYM (_fpCCR),a0 | Signs differ, movel a0@(ROUND),d5 | Rounding mode determines sign subl IMM (ROUND_TO_MINUS),d5 jne Ladddf$pzer | Not round-to-minus, return +0 bset IMM (31),d0 | Return -0 jra Ladddf$pzer Ladddf$b$spec: jeq Ladddf$b$nf | Zero or denorm? movel IMM (0),d3 | Fixup to proper exponent andl IMM (0xfffff),d5 | Strip off the false hidden bit jne Ladddf$2 | Denormalized? tstl d6 | or zero? jne Ladddf$2 jra Ladddf$a | Return a as result | One of the operands is zero. Return the other, but be carefull, as | we still may need to check for Underflow! If the other operand is | non-zero and is tiny, i.e. |value| <= 2^Emin, and the Underflow trap | is enabled, the number needs to be normalized, and sent off to the | Underflow trap routine. Ladddf$b: | Return b (if a is zero) movel d5,d0 | Copy b to result movel d6,d1 movel d3,d4 movel a0,a1 | Use b's sign ' Ladddf$a: movel IMM (0),d2 | Zero the extension words movel d2,d3 jra Ladddf$unnormal | Go finish up! | a is either NaN or +/-INFINITY Ladddf$a$nf: andl IMM (0xfffff),d0 | Strip off the false hidden bit jne 1f | If not zero, must be NaN orl d0,d1 | +/-INFINITY? jeq 2f 1: btst IMM (19),d0 | Signaling NaN? jeq Ld$a$snan | Yep, go signal 2: cmpl d2,d3 | Is b +/-INFINITY, NaN? jeq Ladddf$b$nf | Yep, go check it out tstl d1 | Return jeq Ladddf$infty | +/-INFINITY jra Ld$a | or NaN | b is either NaN or +/-INFINITY, a may be also, though not | a signaling NaN Ladddf$b$nf: addl a0,a1 | Swap signs addl a1,a0 addl a0,a1 andl IMM (0xfffff),d5 | Strip off the false hidden bit jne 1f | If not zero, must be NaN orl d5,d6 | +/-INFINITY? jeq 2f 1: btst IMM (19),d5 | Signaling NaN? jeq Ld$b$snan | Yep, go signal jra Ld$b | Just pass on non-signaling NaN 2: | b is +/-INFINITY, what is a? cmpl d2,d4 | +/-INFINITY, NaN? jne Ladddf$infty | Nope, return b's +/-INFINITY ' tstl d1 | Is a +/-INFINITY? jne Ld$a | Nope, must be NaN | We know that both are +/-INFINITY. It is not valid if their signs | are different cmpl a0,a1 | Compare sign bits jne Ld$inop | Can't subtract INFINITYs! ' Ladddf$infty: jra Ld$infty | +/-INFINITY |============================================================================= | __muldf3 |============================================================================= | double __muldf3(double, double); SYM (__muldf3): link a6,IMM (-24) moveml d2-d7,sp@ pea MULTIPLY+OP1_DOUBLE+OP2_DOUBLE+DST_DOUBLE movel a6@(8),d0 | get first operand movel a6@(12),d1 | movel a6@(16),d2 | get second operand movel a6@(20),d3 | movel d0,a1 | Get a's sign bit in a1 ' bclr IMM (31),d0 | Clear sign from working reg subl d0,a1 | & isolate sign addl d2,a1 | Get sign of result bclr IMM (31),d2 | Clear sign from working reg subl d2,a1 | & isolate sign movel IMM (0xfff00000),d6 | Exponent mask & -hidden bit movel d0,d4 | Get a's exponent in d4 ' andl d6,d4 addl d6,d4 | & adjust bias subl d4,d0 | Clear exponent & set hidden bit movel d2,d5 | Get b's exponent in d5 ' andl d6,d5 addl d6,d5 | & adjust bias subl d5,d2 | Clear exponent & set hidden bit swap d4 | Right justify exponents swap d5 lsrl IMM (4),d4 lsrl IMM (4),d5 movel IMM (0x7fe),d6 | Special case literal cmpl d6,d4 | Zero, denorm, INFINITY, NaN? jcc Lmuldf$a$spec | Yep, go handle Lmuldf$1: cmpl d6,d5 | Zero, denorm, INFINITY, NaN? jcc Lmuldf$b$spec | Yep, go handle Lmuldf$2: addl d5,d4 | add exponents subl IMM (DF_RBIAS),d4 | and subtract bias | We are now ready to do the multiplication. The situation is as follows: | both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were | denormalized to start with!). We keep just 21 bits in the multiplier, | but shift them to the left to make it easier to check the next bit. | We align the multiplicand so its MSB is bit 62. The routine below shifts | the product over after each add, so this leaves the final product with | either bit 19 or 20 of the most significant word as the MSB. At that | point we must either left shift one more place, or increment the | exponent. | If one of the operands has all 32 low bits zero, we choose it to be | the multiplier. This speeds up products of many small numbers. movel d4,a0 | a0 preserves the exponent movel IMM (10),d4 movel d1,d7 | Low word of a zero? jeq 1f movel d3,d5 | Low a non-zero, how about b? jne 2f | Nope, don't switch operands ' movel d0,d6 | Low word of b is zero, movel d1,d5 | Swap operands movel d2,d0 movel d3,d7 | Zero low word of multiplier movel d1,d3 jra 4f 1: movel d3,d5 2: movel d2,d6 4: lsll d4,d0 | Align multiplier lsll d4,d6 | Align b lsll d4,d5 movel IMM (22),d2 lsrl d2,d3 orl d3,d6 | Instead of having a seperate counter, we float a flag bit through | the multiplier word. orl IMM (0x200),d0 | Set "count" bit movel IMM (0),d4 | Literal for carry propagation movel d5,d3 | First bit of multiplier is set - movel d6,d2 | we're normalized! Just do ' movel d4,d1 | copy instead of clear and add lsll IMM (2),d0 | Bump off MSB jra 2f | and go join the loop | Here we get to do the product. 1: jcc 2f | Next bit set? addl d5,d3 | Yep, add in b addxl d6,d2 | and carry to third word addxl d4,d1 | (we're not yet into the 4th) ' 2: addl d3,d3 | Shift left one place addxl d2,d2 addxl d1,d1 addl d0,d0 | Shift multiplier jne 1b | Any bits left? | We've finished the first 21 bits of the multiplier. Now on to the ' | last 32. tstl d7 | More work to do? jeq 3f | Nope, just shift words addxl d7,d7 | Shift out next bit, | and shift in flag bit! 1: jcc 2f | Next bit set? addl d5,d3 | Yep, add in b addxl d6,d2 addxl d4,d1 | and carry to third word addxl d4,d0 | and to forth word 2: addl d3,d3 | Shift left one place addxl d2,d2 addxl d1,d1 addxl d0,d0 addl d7,d7 | Shift multiplier jne 1b | Any bits left? 4: movel a0,d4 | Restore exponent btst IMM (20),d0 | Aligned? jne 1f | Yep, go bump exponent addl d3,d3 | Nope, align addxl d2,d2 addxl d1,d1 addxl d0,d0 | Now bit 20 is set jra Lround$normal | Go finish up! 3: movel d1,d0 | Low 32 bits of multiplier movel d2,d1 | are zero, we just need movel d3,d2 | to left shift 32 bits movel d4,d3 jra 4b 1: addl IMM (1),d4 | Increase exponent jra Lround$normal | Go finish up! Lmuldf$zer: movel IMM (0),d0 | Return a zero movel d0,d1 jra Ld$sign | with proper sign | a is either zero, denormalized, +/-INFINITY, or NaN Lmuldf$a$spec: jeq Lmuldf$a$big | Zero or denorm? movel IMM (0),d4 | Fixup to proper exponent andl IMM (0xfffff),d0 | Strip off the false hidden bit jne Lmuldf$a$den | Denormalized? tstl d1 | or zero? jne Lmuldf$a$den cmpl d6,d5 | Is b +/-INFINITY, NaN? jeq Lmuldf$b$big | Yep, go check it out jra Lmuldf$zer | Nope, then return zero Lmuldf$b$spec: jeq Lmuldf$b$big | Zero or denorm? movel IMM (0),d5 | Fixup to proper exponent andl IMM (0xfffff),d2 | Strip off the false hidden bit jne Lmuldf$b$den | Denormalized? tstl d3 | or zero? jne Lmuldf$b$den jra Lmuldf$zer | Return zero Lmuldf$a$big: bclr IMM (20),d0 | Remove false bit orl d0,d1 | Is a INFINITY? jeq Lmuldf$a$nf btst IMM (19),d0 | Signaling NaN? jeq Ld$a$snan | Yep, go report it cmpl d6,d5 | Is b +/-INFINITY, NaN? jne Ld$a | Nope, pass on NaN bclr IMM (20),d2 | Remove false bit orl d2,d3 | Is b INFINITY? jeq Ld$a | Yep, pass on NaN jra Lmuldf$b$nan Lmuldf$a$nf: cmpl d6,d5 | Is b zero, denorm, INFINITY, NaN? jcs Ld$infty | Nope, return +/-INFINITY jeq Lmuldf$ab$big bclr IMM (20),d2 | Remove false bit orl d2,d3 | Is b zero? jeq Ld$inop | Yep, Illegal operation! jra Ld$infty | Nope, return +/-INFINITY Lmuldf$b$big: bclr IMM (20),d2 | Remove false bit orl d2,d3 | Is b INFINITY? jne Lmuldf$b$nan orl d0,d1 | Is a zero? jne Ld$infty | Nope, return +/-INFINITY jra Ld$inop | Yep, Illegal operation! Lmuldf$ab$big: bclr IMM (20),d2 | Remove false bit orl d2,d3 | Is b INFINITY? jeq Ld$infty | Yep, return +/-INFINITY Lmuldf$b$nan: btst IMM (19),d2 | Signaling NaN? jeq Ld$b$snan | Yep, go report it jra Ld$b | No signaling NaNs, just pass it on | If a number is denormalized we put an exponent of 1 but do not put the | hidden bit back into the fraction; instead we shift left until bit 20 | (the hidden bit) is set, adjusting the exponent accordingly. We do this | to ensure that the product of the fractions is close to 1. Lmuldf$a$den: subl IMM (1),d4 | and adjust exponent addl d1,d1 | shift a left until bit 20 is set addxl d0,d0 | btst IMM (20),d0 jeq Lmuldf$a$den jra Lmuldf$1 | Lmuldf$b$den: subl IMM (1),d5 | and adjust exponent addl d3,d3 | shift b left until bit 20 is set addxl d2,d2 | btst IMM (20),d2 jeq Lmuldf$b$den jra Lmuldf$2 | |============================================================================= | __divdf3 |============================================================================= | double __divdf3(double, double); SYM (__divdf3): link a6,IMM (-24) moveml d2-d7,sp@ pea DIVIDE+OP1_DOUBLE+OP2_DOUBLE+DST_DOUBLE movel a6@(8),d2 | Get a into d2-d3 movel a6@(12),d3 movel a6@(16),d6 | and b into d6-d5 movel a6@(20),d5 movel d2,a1 | Get a's sign bit in a1 ' bclr IMM (31),d2 | Clear sign from working reg subl d2,a1 addl d6,a1 | Get sign of result bclr IMM (31),d6 | Clear sign from working reg subl d6,a1 | Bit 31 is sign of result movel IMM (0xfff00000),d0 | Exponent mask & -hidden bit movel d2,d4 | Get a's exponent in d4 ' andl d0,d4 addl d0,d4 | & adjust bias subl d4,d2 | Clear exponent & set hidden bit movel d6,d1 | Get b's exponent in d1 ' andl d0,d1 addl d0,d1 | & adjust bias subl d1,d6 | Clear exponent & set hidden bit swap d4 | Right justify exponents swap d1 lsrl IMM (4),d4 lsrl IMM (4),d1 movel IMM (0x7fe),d0 | Special case literal cmpl d0,d4 | Zero, denorm, INFINITY, NaN? jcc Ldivdf$a$spec | Yep, go handle Ldivdf$1: cmpl d0,d1 | Zero, denorm, INFINITY, NaN? jcc Ldivdf$b$spec | Yep, go handle Ldivdf$2: subl d1,d4 | Subtract divisor exponent addl IMM (DF_RBIAS),d4 | and add back bias | We add -divisor instead of subtracting as that way we can use the | carry bit directly to build up the quotient. The initial dividend, | and final remainder, is in d2,d3. The quotient is built up in d1. | When the first 21 bits are complete, they are moved to d0, and the | next 32 bits are loaded into d1. Then a 54th bit is created in bit | 31 of d2. | | Instead of using a seperate counter, we float a flag bit through | d1, and loop until it pops into the C flag. We zero d0 initially, | and use that to flag when we've got all 53 bits of the quotient. ' movel d5,a0 | Save plus low word negl d5 | Form minus low word movel d6,d7 negxl d7 | Form minus high word movel IMM (0),d0 | Still need low 32 quotient bits movel IMM (0x800),d1 | Set counter flag addl d5,d3 | Quotient < 1? addxl d7,d2 jcs 5f | Nope, enter subtract loop subl IMM (1),d4 | Yep, adjust exponent | and enter add back loop | Remainder is minus: Add back divisor 1: addl d3,d3 addxl d2,d2 addl a0,d3 addxl d6,d2 jcs 5f 2: addxl d1,d1 jcc 1b tstl d0 | Done 53 bits yet? jne 3f movel d1,d0 | Nope, 32 more to go movel IMM (1), d1 jra 1b 3: addl d3,d3 addxl d2,d2 addl a0,d3 | Get 54th bit for rounding addxl d6,d2 jra 7f | Remainder is plus: Subtract divisor 4: addl d3,d3 addxl d2,d2 addl d5,d3 addxl d7,d2 jcc 2b 5: addxl d1,d1 jcc 4b tstl d0 | Done 53 bits yet? jne 6f movel d1,d0 | Nope, 32 more to go movel IMM (1), d1 jra 4b 6: addl d3,d3 addxl d2,d2 addl d5,d3 | Get 54th bit for rounding addxl d7,d2 7: jcs 1f | 54th bit set? addl a0,d3 | Nope, addxl d6,d2 | correct the remainder jra 2f 1: bset IMM (31),d2 | Set 54th bit 2: jra Lround$normal | Go finish up! | a is either zero, denormalized, +/-INFINITY, or NaN Ldivdf$a$spec: jeq Ldivdf$a$big | Zero or denorm? movel IMM (0),d4 | Fixup to proper exponent andl IMM (0xfffff),d2 | Strip off the false hidden bit jne Ldivdf$a$den | Denormalized? tstl d3 | or zero? jne Ldivdf$a$den cmpl d0,d1 | Is b +/-INFINITY, NaN? jeq Ldivdf$b$big | Yep, go check it out jcs Ldivdf$zer | Normal number, return zero bclr IMM (20),d6 | Strip off the false hidden bit orl d6,d5 | Is b zero? jeq Ld$inop | Yep, illegal operation: 0/0 Ldivdf$zer: movel IMM (0),d0 | Return a zero movel d0,d1 jra Ld$sign | with proper sign Ldivdf$b$spec: jeq Ldivdf$b$big | Zero or denorm? movel IMM (0),d1 | Fixup to proper exponent andl IMM (0xfffff),d6 | Strip off the false hidden bit jne Ldivdf$b$dck | Denormalized? tstl d5 | or zero? jne Ldivdf$b$dck orl d2,d3 | a zero too? jeq Ld$inop | 0/0, go report Invalid operation movel IMM (DF_INFINITY),d0 | Return a signed INFINITY addl a1,d0 movel IMM (0),d1 movel IMM (DIVIDE_BY_ZERO),d7 jra $_exception_handler Ldivdf$a$big: bclr IMM (20),d2 | Remove false bit orl d2,d3 | Is a INFINITY? jeq Ldivdf$a$nf btst IMM (19),d2 | Signaling NaN? jeq Ld$a$snan | Yep, go report it cmpl d0,d1 | Is b zero, denorm, INFINITY, NaN? jne Ld$a | Nope, pass on NaN unchanged bclr IMM (20),d6 | Remove false bit orl d6,d5 | Is b INFINITY? jeq Ld$a | Yep, pass on NaN jra Ldivdf$b$nan Ldivdf$a$nf: cmpl d0,d1 | Is b INFINITY, NaN? jne Ld$infty | Nope, return +/-INFINITY bclr IMM (20),d6 | Remove false bit orl d6,d5 | Is b INFINITY? jeq Ld$inop | Yep, INFINITY/INFINITY is invalid jra Ldivdf$b$nan Ldivdf$b$big: bclr IMM (20),d6 | Remove false bit orl d6,d5 | Is b INFINITY? jeq Ldivdf$zer | Yep, return zero Ldivdf$b$nan: btst IMM (19),d6 | Signaling NaN? jeq Ld$b$snan | Yep, go report it jra Ld$b | No signaling NaNs, just pass it on | If a number is denormalized we put an exponent of 1 but do not put the | hidden bit back into the fraction; instead we shift left until bit 21 | (the hidden bit) is set, adjusting the exponent accordingly. We do this | to ensure that the product of the fractions is close to 1. Ldivdf$a$den: subl IMM (1),d4 | and adjust exponent addl d3,d3 | shift a left until bit 20 is set addxl d2,d2 | btst IMM (20),d2 jeq Ldivdf$a$den jra Ldivdf$1 | Ldivdf$b$dck: movel d2,d7 orl d3,d7 | Is a zero? jeq Ldivdf$zer | Yep, return zero Ldivdf$b$den: subql IMM (1),d1 | and adjust exponent addl d5,d5 | shift b left until bit 20 is set addxl d6,d6 | btst IMM (20),d6 | jeq Ldivdf$b$den jra Ldivdf$2 | |============================================================================= | __negdf2 |============================================================================= | double __negdf2(double); SYM (__negdf2): movel sp@(4),d0 | get number to negate in d0-d1 bchg IMM (31),d0 | negate movel d0,d1 | make a positive copy (for the tests) bclr IMM (31),d1 cmpl IMM (DF_INFINITY),d1 | compare to +INFINITY jge 2f | NaN or +/-INFINITY check further 1: movel sp@(8),d1 | Get lower word rts 2: jne 3f | NaN? movel sp@(8),d1 | Get lower word jeq 1b | Return +/-INFINITY 3: btst IMM (19),d0 | Signaling NaN? jne 1b | Nope, just pass it on link a6,IMM (-24) moveml d2-d7,sp@ pea NEGATE+OP1_DOUBLE+DST_DOUBLE jra Ld$a$snan |============================================================================= | __cmpdf2 |============================================================================= .globl $_cmpdf GREATER = 1 LESS = -1 EQUAL = 0 | int __cmpdf2(double, double); SYM (__cmpdf2): lea GREATER,a1 | Unordered return value jra $_cmpdf | Go perform comparison $_cmpdf: link a6,IMM (-24) moveml d2-d7,sp@ pea COMPARE+OP1_DOUBLE+OP2_DOUBLE+DST_INT movel a6@(8),d0 | get first operand movel a6@(12),d1 | movel a6@(16),d2 | get second operand movel a6@(20),d3 | | First check if a and/or b are (+/-) zero and in that case clear | the sign bit. movel d0,d6 | copy signs into d6 (a) and d5(b) bclr IMM (31),d0 | and clear signs in d0 and d2 movel d2,d5 | bclr IMM (31),d2 | cmpl IMM (DF_INFINITY),d0 | check for a == NaN jcc Lcmpdf$a$big | if d0 >= 0x7ff00000, a is NaN or INFINITY movel d0,d4 | copy into d4 to test for zero orl d1,d4 | jeq Lcmpdf$a$0 | Lcmpdf$0: cmpl IMM (DF_INFINITY),d2 | check for b == NaN jcc Lcmpdf$b$big | if d2 >= 0x7ff00000, b is NaN or INFINITY movel d2,d4 | orl d3,d4 | jeq Lcmpdf$b$0 | Lcmpdf$1: | Check the signs eorl d6,d5 jpl 1f | If the signs are not equal check if a >= 0 tstl d6 jpl Lcmpdf$a$gt$b | if (a >= 0 && b < 0) => a > b jra Lcmpdf$b$gt$a | if (a < 0 && b >= 0) => a < b 1: | If the signs are equal check for < 0 tstl d6 jpl 1f | If both are negative exchange them movel d0,d5 movel d2,d0 movel d5,d2 movel d1,d5 movel d3,d1 movel d5,d3 1: | Now that they are positive we just compare them as longs (does this also | work for denormalized numbers?)(Yep). cmpl d0,d2 jhi Lcmpdf$b$gt$a | |b| > |a| jne Lcmpdf$a$gt$b | |b| < |a| | If we got here d0 == d2, so we compare d1 and d3. cmpl d1,d3 jhi Lcmpdf$b$gt$a | |b| > |a| jne Lcmpdf$a$gt$b | |b| < |a| | If we got here a == b. movel IMM (EQUAL),d0 jra Lcmpdf$ret Lcmpdf$a$gt$b: movel IMM (GREATER),d0 jra Lcmpdf$ret Lcmpdf$b$gt$a: movel IMM (LESS),d0 Lcmpdf$ret: moveml sp@(4),d2-d7 | XXX if frame pointer is ever removed, stack pointer must | be adjusted here. unlk a6 rts Lcmpdf$a$0: movel d0,d6 jra Lcmpdf$0 Lcmpdf$b$0: movel d2,d5 jra Lcmpdf$1 Lcmpdf$a$big: jne Lcmpdf$a$nan | if not equal it's an NaN ' tstl d1 jeq Lcmpdf$0 Lcmpdf$a$nan: btst IMM (19),d0 | Is a a signaling NaN? jeq Lcmpdf$inv | Yep, go signal cmpl IMM (DF_INFINITY),d2 | Is b a signaling NaN? jhi 1f | Maybe jne Lcmpdf$uno | Nope tstl d3 | INFINITY? jeq Lcmpdf$uno | Yep 1: btst IMM (19),d2 | Signaling NaN? jeq Lcmpdf$inv | Yep, go signal jra Lcmpdf$uno Lcmpdf$b$big: jne Lcmpdf$b$nan | if not equal it's an NaN ' tstl d3 jeq Lcmpdf$1 Lcmpdf$b$nan: btst IMM (19),d2 | Is b a signaling NaN? jeq Lcmpdf$inv | Yep, report it Lcmpdf$uno: movel a1,d0 | Report no ordering for NaNs btst IMM (0),d0 | Error on unordered? jne Lcmpdf$ret Lcmpdf$inv: movel a1,d0 | Get unordered code asrl IMM (1),d0 | and dump invalid bit orl IMM (1),d0 jra Ld$inv | and report exception |============================================================================= | Rounding and Exception Handling |============================================================================= | The non-zero value to be rounded and returned is in registers d0-d1-d2-d3, | with the exponent in register d4, the sign in a1, and the operation type on | the top of the stack. The radix point is between bits 19 and 20 of d0, and | the exponent biased by 1022 (adjusted before returning to 1023). | | Three exceptions are checked for: Overflow, Underflow, and Inexact. | | Overflow is checked after rounding. If detected, and trapped, the trap | handler is passed the rounded value, with the exponent adjusted by a bias, | 1536 in the case of DOUBLE. The Inexact flag is set if the rounded value | is not exact. If not trapped, the value is set properly for the rounding | mode, and the Overflow and Inexact flags are set. The Inexact trap is | issued if enabled in this case. | | If the |value| <= 2^Emin (Emin = -1022 for DOUBLE) and the Underflow trap | is enabled, the value is normalize (if it isn't already), rounded, and sent ' | to the trap handler with the exponent bias adjusted. If the trap is not | enabled, the possibly denormalized value is rounded and, if inexact, the | Underflow flag is set. The Inexact flag is set if the result is not exact | in either case, and the Inexact trap taken if enabled, and the Underflow | trap is not. | We've got a normalized number, but need to denormalize it. ' denormal: lea SYM (_fpCCR),a0 btst IMM (1),a0@(TRAPE+3) | Check for Underflow trap enabled jne Lround$common | Yep, remain normalized cmpl IMM (-DF_SIGNIF_DIG-1),d4 | Really tiny? jls degen | Yep, just summarize movel IMM (32),d6 addl d6,d4 | Get left shift count jle 1f | Moving a whole word? subl d4,d6 | Get right shift count orl d2,d3 | Summarize low order bits movel d1,d2 lsll d4,d2 lsrl d6,d1 movel d0,d7 lsll d4,d0 orl d0,d1 movel d7,d0 lsrl d6,d0 movel IMM (0),d4 jra Lround$common 1: jeq 1f | Shift even 32-bit words? addl d6,d4 | Get left shift count subl d4,d6 | and right shift count orl d2,d3 | Summarize lost bits orl d1,d3 | in d3 movel d0,d2 | Get least significant lsll d4,d2 | bits into d2 movel d0,d1 | and the rest into d1 lsrl d6,d1 2: movel IMM (0),d0 | Show zero in d0 movel d0,d4 | and denorm exponent jra Lround$common 1: orl d2,d3 | Summarize in d3 movel d1,d2 | and shift others movel d0,d1 | down one place jra 2b degen: movel IMM (1),d2 | Show we're non-zero ' movel IMM (0),d1 | But all significance is gone jra 2b Lround$normal: | We first check for underflow. tstl d4 | Are we tiny? jlt denormal | Maybe, go check Lround$common: lea SYM (_fpCCR),a0 movel d3,d7 | Anything down here? orl d2,d7 jeq Lround$done | Exact! We're done! ' movel IMM (INEXACT_RESULT),d7 | Set Inexact flag movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lround$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Lround$to$minus jlt Lround$to$zero | jra Lround$to$plus Lround$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lround$done | Yep, truncate jra Lround$inc | Nope, increment Lround$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lround$done | Nope, truncate jra Lround$inc | Yep, increment Lround$to$nearest: | Now round: we do it as follows: after the shifting we can write the | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2. | If delta < 1, do nothing. If delta > 1, add 1 to f. | If delta == 1, we make sure the rounded number will be even (odd?) | (after shifting). addl d3,d3 addxl d2,d2 | Is delta < 1? jcc Lround$done | Yep, do nothing jne Lround$inc | Is delta == 1? movel IMM (1),d3 | If so round to even andl d1,d3 | bit 1 is the last significant bit addl d3,d1 jra 1f Lround$inc: addl IMM (1),d1 | Else add 1 1: movel IMM (0),d2 addxl d2,d0 | Rounding may have changed our exponent. If so, adjust. btst IMM (21),d0 jeq Lround$done lsrl IMM (1),d0 | This only occurs when d1 == 0! addl IMM (1),d4 | Fixup exponent Lround$to$zero: Lround$done: | Now we're properly normalized and rounded, or, if the Underflow trap ' | is not enabled, maybe denormalized and rounded. If we're tiny ' | (|value| <= 2^Emin) the Tiny flag (proposed Underflow flag) is set. | If the rounded value is not exact, the Inexact flag is set. movel a0@(TRAPE),d5 | Get trap enables cmpl IMM (DF_MAX_EXP-1),d4 | Overflow? jge Ltrap$overflow | Nope, go check for Underflow tstl d4 | Check for Underflow jle Ltrap$underflow Ltrap$check: | Check for trapped exceptions. First we put the number back together. swap d4 | Put the exponent in place lsll IMM (4),d4 addl d4,d0 | and add one to exponent tstl d7 | Any exceptions? jne Ltrap$assemble | Yep, go handle addl a1,d0 | Nope, add in sign moveml sp@(4),d2-d7 | We're done ' | XXX if frame pointer is ever removed, stack pointer must | be adjusted here. unlk a6 | and return rts Ltrap$underflow: jlt 2f | Really small (& trap enabled)? cmpl IMM (0x100000),d0 | |value| <= 2^Emin? jgt Ltrap$check | Nope, return jlt 1f | Yep, set Underflow tstl d1 | |value| == 2^Emin? jne Ltrap$check | Nope 1: btst IMM (1),d5 | Underflow trap enabled? jne 2f | Yep, adjust the exponent btst IMM (0),d7 | Result Inexact? jeq Ltrap$check | Nope, won't set Underflow ' jra 3f 2: addl IMM (1536), d4 | Adjust exponent 3: addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Ltrap$check Ltrap$overflow: addl IMM (OVERFLOW),d7 | Set the Overflow flag btst IMM (2),d5 | Overflow trap enabled? jeq 1f | Nope, go adjust value subl IMM (1536),d4 | Adjust exponent jra Ltrap$check | and go trap 1: movel IMM (INEXACT_RESULT+OVERFLOW),d7 | Inexact too movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Loflow$nearest | Yep, return infinity subl IMM (ROUND_TO_PLUS),d6 jhi Loflow$minus jlt Loflow$zero | jra Loflow$plus Loflow$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Loflow$nearest | Nope, go to infinity jra Loflow$zero | Yep, just go big Loflow$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Loflow$nearest | Yep, go to infinity jra Loflow$zero | Nope, just go big Loflow$nearest: movel IMM (0),d0 | Build an infinity value movel d0,d1 movel IMM (DF_MAX_EXP),d4 jra Ltrap$check Loflow$zero: movel IMM (0x1fffff),d0 | Build max value movel IMM (-1),d1 movel IMM (DF_MAX_EXP-2),d4 jra Ltrap$check Ltrap$assemble: addl a1,d0 | Put the sign back jra $_exception_handler Ld$a$snan: movel a6@(8),d0 | Return a movel a6@(12),d1 jra Ld$snan Ld$b$snan: movel a6@(16),d0 | Return b movel a6@(20),d1 Ld$snan: bset IMM (19),d0 | Change to Quiet NaN jra Ld$inv Ld$inop: | Invalid Operation. Return a quiet NaN and set the exception flags movel IMM (Q_NaN),d0 movel IMM (Q_NaN),d1 Ld$inv: movel IMM (INVALID_OPERATION),d7 jra $_exception_handler Ld$infty: | Return a properly signed INFINITY and set the exception flags movel IMM (DF_INFINITY),d0 movel IMM (0),d1 Ld$sign: addl a1,d0 | Sign bit Ld$return: moveml sp@(4),d2-d7 | XXX if frame pointer is ever removed, stack pointer must | be adjusted here. unlk a6 | and return rts Ld$a: movel a6@(8),d0 | Return a movel a6@(12),d1 jra Ld$return Ld$b: movel a6@(16),d0 | Return b movel a6@(20),d1 jra Ld$return #endif /* L_double */ #ifdef L_float | Entry points: .globl SYM (__addsf3) .globl SYM (__subsf3) .globl SYM (__mulsf3) .globl SYM (__divsf3) .globl SYM (__negsf2) .globl SYM (__cmpsf2) .text .even |============================================================================= |============================================================================= | single precision routines |============================================================================= |============================================================================= | A single precision floating point number (float) has the format: | | struct _float { | unsigned int sign : 1; /* sign bit */ | unsigned int exponent : 8; /* exponent, shifted by 126 */ | unsigned int fraction : 23; /* fraction */ | } float; | | Thus sizeof(float) = 4 (32 bits). | | All the routines are callable from C programs, and return the result | in the single register d0. They also preserve all registers except | d0-d1 and a0-a1. |============================================================================= | __subsf3 |============================================================================= | float __subsf3(float, float); SYM (__subsf3): link a6,IMM (-24) moveml d2-d7,sp@ pea SUB+OP1_FLOAT+OP2_FLOAT+DST_FLOAT movel a6@(12),d6 | get second operand bchg IMM (31),d6 | change sign of second operand jra 1f | and fall through |============================================================================= | __addsf3 |============================================================================= | float __addsf3(float, float); SYM (__addsf3): link a6,IMM (-24) moveml d2-d7,sp@ pea ADD+OP1_FLOAT+OP2_FLOAT+DST_FLOAT movel a6@(12),d6 | get second operand 1: movel a6@(8),d0 | get first operand movel d0,a1 | Get a's sign bit in a1 ' movel d6,a0 | Get b's sign in a0 ' bclr IMM (31),d0 | Clear signs from working regs bclr IMM (31),d6 subl d0,a1 subl d6,a0 movel IMM (0xff800000),d2 | Exponent mask & -hidden bit movel d0,d4 | Get a's exponent in d4 ' andl d2,d4 addl d2,d4 | & adjust bias subl d4,d0 | Clear exponent & set hidden bit movel d6,d3 | Get b's exponent in d3 ' andl d2,d3 addl d2,d3 | & adjust bias subl d3,d6 | Clear exponent & set hidden bit swap d4 | Right justify exponents swap d3 lsrl IMM (7),d4 lsrl IMM (7),d3 movel IMM (0xfe),d2 | Special case literal cmpl d2,d4 | Zero, denorm, INFINITY, NaN? jcc Laddsf$a$spec | Yep, go handle Laddsf$1: cmpl d2,d3 | Zero, denorm, INFINITY, NaN? jcc Laddsf$b$spec | Yep, go handle Laddsf$2: | Both operands are finite and non-zero. Here we shift the numbers until | their exponents are the same. subl d4,d3 | compare the exponents jeq Laddsf$3 | If equal we don't need to shift ' jmi 9f | Branch if second exponent is lower addl d3,d4 | Put largest exponent into d4 addl a0,a1 | Swap a & b addl a1,a0 addl a0,a1 movel d0,d2 | Swap fractions movel d6,d0 movel d2,d6 jra 1f | Here we have a's exponent larger than b's, so we have to shift b. 9: negl d3 | Get positive right-shift count | If difference is too large we don't shift, just summarize b's | non-zeroness in d1 1: movel IMM (32), d5 cmpl d5,d3 jge Laddsf$b$small subl d3,d5 | Get left shift count movel d6,d1 lsll d5,d1 lsrl d3,d6 jra Laddsf$4 | If one of the numbers is pretty small (difference of exponents >= | DF_SIGNIF_DIG+1) we just summarize it to show an Inexact operation. Laddsf$b$small: movel IMM (2),d1 | Show b is piddlingly small movel IMM (0),d6 jra Laddsf$4 | Now go do operation & rounding Laddsf$3: movel IMM (0),d1 | Equal exponents; zap extension Laddsf$4: | Now we have the numbers in d0 and d6/d1, the exponent in d4, and | the signs in a1/a0. | Here we have to decide whether to add or subtract the numbers: cmpl a0,a1 | Compare the signs jne Lsubsf$0 | If the signs are different we get | to subtract addl d6,d0 | else we add the numbers | Before rounding normalize so bit #SF_SIGNIF_DIG-1 is set (we consider | the case of denormalized numbers in the rounding routine itself). | As in the addition (not in the subtraction!) we could have set | one more bit we check this: btst IMM (SF_SIGNIF_DIG),d0 jeq Laddsf$unnormal movel IMM (31),d5 lsrl IMM (1),d1 movel d0,d6 lsll d5,d6 orl d6,d1 lsrl IMM (1),d0 addl IMM (1),d4 jra Lround$common | Go finish up! Lsubsf$0: | Here we do the subtraction. negl d1 subxl d6,d0 | jeq Laddsf$ezer | if zero just exit jpl Laddsf$unnormal | if positive skip the following movel a0,a1 | Use b's sign ' negl d0 | We don't need to touch d1 as ' | the exponents are equal. d1 | is therfore zero | Normalization for add/subtract | When the exponents of the operands differ by more than one, there will | at most be one left shift. This means that whenever there is more than | one left shift, the result is exact, and rounding may be skipped. Laddsf$unnormal: movel IMM (0x800000),d6 | MSB position cmpl d6,d0 | Already normalized? jcc Lround$common lea SYM (_fpCCR),a0 btst IMM (1),a0@(TRAPE+3) | Check for Underflow trap enabled jne Laddsf$normal | Yep, don't check for denorm ' tstl d4 | Already denormalized? jle Laddsf$exact | Yep, don't shift ' | First we find the MSB of the value. If the exponents differ by at most one, | no more than the MSB of d1 is set. addl d1,d1 | Perhaps shift last bit out addxl d0,d0 subl IMM (1),d4 | Denormal yet? jeq Laddsf$exact | Yep, but at least it's exact! ' cmpl d6,d0 jcc Lround$common | or perhaps not? | Left shift until it's either normalized, or we've reached the denormal | exponent value 1: addl d0,d0 | Try moving over one place subl IMM (1),d4 | Denormal yet? jeq Laddsf$exact | Yep, go adjust other words cmpl d6,d0 | Normalized yet? jcs 1b | Nope, try another place Laddsf$exact: lsll IMM (7),d4 | All we need to do now is swap d4 | assemble the double float! addl d4,d0 | Put in the exponent Laddsf$ret: addl a1,d0 | and the sign Laddsf$pzer: moveml sp@(4),d2-d7 | We're done! ' | XXX if frame pointer is ever removed, stack pointer must | be adjusted here. unlk a6 rts | As the Underflow trap is enabled, we just normalize the result | without checking for the denormal exponent value. | When the exponents of the operands differ by more than one, there will | at most be one left shift. This means that whenever there is more than | one left shift, the result is exact, and rounding may be skipped. Laddsf$normal: addl d1,d1 | Perhaps shift last bit out addxl d0,d0 subl IMM (1),d4 | Adjust exponent cmpl d6,d0 jcc Lround$common | Normalized yet? | Left shift until it's normalized ' 1: addl d0,d0 | Try moving over one place subl IMM (1),d4 | Adjust exponent cmpl d6,d0 | Normalized yet? jcs 1b | Nope, try another place tstl d4 | Underflow? jgt Laddsf$exact | Nope, assemble & return jlt 2f | Yep, go report cmpl d6,d0 | 2^Emin jne Laddsf$exact 2: addl IMM (192),d4 | Adjust exponent lsll IMM (7),d4 | All we need to do now is swap d4 | assemble the double float! addl d4,d0 | Put in the exponent addl a1,d0 | and the sign movel IMM (UNDERFLOW),d7 | Show error jra $_exception_handler | and go trap | a is either zero, denormalized, +/-INFINITY, or NaN Laddsf$a$spec: jeq Laddsf$a$nf | Zero or denorm? movel IMM (0),d4 | Fixup to proper exponent andl IMM (0x7fffff),d0 | Strip off the false hidden bit jne Laddsf$1 | Denormalized? cmpl d2,d3 | Is b zero, denorm, INFINITY, NaN? jcs Laddsf$b | Nope, then return it jeq Laddsf$b$nf | b zero or denorm? movel IMM (0),d3 | Fixup to proper exponent andl IMM (0x7fffff),d6 | Strip off the false hidden bit jne Laddsf$b | Denormalized? cmpl a0,a1 | Compare signs jeq Laddsf$ret | and return zero if same Laddsf$ezer: lea SYM (_fpCCR),a0 | Signs differ, movel a0@(ROUND),d6 | Rounding mode determines sign subl IMM (ROUND_TO_MINUS),d6 jne Laddsf$pzer | Not round-to-minus, return +0 bset IMM (31),d0 | Return -0 jra Laddsf$pzer Laddsf$b$spec: jeq Laddsf$b$nf | Zero or denorm? movel IMM (0),d3 | Fixup to proper exponent andl IMM (0x7fffff),d6 | Strip off the false hidden bit jne Laddsf$2 | Denormalized? jra Laddsf$a | Return a as result | One of the operands is zero. Return the other, but be carefull, as | we still may need to check for Underflow! If the other operand is | non-zero and is tiny, i.e. |value| <= 2^Emin, and the Underflow trap | is enabled, the number needs to be normalized, and sent off to the | Underflow trap routine. Laddsf$b: | Return b (if a is zero) movel d6,d0 | Copy b to result movel d3,d4 movel a0,a1 | Use b's sign ' Laddsf$a: movel IMM (0),d1 | Zero the extension word jra Laddsf$unnormal | Go finish up! | a is either NaN or +/-INFINITY Laddsf$a$nf: andl IMM (0x7fffff),d0 | Strip off the false hidden bit jeq 1f | +/-INFINITY? btst IMM (22),d0 | Signaling NaN? jeq Lf$a$snan | Yep, go signal 1: cmpl d2,d3 | Is b +/-INFINITY, NaN? jeq Laddsf$b$nf | Yep, go check it out tstl d0 | Return jeq Laddsf$infty | +/-INFINITY jra Lf$a | or NaN | b is either NaN or +/-INFINITY, a may be also, though not | a signaling NaN Laddsf$b$nf: addl a0,a1 | Use b's sign ' addl a1,a0 addl a0,a1 andl IMM (0x7fffff),d6 | Strip off the false hidden bit jeq 1f | +/-INFINITY? btst IMM (22),d6 | Signaling NaN? jeq Lf$b$snan | Yep, go signal jra Lf$b | Just pass on non-signaling NaN | b is +/-INFINITY, what is a? 1: cmpl d2,d4 | +/-INFINITY, NaN? jne Laddsf$infty | Nope, return b's +/-INFINITY ' tstl d0 | Is a +/-INFINITY? jne Lf$a | Nope, must be NaN | We know that both are +/-INFINITY. It is not valid if their signs | are different cmpl a0,a1 | Compare sign bits jne Lf$inop | Can't subtract INFINITYs! ' Laddsf$infty: jra Lf$infty | +/-INFINITY |============================================================================= | __mulsf3 |============================================================================= | float __mulsf3(float, float); SYM (__mulsf3): link a6,IMM (-24) moveml d2-d7,sp@ pea MULTIPLY+OP1_FLOAT+OP2_FLOAT+DST_FLOAT movel a6@(8),d0 | get a into d0 movel a6@(12),d1 | and b into d1 movel d0,a1 | Get a's sign bit in a1 ' bclr IMM (31),d0 | Clear sign from working reg subl d0,a1 addl d1,a1 | Get sign of result bclr IMM (31),d1 | Clear sign from working reg subl d1,a1 | Bit 31 is sign of result movel IMM (0xff800000),d6 | Exponent mask & -hidden bit movel d0,d4 | Get a's exponent in d4 ' andl d6,d4 addl d6,d4 | & adjust bias subl d4,d0 | Clear exponent & set hidden bit movel d1,d5 | Get b's exponent in d5 ' andl d6,d5 addl d6,d5 | & adjust bias subl d5,d1 | Clear exponent & set hidden bit swap d4 | Right justify exponents swap d5 lsrl IMM (7),d4 lsrl IMM (7),d5 movel IMM (0xfe),d6 | Special case literal cmpl d6,d4 | Zero, denorm, INFINITY, NaN? jcc Lmulsf$a$spec | Yep, go handle Lmulsf$1: cmpl d6,d5 | Zero, denorm, INFINITY, NaN? jcc Lmulsf$b$spec | Yep, go handle Lmulsf$2: addl d5,d4 | add exponents subl IMM (SF_RBIAS),d4 | and subtract bias | We are now ready to do the multiplication. The situation is as follows: | both a and b have bit 24 set (even if they were denormalized to start with!). We keep just 21 bits in the multiplier, | but shift them to the left to make it easier to check the next bit. | We align a so its MSB is bit 31. This leaves the final product with | either bit 23 or 24 of the most significant word as the MSB. At that | point we must either left shift one more place, or increment the | exponent. lsll IMM (8),d0 | Align a movel d1, d2 muluw d0, d1 | al * bl Low bits movel d2, d6 swap d6 movel d0, d3 muluw d6, d3 | bu * al swap d0 muluw d0, d2 | au * bl swap d3 muluw d6, d0 | bu * au high bits swap d2 movel d3, d6 clrw d6 eorl d6, d3 addl d6, d1 | Add (bu*al) to low bits addxl d3, d0 | and to high bits movew d2, d3 eorl d3, d2 addl d2, d1 | Add (au*bl) to low bits addxl d3, d0 | and to high bits btst IMM (SF_SIGNIF_DIG-1),d0 | Aligned? jne 1f | Yep, go bump exponent addl d1,d1 | Nope, align addxl d0,d0 | Now bit 23 is set jra Lround$normal | Go finish up! 1: addl IMM (1),d4 | Increase exponent jra Lround$normal | Go finish up! Lmulsf$zer: movel IMM (0),d0 | Return a zero jra Lf$sign | with proper sign | a is either zero, denormalized, +/-INFINITY, or NaN Lmulsf$a$spec: jeq Lmulsf$a$big | Zero or denorm? movel IMM (0),d4 | Fixup to proper exponent andl IMM (0x7fffff),d0 | Strip off the false hidden bit jne Lmulsf$a$den | Denormalized? cmpl d6,d5 | Is b +/-INFINITY, NaN? jne Lmulsf$zer | Nope, then return zero andl IMM (0x7fffff),d1 | Remove false bit jeq Lf$inop | INFINITY * 0 is illegal jra Lmulsf$b$nan Lmulsf$b$spec: jeq Lmulsf$b$big | Zero or denorm? movel IMM (0),d5 | Fixup to proper exponent andl IMM (0x7fffff),d1 | Strip off the false hidden bit jne Lmulsf$b$den | Denormalized? jra Lmulsf$zer | Return zero Lmulsf$a$big: andl IMM (0x7fffff),d0 | Remove false bit jeq Lmulsf$a$nf | Is a INFINITY? btst IMM (22),d0 | Signaling NaN? jeq Lf$a$snan | Yep, go report it cmpl d6,d5 | Is b zero, denorm, INFINITY, NaN? jne Lf$a | Nope, pass on NaN andl IMM (0x7fffff),d1 | Remove false bit jeq Lf$a | Is b INFINITY? jra Lmulsf$b$nan Lmulsf$a$nf: cmpl d6,d5 | Is b zero, denorm, INFINITY, NaN? jcs Lf$infty | Nope, return +/-INFINITY jeq Lmulsf$b$big andl IMM (0x7fffff),d1 | Remove false bit jeq Lf$inop | Is b zero? jra Lf$infty | Nope, return +/-INFINITY Lmulsf$b$big: andl IMM (0x7fffff),d1 | Remove false bit jeq Lf$infty | Is b INFINITY? Lmulsf$b$nan: btst IMM (22),d1 | Signaling NaN? jeq Lf$b$snan | Yep, go report it jra Lf$b | No signaling NaNs, just pass it on | If a number is denormalized we put an exponent of 1 but do not put the | hidden bit back into the fraction; instead we shift left until bit 20 | (the hidden bit) is set, adjusting the exponent accordingly. We do this | to ensure that the product of the fractions is close to 1. Lmulsf$a$den: subl IMM (1),d4 | Adjust exponent addl d0,d0 | shift a left until bit 23 is set btst IMM (SF_SIGNIF_DIG-1),d0 jeq Lmulsf$a$den jra Lmulsf$1 Lmulsf$b$den: subl IMM (1),d5 | Adjust exponent addl d1,d1 | shift b left until bit 23 is set btst IMM (SF_SIGNIF_DIG-1),d1 jeq Lmulsf$b$den jra Lmulsf$2 |============================================================================= | __divsf3 |============================================================================= | float __divsf3(float, float); SYM (__divsf3): link a6,IMM (-24) moveml d2-d7,sp@ pea DIVIDE+OP1_FLOAT+OP2_FLOAT+DST_FLOAT movel a6@(8),d1 | Get a into d1 movel a6@(12),d6 | and b into d6 movel d1,a1 | Get a's sign bit in a1 ' bclr IMM (31),d1 | Clear sign from working reg subl d1,a1 addl d6,a1 | Get sign of result bclr IMM (31),d6 | Clear sign from working reg subl d6,a1 | Bit 31 is sign of result movel IMM (0xff800000),d0 | Exponent mask & -hidden bit movel d1,d4 | Get a's exponent in d4 ' andl d0,d4 addl d0,d4 | & adjust bias subl d4,d1 | Clear exponent & set hidden bit movel d6,d2 | Get b's exponent in d2 ' andl d0,d2 addl d0,d2 | & adjust bias subl d2,d6 | Clear exponent & set hidden bit swap d4 | Right justify exponents swap d2 lsrl IMM (7),d4 lsrl IMM (7),d2 movel IMM (SF_MAX_EXP-1),d0 | Special case literal cmpl d0,d4 | Zero, denorm, INFINITY, NaN? jcc Ldivsf$a$spec | Yep, go handle Ldivsf$1: cmpl d0,d2 | Zero, denorm, INFINITY, NaN? jcc Ldivsf$b$spec | Yep, go handle Ldivsf$2: subl d2,d4 | Subtract divisor exponent addl IMM (SF_RBIAS),d4 | and add back bias | We add -divisor instead of subtracting as that way we can use the | carry bit directly to build up the quotient. The initial dividend, | and final remainder, is in d1. The quotient is built up in d0. | Then a 25th bit is created in bit 31 of d1. | | Instead of using a seperate counter, we float a flag bit through | d0, and loop until it pops into the C flag. movel d6,d5 | Save divisor negl d5 | Form minus divisor movel IMM (0x100),d0 | Set counter flag addl d5,d1 | Quotient < 1? jcs 5f | Nope, enter subtract loop subl IMM (1),d4 | Yep, adjust exponent | and enter add back loop | Remainder is minus: Add back divisor 1: addl d1,d1 addl d6,d1 jcs 5f 2: addxl d0,d0 jcc 1b addl d1,d1 addl d6,d1 | Get 25th bit for rounding jra 7f | Remainder is plus: Subtract divisor 4: addl d1,d1 addl d5,d1 jcc 2b 5: addxl d0,d0 jcc 4b addl d1,d1 addl d5,d1 | Get 25th bit for rounding 7: jcs 1f | 54th bit set? addl d6,d1 | Nope, correct the remainder jra Lround$normal | Go finish up! 1: bset IMM (31),d1 | Set 54th bit jra Lround$normal | Go finish up! | a is either zero, denormalized, +/-INFINITY, or NaN Ldivsf$a$spec: jeq Ldivsf$a$big | Zero or denorm? movel IMM (0),d4 | Fixup to proper exponent andl IMM (0x7fffff),d1 | Strip off the false hidden bit jne Ldivsf$a$den | Denormalized? cmpl d0,d2 | Is b zero, denorm, INFINITY, NaN? jeq Ldivsf$b$big | Yep, go check it out jcs Ldivsf$zer | Normal number, return zero andl IMM (0x7fffff),d6 | Strip off the false hidden bit jeq Lf$inop | 0/0 is illegal! Ldivsf$zer: movel IMM (0),d0 | Return a zero jra Lf$sign | with proper sign Ldivsf$b$spec: jeq Ldivsf$b$big | Zero or denorm? movel IMM (0),d2 | Fixup to proper exponent andl IMM (0x7fffff),d6 | Strip off the false hidden bit jne Ldivsf$b$dck | Denormalized? movel IMM (SF_INFINITY),d0 | Return a +/-INFINITY addl a1,d0 movel IMM (DIVIDE_BY_ZERO),d7 jra $_exception_handler Ldivsf$a$big: andl IMM (0x7fffff),d1 | Remove false bit jeq Ldivsf$a$nf | Is a INFINITY? btst IMM (22),d1 | Signaling NaN? jeq Lf$a$snan | Yep, go report it cmpl d0,d2 | Is b zero, denorm, INFINITY, NaN? jne Lf$a | Nope, pass on NaN andl IMM (0x7fffff),d6 | Remove false bit jeq Lf$a | Is b INFINITY? jra Ldivsf$b$nan Ldivsf$a$nf: cmpl d0,d2 | Is b INFINITY, NaN? jne Lf$infty | Nope, return +/-INFINITY andl IMM (0x7fffff),d6 | Remove false bit jeq Lf$inop | INFINITY/INFINITY is invalid jra Ldivsf$b$nan Ldivsf$b$big: andl IMM (0x7fffff),d6 | Remove false bit jeq Ldivsf$zer | Is b INFINITY? Ldivsf$b$nan: btst IMM (22),d6 | Signaling NaN? jeq Lf$b$snan | Yep, go report it jra Lf$b | No signaling NaNs, just pass it on | If a number is denormalized we put an exponent of 1 but do not put the | hidden bit back into the fraction; instead we shift left until bit 21 | (the hidden bit) is set, adjusting the exponent accordingly. We do this | to ensure that the product of the fractions is close to 1. Ldivsf$a$den: subl IMM (1),d4 | Adjust exponent addl d1,d1 | shift a left until bit 20 is set btst IMM (SF_SIGNIF_DIG-1),d1 jeq Ldivsf$a$den jra Ldivsf$1 Ldivsf$b$dck: tstl d1 | Is a zero? jeq Ldivsf$zer | Yep, return zero Ldivsf$b$den: subql IMM (1),d2 | Adjust exponent addl d6,d6 | shift b left until bit 20 is set btst IMM (SF_SIGNIF_DIG-1),d6 jeq Ldivsf$b$den jra Ldivsf$2 |============================================================================= | __negsf2 |============================================================================= | This is trivial and could be shorter if we didn't bother checking for NaN ' | and +/-INFINITY. | float __negsf2(float); SYM (__negsf2): movel sp@(4),d0 | get number to negate in d0 bchg IMM (31),d0 | negate movel d0,d1 | make a positive copy bclr IMM (31),d1 | cmpl IMM (SF_INFINITY),d1 | compare to +INFINITY jgt 2f 1: rts 2: btst IMM (22),d0 | Signaling NaN? jne 1b | Nope, it's ok then ' link a6,IMM (-24) moveml d2-d7,sp@ pea NEGATE+OP1_FLOAT+DST_FLOAT jra Lf$a$snan | Illegal operaton on signaling NaN |============================================================================= | __cmpsf2 |============================================================================= .globl $_cmpsf GREATER = 1 LESS = -1 EQUAL = 0 | int __cmpsf2(float, float); SYM (__cmpsf2): lea GREATER,a1 | Unordered return value jra $_cmpsf | Go perform comparison $_cmpsf: link a6,IMM (-24) moveml d2-d7,sp@ pea COMPARE+OP1_FLOAT+OP2_FLOAT+DST_INT movel a6@(8),d0 | get first operand movel a6@(12),d1 | get second operand | Check if either is NaN, and in that case return garbage and signal | INVALID_OPERATION. Check also if either is zero, and clear the signs | if necessary. movel d0,d6 andl IMM (0x7fffffff),d0 jeq Lcmpsf$a$0 cmpl IMM (SF_INFINITY),d0 jhi Lcmpsf$a$nan Lcmpsf$1: movel d1,d5 andl IMM (0x7fffffff),d1 jeq Lcmpsf$b$0 cmpl IMM (SF_INFINITY),d1 jhi Lcmpsf$b$nan Lcmpsf$2: | Check the signs eorl d6,d5 jpl 1f | If the signs are not equal check if a >= 0 tstl d6 jpl Lcmpsf$a$gt$b | if (a >= 0 && b < 0) => a > b jmi Lcmpsf$b$gt$a | if (a < 0 && b >= 0) => a < b 1: | If the signs are equal check for < 0 tstl d6 jpl 1f | If both are negative exchange them movel d0,d5 movel d1,d0 movel d5,d1 1: | Now that they are positive we just compare them as longs (does this also | work for denormalized numbers?). cmpl d0,d1 jhi Lcmpsf$b$gt$a | |b| > |a| jne Lcmpsf$a$gt$b | |b| < |a| | If we got here a == b. movel IMM (EQUAL),d0 Lcmpsf$ret: moveml sp@(4),d2-d7 unlk a6 rts Lcmpsf$a$gt$b: movel IMM (GREATER),d0 jra Lcmpsf$ret Lcmpsf$b$gt$a: movel IMM (LESS),d0 jra Lcmpsf$ret Lcmpsf$a$0: bclr IMM (31),d6 jra Lcmpsf$1 Lcmpsf$b$0: bclr IMM (31),d5 jra Lcmpsf$2 Lcmpsf$a$nan: btst IMM (22),d0 | Is a a signaling NaN? jeq Lcmpsf$inv | Yep, go signal bclr IMM (31),d1 cmpl IMM (SF_INFINITY),d1 jls Lcmpsf$uno Lcmpsf$b$nan: btst IMM (22),d1 | Is B a signaling NaN? jeq Lcmpsf$inv Lcmpsf$uno: movel a1,d0 | Report no ordering for NaNs btst IMM (0),d0 | Error on unordered? jne Lcmpsf$ret Lcmpsf$inv: movel a1,d0 | Get unordered code asrl IMM (1),d0 | and dump invalid bit orl IMM (1),d0 jra Lf$inv | and report exception |============================================================================= | Rounding and Exception Handling |============================================================================= | The non-zero value to be rounded and returned is in registers d0-d1, | with the exponent in register d4, the sign in a1, and the operation type on | the top of the stack. The radix point is between bits 22 and 23 of d0, and | the exponent biased by 126 (adjusted before returning to 127). | | Three exceptions are checked for: Overflow, Underflow, and Inexact. | | Overflow is checked after rounding. If detected, and trapped, the trap | handler is passed the rounded value, with the exponent adjusted by a bias, | 192 in the case of FLOAT. The Inexact flag is set if the rounded value | is not exact. If not trapped, the value is set properly for the rounding | mode, and the Overflow and Inexact flags are set. The Inexact trap is | issued if enabled in this case. | | If the |value| <= 2^Emin (Emin = -126 for FLOAT) and the Underflow trap | is enabled, the value is normalize (if it isn't already), rounded, and sent ' | to the trap handler with the exponent bias adjusted. If the trap is not | enabled, the possibly denormalized value is rounded and, if inexact, the | Underflow flag is set. The Inexact flag is set if the result is not exact | in either case, and the Inexact trap taken if enabled, and the Underflow | trap is not. | We've got a normalized number, but need to denormalize it. ' denormal: lea SYM (_fpCCR),a0 btst IMM (1),a0@(TRAPE+3) | Check for Underflow trap enabled jne Lround$common | Yep, remain normalized cmpl IMM (-SF_SIGNIF_DIG-1),d4 | Really tiny? jls degen | Yep, just summarize movel IMM (32),d6 addl d6,d4 | Get left shift count subl d4,d6 | Get right shift count movel IMM (0),d2 | Summarize low order bits movew d1,d2 clrw d1 swap d1 orl d1,d2 movel d0,d1 lsll d4,d1 orl d2,d1 | Get back summary bits lsrl d6,d0 movel IMM (0),d4 | Show denormalized jra Lround$common degen: movel IMM (2),d1 | Show we're non-zero ' movel IMM (0),d0 | But all significance is gone movel IMM (0),d4 | Show denormalized jra Lround$common Lround$normal: | We first check for underflow. tstl d4 | Are we tiny? jlt denormal | Maybe, go check Lround$common: lea SYM (_fpCCR),a0 movel d1,d7 | Anything down here? jeq Lround$done | Exact! We're done! ' movel IMM (INEXACT_RESULT),d7 | Set Inexact flag movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lround$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Lround$to$minus jlt Lround$to$zero | jra Lround$to$plus Lround$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lround$done | Yep, truncate jra Lround$inc | Nope, increment Lround$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lround$done | Nope, truncate jra Lround$inc | Yep, increment Lround$to$nearest: | Now round: we do it as follows: after the shifting we can write the | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2. | If delta < 1, do nothing. If delta > 1, add 1 to f. | If delta == 1, we make sure the rounded number will be even (odd?) | (after shifting). addl d1,d1 | Is delta < 1? jcc Lround$done | Yep, do nothing jne Lround$inc | Is delta == 1? movel IMM (1),d3 | If so round to even andl d0,d3 | bit 1 is the last significant bit addl d3,d0 jra Lround$finish Lround$inc: addl IMM (1),d0 | Else add 1 Lround$finish: | Rounding may have changed our exponent. If so, adjust. btst IMM (SF_SIGNIF_DIG),d0 jeq Lround$done lsrl IMM (1),d0 addl IMM (1),d4 | Fixup exponent Lround$to$zero: Lround$done: | Now we're properly normalized and rounded, or, if the Underflow trap ' | is not enabled, maybe denormalized and rounded. If we're tiny ' | (|value| <= 2^Emin) the Tiny flag (proposed Underflow flag) is set. | If the rounded value is not exact, the Inexact flag is set. movel a0@(TRAPE),d5 | Get trap enables cmpl IMM (SF_MAX_EXP-1),d4 | Overflow? jge Ltrap$overflow | Nope, go check for Underflow tstl d4 | Check for Underflow jle Ltrap$underflow Ltrap$check: | Check for trapped exceptions. First we put the number back together. swap d4 | Put the exponent in place lsll IMM (7),d4 addl d4,d0 | and add one to exponent tstl d7 | Any exceptions? jne Ltrap$assemble | Yep, go handle addl a1,d0 | Nope, add in sign moveml sp@(4),d2-d7 | We're done ' | XXX if frame pointer is ever removed, stack pointer must | be adjusted here. unlk a6 | and return rts Ltrap$underflow: jlt 2f | Really small (& trap enabled)? cmpl IMM (0x800000),d0 | |value| <= 2^Emin? jgt Ltrap$check | Nope, return btst IMM (1),d5 | Underflow trap enabled? jne 2f | Yep, adjust the exponenet btst IMM (0),d7 | Result Inexact? jeq Ltrap$check | Nope, won't set Underflow ' jra 3f 2: addl IMM (192), d4 | Adjust exponent 3: addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Ltrap$check Ltrap$overflow: addl IMM (OVERFLOW),d7 | Set the Overflow flag btst IMM (2),d5 | Overflow trap enabled? jeq 1f | Nope, go adjust value subl IMM (192),d4 | Adjust exponent jra Ltrap$check | and go trap 1: movel IMM (INEXACT_RESULT+OVERFLOW),d7 | Inexact too movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Loflow$nearest | Yep, return infinity subl IMM (ROUND_TO_PLUS),d6 jhi Loflow$minus jlt Loflow$zero | jra Loflow$plus Loflow$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Loflow$nearest | Nope, go to infinity jra Loflow$zero | Yep, just go big Loflow$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Loflow$nearest | Yep, go to infinity jra Loflow$zero | Nope, just go big Loflow$nearest: movel IMM (0),d0 | Build an infinity value movel d0,d1 movel IMM (SF_MAX_EXP),d4 jra Ltrap$check Loflow$zero: movel IMM (0xffffff),d0 | Build max value movel IMM (SF_MAX_EXP-2),d4 jra Ltrap$check Ltrap$assemble: addl a1,d0 | Put the sign back jra $_exception_handler Lf$a$snan: movel a6@(8),d0 | Return a jra Lf$snan Lf$b$snan: movel a6@(12),d0 | Return b Lf$snan: bset IMM (22),d0 | Change to Quiet NaN jra Lf$inv Lf$inop: | Invalid Operation. Return a quiet NaN and set the exception flags movel IMM (Q_NaN),d0 Lf$inv: movel IMM (INVALID_OPERATION),d7 jra $_exception_handler Lf$infty: | Return a properly signed and set the exception flags movel IMM (SF_INFINITY),d0 Lf$sign: addl a1,d0 | Sign bit Lf$return: moveml sp@(4),d2-d7 | XXX if frame pointer is ever removed, stack pointer must | be adjusted here. unlk a6 | and return rts Lf$a: movel a6@(8),d0 | Return a jra Lf$return Lf$b: movel a6@(12),d0 | Return b jra Lf$return #endif /* L_float */ | gcc expects the routines __eqdf2, __nedf2, __gtdf2, __gedf2, | __ledf2, __ltdf2 to all return the same value as a direct call to | __cmpdf2 would. | Conditionals now return the same values as the routines in fp-bit.c | when the operands are unordered. They do this by passing in a1 the | unordered return value to the internal compare routine $_cmpdf. | If this value is *2 (bit 0 is 0) and the result is unordered, an | Invalid Operation exception occurs. #ifdef L_eqdf2 .text .proc .globl $_cmpdf .globl SYM (__eqdf2) SYM (__eqdf2): lea 1,a1 | Unordered return value jra $_cmpdf | Go perform comparison #endif /* L_eqdf2 */ #ifdef L_nedf2 .text .proc .globl $_cmpdf .globl SYM (__nedf2) SYM (__nedf2): lea 1,a1 | Unordered return value jra $_cmpdf | Go perform comparison #endif /* L_nedf2 */ #ifdef L_gtdf2 .text .proc .globl $_cmpdf .globl SYM (__gtdf2) SYM (__gtdf2): lea -2,a1 | Unordered return value jra $_cmpdf | Go perform comparison #endif /* L_gtdf2 */ #ifdef L_gedf2 .text .proc .globl $_cmpdf .globl SYM (__gedf2) SYM (__gedf2): lea -2,a1 | Unordered return value jra $_cmpdf | Go perform comparison #endif /* L_gedf2 */ #ifdef L_ltdf2 .text .proc .globl $_cmpdf .globl SYM (__ltdf2) SYM (__ltdf2): lea 2,a1 | Unordered return value jra $_cmpdf | Go perform comparison #endif /* L_ltdf2 */ #ifdef L_ledf2 .text .proc .globl $_cmpdf .globl SYM (__ledf2) SYM (__ledf2): lea 2,a1 | Unordered return value jra $_cmpdf | Go perform comparison #endif /* L_ledf2 */ | The comments above about __eqdf2, et. al., also apply to __eqsf2, | et. al., except that the latter call __cmpsf2 rather than __cmpdf2. #ifdef L_eqsf2 .text .proc .globl $_cmpsf .globl SYM (__eqsf2) SYM (__eqsf2): lea 1,a1 | Unordered return value jra $_cmpsf | Go perform comparison #endif /* L_eqsf2 */ #ifdef L_nesf2 .text .proc .globl $_cmpsf .globl SYM (__nesf2) SYM (__nesf2): lea 1,a1 | Unordered return value jra $_cmpsf | Go perform comparison #endif /* L_nesf2 */ #ifdef L_gtsf2 .text .proc .globl $_cmpsf .globl SYM (__gtsf2) SYM (__gtsf2): lea -2,a1 | Unordered return value jra $_cmpsf | Go perform comparison #endif /* L_gtsf2 */ #ifdef L_gesf2 .text .proc .globl $_cmpsf .globl SYM (__gesf2) SYM (__gesf2): lea -2,a1 | Unordered return value jra $_cmpsf | Go perform comparison #endif /* L_gesf2 */ #ifdef L_ltsf2 .text .proc .globl $_cmpsf .globl SYM (__ltsf2) SYM (__ltsf2): lea 2,a1 | Unordered return value jra $_cmpsf | Go perform comparison #endif /* L_ltsf2 */ #ifdef L_lesf2 .text .proc .globl $_cmpsf .globl SYM (__lesf2) SYM (__lesf2): lea 2,a1 | Unordered return value jra $_cmpsf | Go perform comparison #endif /* L_lesf2 */ #ifdef L_extended | Entry points: .globl SYM (__addxf3) .globl SYM (__subxf3) .globl SYM (__mulxf3) .globl SYM (__divxf3) .globl SYM (__negxf2) .globl SYM (__cmpxf2) .globl SYM (__eqxf2) .globl SYM (__nexf2) .globl SYM (__gtxf2) .globl SYM (__gexf2) .globl SYM (__ltxf2) .globl SYM (__lexf2) .globl SYM (__extenddfxf2) .globl SYM (__truncxfdf2) .globl SYM (__extendsfxf2) .globl SYM (__truncxfsf2) .globl SYM (__floatsixf) .globl SYM (__fixxfsi) .text .even |============================================================================= |============================================================================= | extended precision routines |============================================================================= |============================================================================= | An extended precision floating point number (long double) has the format: | | struct _long_double { | unsigned int sign : 1; /* sign bit */ | unsigned int exponent : 15; /* exponent, shifted by 16383 */ | unsigned int z : 16; /* unused bits */ | unsigned int fraction : 64; /* fraction */ | } long_double; | | Thus sizeof(long_double) = 12 (96 bits). | | All the routines are callable from C programs, and return the result | at the location addressed by a0. They also preserve all registers except | d0-d1 and a0-a1. |============================================================================= | __subxf3 |============================================================================= | double __subdf3(double, double); SYM (__subxf3): link a6,IMM (-24) moveml d2-d7,sp@ pea SUB+OP1_LDOUBLE+OP2_LDOUBLE+DST_LDOUBLE movel a2,sp@- movel a6@(20),d6 | Get second operand's exponent & sign ' bchg IMM (31),d6 | Ghange sign of second operand jra 1f | and fall through, so we always add |============================================================================= | __addxf3 |============================================================================= | long double __addxf3(long double, long double); SYM (__addxf3): link a6,IMM (-24) moveml d2-d7,sp@ pea ADD+OP1_LDOUBLE+OP2_LDOUBLE+DST_LDOUBLE movel a2,sp@- movel a6@(20),d6 | Get second operand's exponent & sign ' 1: movel a6@(8),d0 | Get first operand movel a6@(12),d1 movel a6@(16),d2 movel a6@(24),d4 | Get second operand movel a6@(28),d5 clrw d0 | Zap unused bits clrw d6 movel d0,a1 | Get a's sign bit in a1 ' movel d6,a2 | Get b's sign in a2 ' bclr IMM (31),d0 | Clear signs from working regs bclr IMM (31),d6 subl d0,a1 subl d6,a2 swap d0 | Right justify exponents swap d6 movel IMM (XF_MAX_EXP),d3 | Special case literal cmpl d3,d0 | +/-INFINITY, NaN? jeq Laddxf$a$nf | Yep, go handle cmpl d3,d6 | +/-INFINITY, NaN? jeq Laddxf$b$nf | Yep, go handle tstl d1 jeq Laddxf$a$mbz | May be zero? Laddxf$1: tstl d4 jeq Laddxf$b$mbz | May be zero? Laddxf$2: | Both operands are finite and non-zero. Here we shift the numbers until | their exponents are the same. If the operand with the greater exponent | is unnormalized, it is normalized, and the exponents are compared again. cmpl d0,d6 | compare the exponents jeq Laddxf$3 | If equal we don't need to shift ' jmi 9f | Branch if second exponent is lower tstl d4 | Is b normalized? jpl Laddxf$b$norm | Nope, go normalize it subl d0,d6 | Put largest exponent into d0 addl d6,d0 addl a2,a1 | Swap a & b addl a1,a2 addl a2,a1 movel d1,d3 | Swap high words movel d4,d1 movel d3,d4 movel d2,d3 | Swap low words movel d5,d2 movel d3,d5 jra 1f | Here we have a's exponent larger than b's, so we have to shift b (maybe). 9: tstl d1 | Is a normalized? jpl Laddxf$a$norm | Nope, go normalize it subl d0,d6 negl d6 | Get positive right-shift count | If difference is too large we don't shift, just summarize b's | non-zeroness in d6 1: addl d0,a2 | Save exponent movel IMM (32), d0 cmpl d0,d6 | If difference >= 32, shift by longs jge 5f subl d6,d0 | Get left shift count movel d5,d3 lsll d0,d3 lsrl d6,d5 movel d4,d7 lsll d0,d7 orl d7,d5 lsrl d6,d4 jra 3f 5: subl d0,d6 | Get right-shift count cmpl d0,d6 | Shift still more? jge 6f subl d6,d0 | Get left-shift count movel d5,d7 movel d5,d3 lsrl d6,d3 movel d4,d5 lsll d0,d4 orl d4,d3 lsrl d6,d5 lsll d0,d7 jeq 2f bset IMM (0),d3 | Set sticky bit jra 2f 6: subl d0,d6 | Get right-shift count cmpl d0,d6 | Really small? jge Laddxf$b$small | Yep, just summarize it subl d6,d0 | Get left-shift count movel d4,d3 lsrl d6,d3 | This is all the significance left lsll d0,d4 | Anything down here? jne 1f | Yep, go flag it tstl d5 | How about here? jeq 2f | Nope, leave it alone 1: bset IMM (0),d3 | Set sticky bit movel IMM (0),d5 | Zap vacated words 2: movel IMM (0),d4 3: movew a2,d0 | Restore exponent subl d0,a2 jra Laddxf$4 | If one of the numbers is pretty small (difference of exponents >= | DF_SIGNIF_DIG+1) we just summarize it to show an Inexact operation. Laddxf$b$small: movel IMM (1),d3 | Show b is piddling movel IMM (0),d4 movel d4,d5 movew a2,d0 | Restore exponent subl d0,a2 jra Laddxf$4 | Now go do operation & rounding Laddxf$3: movel IMM (0),d3 | Equal exponents; zap extension word Laddxf$4: | Now we have the numbers in d1-d2 and d4-d5/d3, the exponent in d0, and | the signs in a0,a1. Either d1-d2 is normalized, or the exponents are | equal. | Here we have to decide whether to add or subtract the numbers: cmpl a2,a1 | Compare signs jne Lsubxf$0 | If the signs are different we get | to subtract addl d5,d2 | else we add the numbers addxl d4,d1 | | Before rounding normalize so bit #DF_SIGNIF_DIG-1 is set (we will consider | the case of denormalized numbers in the rounding routine itself). | As in the addition (not in the subtraction!) we could have set | one more bit we check this: jcc Laddxf$unnormal movel IMM (31),d5 lsrl IMM (1),d3 subxl d6,d6 | Preserve sticky bit negl d6 orl d6,d3 movel d2,d4 lsll d5,d4 orl d4,d3 lsrl IMM (1),d2 movel d1,d4 lsll d5,d4 orl d4,d2 lsrl IMM (1),d1 bset d5,d1 addl IMM (1),d0 jra Lround$common | Go finish up! Lsubxf$0: | Here we do the subtraction. negl d3 subxl d5,d2 subxl d4,d1 jeq Laddxf$ezer | if zero just exit jcc Laddxf$unnormal | if positive skip the following movel a2,a1 | Use b's sign ' negl d2 | and negate result negxl d1 | We don't need to touch d3 as ' | the exponents are equal. d3 | is therfore zero | Normalization for add/subtract | When the exponents of the operands differ by more than one, there will | at most be one left shift. This means that whenever there is more than | one left shift, the result is exact, and rounding may be skipped. Laddxf$unnormal: tstl d1 | Already normalized? jmi Lround$common lea SYM (_fpCCR),a2 btst IMM (1),a2@(TRAPE+3) | Check for Underflow trap enabled jne Laddxf$normal | Yep, don't check for denorm ' tstl d0 | Already min exponent? jeq Laddxf$exact | Yep, don't shift ' | Find the MSB of the value. At most one bit will be shifted from d3 into d2. subl IMM (1),d0 | Min exponent yet? jeq Laddxf$exact1 | Yep, but at least it's exact! ' addl d3,d3 | Perhaps shift last bit out addxl d2,d2 addxl d1,d1 jmi Lround$common | or perhaps not? | Left shift until it's either normalized, or we've reached the denormal | exponent value 1: subl IMM (1),d0 | Denormal yet? jeq Laddxf$exact2 | Yep, go adjust other words | (X is zero if branch is taken) addl d2,d2 | Try moving over one place addxl d1,d1 | Normalized yet? jpl 1b | Nope, try another place Laddxf$exact: swap d0 | All we need to do now is | assemble the extended float! Laddxf$ret: addl a1,d0 | Put in the sign Laddxf$pzer: movel d0,a0@ movel d1,a0@(4) movel d2,a0@(8) movel a0,d0 movel sp@+,a2 moveml sp@(4),d2-d7 | We're done! ' unlk a6 rts Laddxf$exact1: addl d3,d3 | Perhaps shift last bit out Laddxf$exact2: addxl d2,d2 | Try moving over one place addxl d1,d1 | Normalized yet? jra Laddxf$exact | As the Underflow trap is enabled, we just normalize the result | without checking for the denormal exponent value. Laddxf$normal: subl IMM (1),d0 | Adjust exponent addl d3,d3 | Perhaps shift last bit out addxl d2,d2 addxl d1,d1 jmi Lround$common | Normalized yet? | Left shift until it's normalized ' 1: subl IMM (1),d0 | Adjust exponent addl d2,d2 | Try moving over one place addxl d1,d1 | Normalized yet? jpl 1b | Nope, try another place tstl d0 | Underflow? jgt Laddxf$exact | Nope, assemble & return jlt 2f | Yep, go report movel d1,d4 addl d4,d4 orl d2,d4 | 2^Emin? jne Laddxf$exact 2: addl IMM (24576),d0 | Adjust exponent | All we need to do now is swap d0 | assemble the double float! addl a1,d0 | Put in the exponent and the sign movel IMM (UNDERFLOW),d7 | Show error movel sp@+,a2 jra $_exception_handler | and go trap | a may be zero Laddxf$a$mbz: tstl d2 | Really zero? jne Laddxf$1 | Nope, continue movel a6@(24),d4 | Get second operand movel a6@(28),d2 movel d4,d1 orl d2,d4 | Zero also? jne Laddxf$b | Nope, then return it movel IMM (0),d0 | Fixup to proper exponent cmpl a2,a1 | Compare signs jeq Laddxf$ret | and return zero if same Laddxf$ezer: movel IMM (0),d0 lea SYM (_fpCCR),a2 | Signs differ, movel a2@(ROUND),d4 | Rounding mode determines sign subl IMM (ROUND_TO_MINUS),d4 jne Laddxf$pzer | Not round-to-minus, return +0 bset IMM (31),d0 | Return -0 jra Laddxf$pzer | b may be zero Laddxf$b$mbz: tstl d5 | Really zero? jne Laddxf$2 | Nope, continue jra Laddxf$a | Return a as result | One of the operands is zero. Return the other, but be carefull, as | we still may need to check for Underflow! If the other operand is | non-zero and is tiny, i.e. |value| <= 2^Emin, and the Underflow trap | is enabled, the number needs to be normalized, and sent off to the | Underflow trap routine. Laddxf$b: | Return b (if a is zero) movel a2,a1 | Use b's sign ' movel d6,d0 | and exponent Laddxf$a: movel IMM (0),d3 | Zero the extension word jra Laddxf$unnormal | Go finish up! | a is either NaN or +/-INFINITY Laddxf$a$nf: movel d1,d7 addl d7,d7 orl d7,d2 | +/-INFINITY? jeq 1f tstl d7 | Signaling NaN? jpl Lx$a$snan | Yep, go signal 1: cmpl d3,d6 | Is b +/-INFINITY, NaN? jeq Laddxf$b$nf | Yep, go check it out tstl d2 | Return jeq Laddxf$infty | +/-INFINITY jra Lx$a | or NaN | b is either NaN or +/-INFINITY, a may be also, though not | a signaling NaN Laddxf$b$nf: addl a2,a1 | Swap signs addl a1,a2 addl a2,a1 movel d4,d7 addl d7,d7 orl d7,d5 | +/-INFINITY? jeq 1f tstl d7 | Signaling NaN? jpl Lx$b$snan | Yep, go signal jra Lx$b | Just pass on non-signaling NaN | b is +/-INFINITY, what is a? 1: cmpl d3,d0 | +/-INFINITY, NaN? jne Laddxf$infty | Nope, return b's +/-INFINITY ' tstl d2 | Is a +/-INFINITY? jne Lx$a | Nope, must be NaN | We know that both are +/-INFINITY. It is not valid if their signs | are different cmpl a2,a1 | Compare sign bits jne Lx$inop | Can't subtract INFINITYs! ' Laddxf$infty: jra Lx$infty | +/-INFINITY | Normalize a Laddxf$a$norm: subl IMM (1),d0 | Adjust exponent jeq 2f | To denormal boundary? addl d2,d2 | Shift a over one place addxl d1,d1 jpl Laddxf$a$norm | Normalized yet? jra Laddxf$2 | Go compare exponents again 2: addl d2,d2 | Shift a over one more place addxl d1,d1 jra Laddxf$2 | Normalize b Laddxf$b$norm: subl IMM (1),d6 | Adjust exponent jeq 2f | To denormal boundary? addl d5,d5 | Shift b over one place addxl d4,d4 jpl Laddxf$b$norm | Normalized yet? jra Laddxf$2 | Go compare exponents again 2: addl d5,d5 | Shift b over one more place addxl d4,d4 jra Laddxf$2 |============================================================================= | __mulxf3 |============================================================================= | long double __mulxf3(long double, long double); SYM (__mulxf3): link a6,IMM (-24) moveml d2-d7,sp@ pea MULTIPLY+OP1_LDOUBLE+OP2_LDOUBLE+DST_LDOUBLE movel a2,sp@- movel a6@(8),d0 | Get first operand's sign & exponent ' movel a6@(20),d5 | Get second operand's sign & exponent ' clrw d0 | Zap unused bits clrw d5 movel d0,a1 | Get a's sign bit in a1 ' bclr IMM (31),d0 | Clear sign from working reg subl d0,a1 | & isolate sign addl d5,a1 | Get sign of result bclr IMM (31),d5 | Clear sign from working reg subl d5,a1 | & isolate sign swap d0 | Right justify exponents swap d5 movel IMM (XF_MAX_EXP),d6 | Special case literal cmpl d6,d0 | Zero, denorm, INFINITY, NaN? jeq Lmulxf$a$big | Yep, go handle cmpl d6,d5 | Zero, denorm, INFINITY, NaN? jeq Lmulxf$b$big | Yep, go handle movel a6@(16),d2 | Get first operand movel a6@(12),d1 | jpl Lmulxf$a$mbz | Unnormalized? Lmulxf$1: movel a6@(28),d4 | Get second operand movel a6@(24),d3 jpl Lmulxf$b$mbz | Unnormalized? Lmulxf$2: addl d5,d0 | add exponents subl IMM (XF_BIAS),d0 | and subtract bias | We are now ready to do the multiplication. The situation is as follows: | both a and b have bit 63 ( bit 31 of d1 and d3) set (even if they were | unnormalized to start with!). The final product has | either bit 30 or 31 of the most significant word as the MSB. At that | point we must either left shift one more place, or increment the | exponent. | If one of the operands has all 32 low bits zero, we choose it to be | the multiplier. This speeds up products of many small numbers. movel d0,a2 | a2 preserves the exponent movel d2,d7 | Low word of a zero? jeq 1f movel d4,d5 | Low a non-zero, how about b? jne 2f | Nope, don't switch operands ' movel d1,d6 | Low word of b is zero, movel d2,d5 | Swap operands movel d3,d1 movel d4,d7 | Zero in low word of multiplier movel d6,d3 movel d2,d4 jra 4f 1: movel d4,d5 2: movel d3,d6 4: | Instead of having a seperate counter, we float a flag bit through | the multiplier word. movel IMM (0),d0 | Literal for carry propagation movel d0,d2 | Zero upper word of product so far addl d1,d1 | Bump off MSB, addl IMM (1),d1 | set "count" bit addl d1,d1 | Shift multiplier | Here we get to do the product. 1: jcc 2f | Next bit set? addl d4,d4 | Shift left one place addxl d3,d3 addxl d2,d2 addl d5,d4 | and add in b addxl d6,d3 | and carry to third word addxl d0,d2 | (we're not yet into the 4th) ' addl d1,d1 | Shift multiplier jne 1b | Any bits left? jra 3f 2: addl d4,d4 | Shift left one place addxl d3,d3 addxl d2,d2 addl d1,d1 | Shift multiplier jne 1b | Any bits left? | We've finished the first 32 bits of the multiplier. Now on to the ' | last 32. 3: tstl d7 | More work to do? jeq 3f | Nope, just shift words addxl d7,d7 | Shift out next bit, | and shift in flag bit! 1: jcc 2f | Next bit set? addl d4,d4 | Shift left one place addxl d3,d3 addxl d2,d2 addxl d1,d1 addl d5,d4 | and add in b addxl d6,d3 addxl d0,d2 | and carry to third word addxl d0,d1 | and to forth word addl d7,d7 | Shift multiplier jne 1b | Any bits left? jra 4f 2: addl d4,d4 | Shift left one place addxl d3,d3 addxl d2,d2 addxl d1,d1 addl d7,d7 | Shift multiplier jne 1b | Any bits left? 4: movel a2,d0 | Restore exponent tstl d1 | Aligned? jmi 1f | Yep, go bump exponent addl d4,d4 | Nope, align addxl d3,d3 addxl d2,d2 addxl d1,d1 | Now bit 31 is set 5: movel d4,d6 | Now set the sticky bit negl d6 orl d6,d4 lsrl IMM (1),d4 orl d4,d3 jra Lround$normal | Go finish up! 3: movel d2,d1 | Low 32 bits of multiplier movel d3,d2 | are zero, we just need movel d4,d3 | to left shift 32 bits movel d0,d4 jra 4b 1: addl IMM (1),d0 | Increase exponent jra 5b | Go finish up! Lmulxf$zer: movel IMM (0),d1 | Return a zero movel d1,d2 movel d1,d0 jra Lx$sign | with proper sign | a may be zero Lmulxf$a$mbz: jne Lmulxf$a$den | Unnormalized? tstl d2 | or zero? jeq Lmulxf$zer | Yep, return zero Lmulxf$a$den: subl IMM (1),d0 | Adjust exponent addl d2,d2 | shift a left until bit 31 is set addxl d1,d1 | jpl Lmulxf$a$den jra Lmulxf$1 | Lmulxf$b$mbz: jne Lmulxf$b$den | Unnormalized? tstl d4 | or zero? jeq Lmulxf$zer | Yep, return zero Lmulxf$b$den: subl IMM (1),d5 | and adjust exponent addl d4,d4 | shift b left until bit 31 is set addxl d3,d3 | jpl Lmulxf$b$den jra Lmulxf$2 | Lmulxf$a$big: movel a6@(12),d1 | Get first operand addl d1,d1 movel d1,d2 orl a6@(16),d2 | Is a +/-INFINITY? jeq Lmulxf$a$nf tstl d1 | Signaling NaN? jpl Lx$a$snan | Yep, go report it cmpl d6,d5 | Is b +/-INFINITY, NaN? jne Lx$a | Nope, pass on NaN movel a6@(24),d3 addl d3,d3 movel d3,d4 orl a6@(28),d4 | Is b INFINITY? jeq Lx$a | Yep, pass on NaN jra Lmulxf$b$nan Lmulxf$a$nf: movel a6@(24),d3 | Get b cmpl d6,d5 | Is b +/-INFINITY, NaN? jeq Lmulxf$ab$big | Yep, go check it out orl a6@(28),d3 | Is b zero? jeq Lx$inop | Yep, Illegal operation! jra Lx$infty | Nope, return +/-INFINITY Lmulxf$b$big: movel a6@(24),d3 addl d3,d3 movel d3,d4 orl a6@(28),d4 | Is b INFINITY? jne Lmulxf$b$nan orl d1,d2 | Is a zero? jne Lx$infty | Nope, return +/-INFINITY jra Lx$inop | Yep, Illegal operation! Lmulxf$ab$big: addl d3,d3 movel d3,d4 orl a6@(28),d4 | Is b INFINITY? jeq Lx$infty | Yep, return +/-INFINITY Lmulxf$b$nan: tstl d3 | Signaling NaN? jpl Lx$b$snan | Yep, go report it jra Lx$b | No signaling NaNs, just pass it on |============================================================================= | __divxf3 |============================================================================= | long double __divdf3(long double, long double); SYM (__divxf3): link a6,IMM (-24) moveml d2-d7,sp@ pea DIVIDE+OP1_LDOUBLE+OP2_LDOUBLE+DST_LDOUBLE movel a2,sp@- movel a6@(8),d0 | Get a's sign and exponent ' movel a6@(20),d2 | and b's sign and exponent ' clrw d0 | Zap unused bits clrw d2 movel d0,a1 | Get a's sign bit in a1 ' bclr IMM (31),d0 | Clear sign from working reg subl d0,a1 addl d2,a1 | Get sign of result bclr IMM (31),d2 | Clear sign from working reg subl d2,a1 | Bit 31 is sign of result swap d0 | Right justify exponents swap d2 movel IMM (XF_MAX_EXP),d1 | Special case literal cmpl d1,d0 | +/-INFINITY, NaN? jeq Ldivxf$a$big | Yep, go handle cmpl d1,d2 | +/-INFINITY, NaN? jeq Ldivxf$b$big | Yep, go handle movel a6@(16),d4 | Get a in d3-d4 movel a6@(12),d3 jpl Ldivxf$a$mbz | Unnormalized? Ldivxf$1: movel a6@(28),d7 | Get b in d6-d7 movel a6@(24),d6 jpl Ldivxf$b$mbz | Unnormalized? Ldivxf$2: subl d2,d0 | Subtract divisor exponent addl IMM (XF_BIAS),d0 | and add back bias | We add -divisor instead of subtracting as that way we can use the | carry bit directly to build up the quotient. The initial dividend, | and final remainder, is in d3,d4. The quotient is built up in d2. | When the first 32 bits are complete, they are moved to d1, and the | next 32 bits are loaded into d2. Then a 65th bit is created in bit | 31 of d3. | | Instead of using a seperate counter, we float a flag bit through | d2, and loop until it pops into the C flag. We zero d1 initially, | and use that to flag when we've got all 64 bits of the quotient. ' subl d7,d4 | Trial subtract to find first subxl d6,d3 | set bit subxl d2,d2 movel IMM (0),d1 | Still need low 32 quotient bits lsrl IMM (1),d7 | Halve the divisor jcs Ldivxf$odd | Odd? | Divisor is even lsrl IMM (1),d6 jcc 1f bset IMM (31),d7 1: movel d7,a2 | Save plus low word negl d7 | Form minus low word movel d6,d5 negxl d5 | Form minus high word tstl d2 | First quotient bit yet? jeq 1f movel IMM (1),d2 | Set counter flag subl IMM (1),d0 | Adjust exponent jra 11f | and enter add back loop 1: movel IMM (3),d2 | Set counter flat & first bit jra 12f | and go enter subract loop | Remainder is minus: Add back divisor 1: addl d4,d4 | Shift remainder, shift in zero addxl d3,d3 11: addl a2,d4 | Add back divisor addxl d6,d3 jcs 5f | This quotient bit set? 2: addxl d2,d2 | Nope, do another add back jcc 1b | Done with this word? tstl d1 | Done 64 bits yet? jne 3f movel d2,d1 | Nope, 32 more to go movel IMM (1), d2 jra 1b 3: addl d4,d4 addxl d3,d3 addl a2,d4 | Get 65th bit for rounding addxl d6,d3 jra 7f | Remainder is plus: Subtract divisor 4: addl d4,d4 | Shift remainder, shift in zero addxl d3,d3 12: addl d7,d4 | Subtract divisor addxl d5,d3 jcc 2b | This quotient bit set? 5: addxl d2,d2 | Nope, do another subtract jcc 4b | Done with this word? tstl d1 | Done 64 bits yet? jne 6f movel d2,d1 | Nope, 32 more to go movel IMM (1), d2 jra 4b 6: addl d4,d4 addxl d3,d3 addl d7,d4 | Get 65th bit for rounding addxl d5,d3 7: jcs 1f | 65th bit set? Ldivxf$remainder: subl d7,d4 | Nope, subxl d5,d3 | correct the remainder jra 2f 1: bset IMM (31),d3 | Set 65th bit 2: movel d4,d5 | Set sticky bit negl d5 orl d5,d4 lsrl IMM (1),d4 orl d4,d3 jra Lround$normal | Go finish up! | Divisor is odd Ldivxf$odd: lsrl IMM (1),d6 jcc 1f bset IMM (31),d7 1: movel d7,a2 | Save plus low word notl d7 | Form minus low word movel d6,d5 notl d5 | Form minus high word tstl d2 | To first set quotient bit yet? jeq 1f movel IMM (-2),d2 | Nope, set counter flag subl IMM (1),d0 | Adjust exponent jra 12f | and enter add back loop 1: movel IMM (-3),d2 | Yep, set bit & counter flag jra 13f | & enter subtract loop 11: addl d2,d2 | (Set X flag) | Remainder is minus: Add back divisor 1: addxl d4,d4 | Shift remainder, shift in one addxl d3,d3 12: addl a2,d4 | Add back divisor addxl d6,d3 jcs 5f | This quotient bit set? 2: addxl d2,d2 | Nope, do another add back jcs 1b | Done with this word? tstl d1 | Done 64 bits yet? jne 3f movel d2,d1 | Nope, 32 more to go movel IMM (-1), d2 jra 11b 3: addl d4,d4 addxl d3,d3 addl IMM (1),d4 addl a2,d4 | Get 65th bit for rounding addxl d6,d3 jra 7f | Remainder is plus: Subtract divisor 4: addxl d4,d4 | Shift remainder, shift in one addxl d3,d3 13: addl d7,d4 | Subtract divisor addxl d5,d3 jcc 2b | This quotient bit set? 5: addxl d2,d2 | Nope, do another subtract jcs 4b | Done with this word? tstl d1 | Done 64 bits yet? jne 6f movel d2,d1 | Nope, 32 more to go movel IMM (-1), d2 addl d2,d2 | (Set X flag) jra 4b 6: addl d4,d4 addxl d3,d3 addl IMM (1),d4 addl d7,d4 | Get 65th bit for rounding addxl d5,d3 7: jcc Ldivxf$remainder | 65th bit set? movel IMM (0xc0000000),d3 | Set 65th bit, n/z remainder jra Lround$normal | Go finish up! | a may be zero Ldivxf$a$mbz: jne Ldivxf$a$den | Unnormalized? tstl d4 | or zero? jne Ldivxf$a$den movel a6@(16),d6 | Get b orl a6@(20),d6 | Is it zero? jeq Lx$inop | Yep, illegal operation: 0/0 Ldivxf$zer: movel IMM (0),d1 | Return a zero movel d1,d2 movel d1,d0 jra Lx$sign | with proper sign Ldivxf$a$den: subl IMM (1),d0 | and adjust exponent addl d4,d4 | shift a left until bit 31 is set addxl d3,d3 jpl Ldivxf$a$den jra Ldivxf$1 Ldivxf$b$mbz: jne Ldivxf$b$den | Unnormalized? tstl d7 | or zero? jeq Ldivxf$b$zer Ldivxf$b$den: subl IMM (1),d2 | Adjust exponent addl d7,d7 | Shift b left until bit 31 is set addxl d6,d6 | jpl Ldivxf$b$den jra Ldivxf$2 | Ldivxf$b$zer: movel IMM (XF_INFINITY),d0 | Return a signed INFINITY addl a1,d0 movel IMM (0),d1 movel d1,d2 movel IMM (DIVIDE_BY_ZERO),d7 movel sp@+,a2 jra $_exception_handler Ldivxf$a$big: movel a6@(12),d3 | Get a addl d3,d3 movel d3,d4 orl a6@(16),d4 | Is a +/-INFINITY? jeq Ldivxf$a$nf tstl d3 | Signaling NaN? jpl Lx$a$snan | Yep, go report it cmpl d1,d2 | Is b +/-INFINITY, NaN? jne Lx$a | Nope, pass on NaN unchanged movel a6@(24),d6 | Get b addl d6,d6 movel d6,d7 orl a6@(28),d7 | Is b INFINITY? jeq Lx$a | Yep, pass on NaN jra Ldivxf$b$nan Ldivxf$a$nf: cmpl d1,d2 | Is b +/-INFINITY, NaN? jne Lx$infty | Nope, return +/-INFINITY movel a6@(24),d6 | Get b addl d6,d6 movel d6,d7 orl a6@(28),d7 | Is b INFINITY? jeq Lx$inop | Yep, INFINITY/INFINITY is invalid jra Ldivxf$b$nan Ldivxf$b$big: movel a6@(24),d6 | Get b addl d6,d6 movel d6,d7 orl a6@(28),d7 | Is b INFINITY? jeq Ldivxf$zer | Yep, return zero Ldivxf$b$nan: tstl d6 | Signaling NaN? jpl Lx$b$snan | Yep, go report it jra Lx$b | No signaling NaNs, just pass it on |============================================================================= | __negxf2 |============================================================================= | long double __negxf2(long double); SYM (__negxf2): movel sp@(4),d0 | Get number to negate in d0-d1 clrw d0 | Zap unused bits bchg IMM (31),d0 | Negate movel d0,d1 | Make a positive copy (for the tests) bclr IMM (31),d1 cmpl IMM (XF_INFINITY),d1 | Compare to +INFINITY jge 2f | NaN or +/-INFINITY, check further 1: movel sp@(8),d1 | Get middle word movel d0,a0@ movel d1,a0@(4) movel sp@(12),d1 movel d1,a0@(8) movel a0,d0 rts 2: movel sp@(8),d1 | Get middle word addl d1,d1 jne 3f | Go check NaN type tstl sp@(12) jeq 1b | Return +/-INFINITY tstl d1 | Signaling NaN? 3: jmi 1b | Nope, just pass it on link a6,IMM (-24) moveml d2-d7,sp@ pea NEGATE+OP1_LDOUBLE+DST_LDOUBLE movel a2,sp@- jra Lx$a$snan |============================================================================= | __cmpxf2 |============================================================================= .globl $_cmpxf GREATER = 1 LESS = -1 EQUAL = 0 | int __cmpxf2(long double, long double); SYM (__cmpxf2): lea GREATER,a1 | Unordered return value jra $_cmpxf | Go perform comparison SYM (__eqxf2): lea GREATER,a1 | Unordered return value jra $_cmpxf | Go perform comparison SYM (__nexf2): lea GREATER,a1 | Unordered return value jra $_cmpxf | Go perform comparison SYM (__gtxf2): lea LESS+LESS,a1 | Unordered return value jra $_cmpxf | Go perform comparison SYM (__gexf2): lea LESS+LESS,a1 | Unordered return value jra $_cmpxf | Go perform comparison SYM (__ltxf2): lea GREATER+GREATER,a1 | Unordered return value jra $_cmpxf | Go perform comparison SYM (__lexf2): lea GREATER+GREATER,a1 | Unordered return value jra $_cmpxf | Go perform comparison $_cmpxf: link a6,IMM (-24) moveml d2-d7,sp@ pea COMPARE+OP1_LDOUBLE+OP2_LDOUBLE+DST_INT movel a2,sp@- movel a6@(8),d4 | Get first operand movel a6@(12),d0 movel a6@(16),d1 movel a6@(20),d5 | Get second operand movel a6@(24),d2 movel a6@(28),d3 clrw d4 | Zap unused bits clrw d5 | First check if a and/or b are (+/-) zero and in that case clear | the sign bit. movel d4,d6 | copy signs into d6 (a) and a2 (b) bclr IMM (31),d4 | and clear signs in d0 and d2 movel d5,a2 bclr IMM (31),d5 subl d4,d6 | Extract the signs subl d5,a2 cmpl IMM (XF_INFINITY),d4 | check for a == NaN jeq Lcmpxf$a$big | if d4 == 0x7fff0000, a is NaN or INFINITY tstl d0 | Is a normalized? jpl Lcmpxf$a$mbz Lcmpxf$0: cmpl IMM (XF_INFINITY),d5 | check for b == NaN jeq Lcmpxf$b$big | if d5 == 0x7fff0000, b is NaN or INFINITY tstl d2 | Is b normalized? jpl Lcmpxf$b$mbz Lcmpxf$1: | Check the signs addl a2,d6 jpl 1f | If the signs are not equal check if a >= 0 movel a2,d6 jmi Lcmpxf$a$gt$b | if (a >= 0 && b < 0) => a > b jra Lcmpxf$b$gt$a | if (a < 0 && b >= 0) => a < b 1: | If the signs are equal check for < 0 movel a2,d6 jpl 1f | If both are negative exchange them movel d0,d7 movel d2,d0 movel d7,d2 movel d1,d7 movel d3,d1 movel d7,d3 movel d4,d7 movel d5,d4 movel d7,d5 1: | Now that they are positive we just compare them as longs (does this also | work for denormalized numbers?)(Yep). cmpl d4,d5 jgt Lcmpxf$b$gt$a | |b| > |a| jne Lcmpxf$a$gt$b | |b| < |a| | If we got here d4 == d5, so we compare d0 and d2. cmpl d0,d2 jhi Lcmpxf$b$gt$a | |b| > |a| jne Lcmpxf$a$gt$b | |b| < |a| | If we got here d0 == d2, so we compare d1 and d3. cmpl d1,d3 jhi Lcmpxf$b$gt$a | |b| > |a| jne Lcmpxf$a$gt$b | |b| < |a| | If we got here a == b. movel IMM (EQUAL),d0 jra Lcmpxf$ret Lcmpxf$a$gt$b: movel IMM (GREATER),d0 jra Lcmpxf$ret Lcmpxf$b$gt$a: movel IMM (LESS),d0 Lcmpxf$ret: movel sp@+,a2 moveml sp@(4),d2-d7 unlk a6 rts Lcmpxf$a$mbz: jne Lcmpxf$a$den | a isn't normal, go fix ' tstl d1 | Really zero? jne Lcmpxf$a$den movel d0,d6 | Make a +0 movel d0,d4 | & zap exponent jra Lcmpxf$0 Lcmpxf$a$den: tstl d4 | Is a a denormal? jeq Lcmpxf$0 | Yep, it's ok then ' 1: subl IMM (1), d4 | To denormal boundary yet? jeq 2f | Yep, go finish up addl d1,d1 | Shift another place addxl d0,d0 | Normalized yet? jpl 1b | Nope, do another step jra Lcmpxf$0 2: addl d1,d1 | Finish what we started addxl d0,d0 jra Lcmpxf$0 Lcmpxf$b$mbz: jne Lcmpxf$b$den | b isn't normal, go fix ' tstl d3 | Really zero? jne Lcmpxf$b$den movel d2,a2 | Make b +0 movel d2,d5 | & zap exponent jra Lcmpxf$1 Lcmpxf$b$den: tstl d5 | Is b a denormal? jeq Lcmpxf$1 | Yep, it's ok then ' 1: subl IMM (1), d5 | To denormal boundary yet? jeq 2f | Yep, go finish up addl d3,d3 | Shift another place addxl d2,d2 | Normalized yet? jpl 1b | Nope, do another step jra Lcmpxf$1 2: addl d3,d3 | Finish what we started addxl d2,d2 jra Lcmpxf$1 Lcmpxf$a$big: addl d0,d0 orl d0,d1 jeq Lcmpxf$0 | If equal it's +/-INFINITY ' tstl d0 | Is a a signaling NaN? jpl Lcmpxf$inv | Yep, go signal cmpl IMM (XF_INFINITY),d5 | Is b a signaling NaN? jne Lcmpxf$uno | Nope addl d2,d2 orl d2,d3 | INFINITY? jeq Lcmpxf$uno | Yep tstl d2 | Signaling NaN? jpl Lcmpxf$inv | Yep, go signal jra Lcmpxf$uno Lcmpxf$b$big: addl d2,d2 orl d2,d3 jeq Lcmpxf$1 | If equal it's +/-INFINITY ' tstl d2 | Is b a signaling NaN? jpl Lcmpxf$inv | Yep, report it Lcmpxf$uno: movel a1,d0 | Report no ordering for NaNs btst IMM (0),d0 | Error on unordered? jne Lcmpxf$ret Lcmpxf$inv: movel a1,d0 | Get unordered code asrl IMM (1),d0 | and dump invalid bit orl IMM (1),d0 movel IMM (INVALID_OPERATION),d7 movel sp@+,a2 jra $_exception_handler | and report exception |============================================================================= | Rounding and Exception Handling |============================================================================= | The non-zero value to be rounded and returned is in registers d1-d2-d3, | with the exponent in register d0, the sign in a1, and the operation type on | the top of the stack. The radix point is between bits 31 and 30 of d1, and | the exponent biased by 16383. | | Three exceptions are checked for: Overflow, Underflow, and Inexact. | | Overflow is checked after rounding. If detected, and trapped, the trap | handler is passed the rounded value, with the exponent adjusted by a bias, | 24576 in the case of EXTENED. The Inexact flag is set if the rounded value | is not exact. If not trapped, the value is set properly for the rounding | mode, and the Overflow and Inexact flags are set. The Inexact trap is | issued if enabled in this case. | | If the |value| <= 2^Emin (Emin = -16383 for EXTENDED) and the Underflow trap | is enabled, the value is normalize (if it isn't already), rounded, and sent ' | to the trap handler with the exponent bias adjusted. If the trap is not | enabled, the possibly denormalized value is rounded and, if inexact, the | Underflow flag is set. The Inexact flag is set if the result is not exact | in either case, and the Inexact trap taken if enabled, and the Underflow | trap is not. | We've got a normalized number, but need to denormalize it. ' denormal: lea SYM (_fpCCR),a2 btst IMM (1),a2@(TRAPE+3) | Check for Underflow trap enabled jne Lround$common | Yep, remain normalized cmpl IMM (-XF_SIGNIF_DIG-1),d0 | Really tiny? jls degen | Yep, just summarize movel IMM (32),d6 addl d6,d0 | Get left shift count jle 1f | Moving a whole word? subl d0,d6 | Get right shift count asrl IMM (1),d3 | Summarize low order bits subxl d4,d4 asrl IMM (1),d4 orl d4,d3 movel d2,d4 lsll d0,d4 orl d4,d3 lsrl d6,d2 movel d1,d7 lsll d0,d1 orl d1,d2 movel d7,d1 lsrl d6,d1 movel IMM (0),d0 jra Lround$common 1: jeq 1f | Shift even 32-bit words? addl d6,d0 | Get left shift count subl d0,d6 | and right shift count orl d2,d3 | Summarize lost bits asrl IMM (1),d3 subxl d4,d4 asrl IMM (1),d4 orl d4,d3 movel d1,d4 | Get least significant lsll d0,d4 orl d4,d3 | bits into d3 movel d1,d2 | and the rest into d2 lsrl d6,d2 2: movel IMM (0),d1 | Show zero in d1 movel d1,d0 | and denorm exponent jra Lround$common 1: asrl IMM (1),d3 | Summarize low order bits subxl d4,d4 asrl IMM (1),d4 orl d4,d3 orl d2,d3 | and shift others movel d1,d2 | down one place jra 2b degen: movel IMM (1),d3 | Show we're non-zero ' movel IMM (0),d2 | But all significance is gone jra 2b Lround$normal: | We first check for underflow. tstl d0 | Are we tiny? jlt denormal | Maybe, go check Lround$common: lea SYM (_fpCCR),a2 movel d3,d7 | Anything down here? jeq Lround$done | Exact! We're done! ' movel IMM (INEXACT_RESULT),d7 | Set Inexact flag movel a2@(ROUND),d6 | Get rounding mode in d6 jeq Lround$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Lround$to$minus jlt Lround$to$zero | jra Lround$to$plus Lround$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lround$done | Yep, truncate jra Lround$inc | Nope, increment Lround$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lround$done | Nope, truncate jra Lround$inc | Yep, increment Lround$to$nearest: | Now round: we do it as follows: after the shifting we can write the | fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2. | If delta < 1, do nothing. If delta > 1, add 1 to f. | If delta == 1, we make sure the rounded number will be even (odd?) | (after shifting). addl d3,d3 | Is delta < 1? jcc Lround$done | Yep, do nothing jne Lround$inc | Is delta == 1? movel IMM (1),d4 | If so round to even andl d2,d4 | bit 1 is the last significant bit addl d4,d2 jra 1f Lround$inc: addl IMM (1),d2 | Else add 1 1: movel IMM (0),d3 addxl d3,d1 | Rounding may have changed our exponent. If so, adjust. jcc Lround$done bset IMM (31),d1 | This only occurs when d2 == 0! addl IMM (1),d0 | Fixup exponent Lround$to$zero: Lround$done: | Now we're properly normalized and rounded, or, if the Underflow trap ' | is not enabled, maybe denormalized and rounded. If we're tiny ' | (|value| <= 2^Emin) the Tiny flag (proposed Underflow flag) is set. | If the rounded value is not exact, the Inexact flag is set. movel a2@(TRAPE),d5 | Get trap enables cmpl IMM (XF_MAX_EXP),d0 | Overflow? jge Ltrap$overflow | Nope, go check for Underflow tstl d0 | Check for Underflow jle Ltrap$underflow Ltrap$check: | Check for trapped exceptions. First we put the number back together. swap d0 | Put the exponent in place addl a1,d0 | Add in sign movel sp@+,a2 tstl d7 | Any exceptions? jne $_exception_handler | Yep, go handle Lx$return: movel d0,a0@ movel d1,a0@(4) movel d2,a0@(8) movel a0,d0 moveml sp@(4),d2-d7 | We're done ' unlk a6 rts Ltrap$underflow: jlt 2f | Really small (& trap enabled)? movel d1,d4 | |value| <= 2^Emin? addl d4,d4 jcc 3f | Yep, set Underflow jne Ltrap$check | Nope, return tstl d2 | |value| == 2^Emin? jne Ltrap$check | Nope 1: btst IMM (1),d5 | Underflow trap enabled? jne 2f | Yep, adjust the exponent btst IMM (0),d7 | Result Inexact? jeq Ltrap$check | Nope, won't set Underflow ' jra 3f 2: addl IMM (24576), d0 | Adjust exponent 3: addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Ltrap$check Ltrap$overflow: addl IMM (OVERFLOW),d7 | Set the Overflow flag btst IMM (2),d5 | Overflow trap enabled? jeq 1f | Nope, go adjust value subl IMM (24576),d0 | Adjust exponent jra Ltrap$check | and go trap 1: movel IMM (INEXACT_RESULT+OVERFLOW),d7 | Inexact too movel a2@(ROUND),d6 | Get rounding mode in d6 jeq Loflow$nearest | Yep, return infinity subl IMM (ROUND_TO_PLUS),d6 jhi Loflow$minus jlt Loflow$zero | jra Loflow$plus Loflow$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Loflow$nearest | Nope, go to infinity jra Loflow$zero | Yep, just go big Loflow$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Loflow$nearest | Yep, go to infinity jra Loflow$zero | Nope, just go big Loflow$nearest: movel IMM (0),d1 | Build an infinity value movel d1,d2 movel IMM (XF_MAX_EXP),d0 jra Ltrap$check Loflow$zero: movel IMM (0xffffffff),d1 | Build max value movel d1,d2 movel IMM (XF_MAX_EXP-1),d0 jra Ltrap$check Lx$a$snan: movel a6@(8),d0 | Return a movel a6@(12),d1 movel a6@(16),d2 jra Lx$snan Lx$b$snan: movel a6@(20),d0 | Return b movel a6@(24),d1 movel a6@(28),d2 Lx$snan: bset IMM (30),d1 | Change to Quiet NaN clrw d0 | Zap unused bits jra Lx$inv Lx$inop: | Invalid Operation. Return a quiet NaN and set the exception flags movel IMM (Q_NaN),d1 movel IMM (Q_NaN),d2 movel IMM (0xffff0000),d0 Lx$inv: movel IMM (INVALID_OPERATION),d7 movel sp@+,a2 jra $_exception_handler Lx$infty: | Return a properly signed INFINITY and set the exception flags movel IMM (XF_INFINITY),d0 movel IMM (0),d1 movel d1,d2 Lx$sign: addl a1,d0 | Sign bit movel sp@+,a2 jra Lx$return Lx$a: movel a6@(8),d0 | Return a movel a6@(12),d1 movel a6@(16),d2 clrw d0 movel sp@+,a2 jra Lx$return Lx$b: movel a6@(20),d0 | Return b movel a6@(24),d1 movel a6@(28),d2 clrw d0 movel sp@+,a2 jra Lx$return |============================================================================= | __extenddfxf2 |============================================================================= | long double __extenddfxf2 (double) SYM (__extenddfxf2): link a6,IMM (-24) moveml d2-d7,sp@ pea EXTENDDFXF+OP1_DOUBLE+DST_LDOUBLE movel a6@(8),d1 | Get operand movel a6@(12),d2 movel d1,a1 | Save sign away bclr IMM (31),d1 | & seperate from exponent subl d1,a1 movel IMM (DF_INFINITY),d0 andl d1,d0 | Zero or denormal? jeq 3f subl d0,d1 | Seperate from significand cmpl IMM (DF_INFINITY),d0 | +/-INFINITY or NaN? jeq 6f bset IMM (20),d1 | Set formerly hidden bit 1: asrl IMM (4),d0 | Put exponent in place addl IMM ((XF_BIAS-DF_BIAS)*0x10000),d0 | and adjust 8: movel IMM (11),d5 | Left shift count movel IMM (21),d6 | and right shift count lsll d5,d1 | Align upper bits movel d2,d4 lsrl d6,d4 orl d4,d1 | And carry over from lower bits lsll d5,d2 | Align lower bits 7: addl a1,d0 | Combine sign with exponent movel d0,a0@ | Save away result movel d1,a0@(4) movel d2,a0@(8) movel a0,d0 moveml sp@(4),d2-d7 unlk a6 rts | Here we have a denormalized float. Normalize it & correct the exponent 3: movel d1,d4 | Zero? orl d2,d4 jeq 7b | Yep, pass it on wider movel IMM (0x100000),d4 jra 5f 4: subl d4,d0 | Adjust exponent 5: addl d2,d2 | and shift significand addxl d1,d1 cmpl d4,d1 | Normalized yet? jcs 4b jra 1b | We've got a +/-INFINITY or a NaN ' 6: movel d1,d4 orl d2,d4 jne 1f | NaN? movel IMM (XF_INFINITY),d0 | Nope, return long double +/-INFINITY jra 7b |We've a NaN here. Is it a signaling one? ' 1: movel IMM (XF_INFINITY),d0 bset IMM (19),d1 jne 8b | Nope, just pass it on movel IMM (11),d5 | Left shift count movel IMM (21),d6 | and right shift count lsll d5,d1 | Align upper bits movel d2,d4 lsrl d6,d4 orl d4,d1 | And carry over from lower bits lsll d5,d2 | Align lower bits addl a1,d0 | Combine sign with exponent movel IMM (INVALID_OPERATION),d7 jra $_exception_handler |============================================================================= | __truncxfdf2 |============================================================================= | double __truncxfdf2 (long double) SYM (__truncxfdf2): link a6,IMM (-24) moveml d2-d7,sp@ pea TRUNCXFDF+OP1_LDOUBLE+DST_DOUBLE lea SYM (_fpCCR),a0 movel a6@(8),d4 | Get exponent & sign movel a6@(12),d0 movel a6@(16),d1 movel d4,a1 | Isolate sign bclr IMM (31),d4 subl d4,a1 clrw d4 swap d4 | Put exponent in place cmpl IMM (XF_MAX_EXP),d4 | Special case? jeq Lxfdf$nf | Yep, go handle tstl d0 jpl Lxfdf$mbz | Unnormalized/zero? | We've got a canonical finite number here. Check for underflow/overflow. ' Lxfdf$0: cmpl IMM (XF_BIAS-DF_BIAS),d4 jle Lxfdf$under | Have underflow, go handle | Round & check for overflow Lxfdf$1: movel IMM (0x7ff),d7 andl d1,d7 | Anything down here? jeq Lxfdf$to$zero movel d7,d3 movel IMM (INEXACT_RESULT),d7 | Set Inexact flag movel IMM (0x800),d2 | Increment constant movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lxfdf$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Lxfdf$to$minus jlt Lxfdf$to$zero | jra Lxfdf$to$plus Lxfdf$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lxfdf$zero | Yep, truncate jra Lxfdf$inc | Nope, increment Lxfdf$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lxfdf$zero | Nope, truncate jra Lxfdf$inc | Yep, increment Lxfdf$to$nearest: cmpl IMM (0x400),d3 | Which way is nearest? jlt Lxfdf$to$zero jne Lxfdf$inc andl d1,d2 | Half way; round to low bit zero Lxfdf$inc: addl d2,d1 | Increment low bits movel IMM (0),d2 addxl d2,d0 | and carry to high bits jcc Lxfdf$to$zero | Need to renormalize? bset IMM (31),d0 | Yep, we're power of 2 then! ' addl IMM (1),d4 Lxfdf$to$zero: movel a0@(TRAPE),d5 | Get trap enables addl IMM (DF_RBIAS-XF_BIAS),d4 | Adjust exponent jle Lxfdf$underflow | Underflow? cmpl IMM (DF_MAX_EXP-1),d4 | Overflow? jge Lxfdf$overflow | We're in range here. The only trap we can have is inexact. ' Lxfdf$check: movel IMM (11),d2 | Realign significand lsrl d2,d1 | to double format movel d0,d3 lsrl d2,d0 movel IMM (21),d2 lsll d2,d3 orl d3,d1 swap d4 | Put the exponent in place lsll IMM (4),d4 addl d4,d0 | and add one to exponent | Check for trapped exceptions. First we put the number back together. addl a1,d0 | Add in sign tstl d7 | Any exceptions? jne $_exception_handler | Yep, go handle Lxfdf$ret: moveml sp@(4),d2-d7 | We're done; ' unlk a6 | return rts Lxfdf$underflow: jlt 1f | Really small (& trap enabled)? movel d0,d2 | |value| < 2^Emin? jpl 1f | Yep, set Underflow addl d2,d2 orl d1,d2 | |value| == 2^Emin? jne Lxfdf$check | Nope 1: btst IMM (1),d5 | Underflow trap enabled? jne 2f | Yep, return rounded DOUBLE btst IMM (0),d7 | Result Inexact? jeq Lxfdf$check | Nope, won't set Underflow ' addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Lxfdf$check 2: addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Lxfdf$dtrap Lxfdf$overflow: addl IMM (OVERFLOW),d7 | Set the Overflow flag btst IMM (2),d5 | Overflow trap enabled? jne Lxfdf$dtrap | Yep, return rounded DOUBLE movel IMM (INEXACT_RESULT+OVERFLOW),d7 | Inexact too movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lxfdf$nearest | Yep, return infinity subl IMM (ROUND_TO_PLUS),d6 jgt Lxfdf$minus jlt Lxfdf$zero | jra Lxfdf$plus Lxfdf$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lxfdf$nearest | Nope, go to infinity jra Lxfdf$zero | Yep, just go big Lxfdf$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lxfdf$nearest | Yep, go to infinity jra Lxfdf$zero | Nope, just go big Lxfdf$nearest: movel IMM (0),d0 | Build an infinity value movel d0,d1 movel IMM (DF_MAX_EXP),d4 jra Lxfdf$check Lxfdf$zero: movel IMM (-1),d0 | Build max value movel IMM (-1),d1 movel IMM (DF_MAX_EXP-2),d4 jra Lxfdf$check | We've got an Underflow or Overflow trap. Pass a rounded LONG DOUBLE ' | to the trap handler. It had better convert it to a DOUBLE as that's ' | what the program wants, not the first two words of a LONG DOUBLE! | (Come back some day & fix for case when rounding causes overflow | of long double format!!!!!) Lxfdf$dtrap: movel IMM (0xfffff800),d2 | Strip off low bits andl d1,d2 | & put in correct registers movel d0,d1 subl IMM (DF_RBIAS-XF_BIAS),d4 | Restore the exponent swap d4 movel d4,d0 addl a1,d0 | Combine with sign jra $_exception_handler | Significand is either zero or unnormalized Lxfdf$mbz: movel d1,d2 orl d0,d2 | Zero? jeq Lxfdf$sign | Shift left until normalized or at the denormal boundary tstl d4 jeq Lxfdf$under 1: subl IMM (1),d4 jeq 3f addl d1,d1 addxl d0,d0 jpl 1b jra Lxfdf$0 3: addl d1,d1 addxl d0,d0 | We have an underflow. If it's not too far out, and the underflow ' | trap isn't enabled, we may return a DOUBLE denormal ' Lxfdf$under: btst IMM (1),a0@(TRAPE+3) | Underflow trap enabled? jne Lxfdf$1 | Yep, go round whatever | we've got & trap ' addl IMM (DF_RBIAS-XF_BIAS),d4 | Adjust the bias cmpl IMM (-DF_SIGNIF_DIG-1),d4 | Anything going to be left? jlt 4f | Nope, just summarize movel IMM (31),d7 tstl d4 | Denormal yet? jeq 2f 1: lsrl IMM (1),d1 subxl d2,d2 negl d2 orl d2,d1 | Retain sticky bit lsrl IMM (1),d0 subxl d2,d2 lsll d7,d2 orl d2,d1 addl IMM (1),d4 jlt 1b 2: subl IMM (DF_RBIAS-XF_BIAS),d4 | Adjust the bias jra Lxfdf$1 | Go round the denormal 4: movel IMM (1),d1 movel IMM (0),d0 movel IMM (XF_BIAS-DF_RBIAS),d4 jra Lxfdf$1 | +/-INFINITY or NaN Lxfdf$nf: bclr IMM (31),d0 | Dump the integer bit movel d0,d2 orl d1,d2 | +/-INFINITY? jne Lxfdf$nan | Nope, go handle NaN movel IMM (DF_INFINITY),d0 | Make a DOUBLE +/-INFINITY Lxfdf$sign: addl a1,d0 | Add in sign jra Lxfdf$ret Lxfdf$nan: movel IMM (11),d2 lsrl d2,d1 | Preserve pattern movel d0,d3 lsrl d2,d0 movel IMM (21),d2 lsll d2,d3 orl d3,d1 orl IMM (DF_INFINITY),d0 bset IMM (19),d0 | Signaling NaN? jne Lxfdf$sign movel IMM (INVALID_OPERATION),d7 jra $_exception_handler |============================================================================= | __extendsfxf2 |============================================================================= | long double __extendsfxf2 (float) SYM (__extendsfxf2): link a6,IMM (-24) moveml d2-d7,sp@ pea EXTENDSFXF+OP1_FLOAT+DST_LDOUBLE movel a6@(8),d1 | Get operand movel IMM (0),d2 | Zap lower bits movel d1,a1 | Save sign away bclr IMM (31),d1 | & seperate from exponent subl d1,a1 movel IMM (SF_INFINITY),d0 andl d1,d0 | Zero or denormal? jeq 3f subl d0,d1 | Seperate from significand cmpl IMM (SF_INFINITY),d0 | +/-INFINITY or NaN? jeq 6f bset IMM (23),d1 | Set formerly hidden bit 1: asrl IMM (7),d0 | Put exponent in place addl IMM ((XF_BIAS-SF_BIAS)*0x10000),d0 | and adjust 8: lsll IMM (8),d1 | Align upper bits 7: addl a1,d0 | Combine sign with exponent movel d0,a0@ | Save away result movel d1,a0@(4) movel d2,a0@(8) movel a0,d0 moveml sp@(4),d2-d7 unlk a6 rts | Here we have a denormalized float. Normalize it & correct the exponent 3: tstl d1 | Zero? jeq 7b | Yep, pass it on wider movel IMM (0x800000),d4 jra 5f 4: subl d4,d0 | Adjust exponent 5: addl d1,d1 | and shift significand cmpl d4,d1 | Normalized yet? jcs 4b jra 1b | We've got a +/-INFINITY or a NaN ' 6: tstl d1 jne 1f | NaN? movel IMM (XF_INFINITY),d0 | Nope, return long double +/-INFINITY jra 7b |We've a NaN here. Is it a signaling one? ' 1: movel IMM (XF_INFINITY),d0 bset IMM (22),d1 jne 8b | Nope, just pass it on lsll IMM (8),d1 | Align upper bits addl a1,d0 | Combine sign with exponent movel IMM (INVALID_OPERATION),d7 jra $_exception_handler |============================================================================= | __truncxfsf2 |============================================================================= | float __truncxfsf2 (long double) SYM (__truncxfsf2): link a6,IMM (-24) moveml d2-d7,sp@ pea TRUNCXFSF+OP1_LDOUBLE+DST_FLOAT lea SYM (_fpCCR),a0 movel a6@(8),d4 | Get exponent & sign movel a6@(12),d0 movel a6@(16),d1 movel d4,a1 | Isolate sign bclr IMM (31),d4 subl d4,a1 clrw d4 swap d4 | Put exponent in place cmpl IMM (XF_MAX_EXP),d4 | Special case? jeq Lxfsf$nf | Yep, go handle tstl d0 jpl Lxfsf$mbz | Unnormalized/zero? | We've got a canonical finite number here. Check for underflow/overflow. ' Lxfsf$0: cmpl IMM (XF_BIAS-SF_RBIAS),d4 jle Lxfsf$under | Have underflow, go handle | Round & check for overflow Lxfsf$1: tstl d1 jeq 1f orl IMM (1),d0 | Set sticky bit 1: movel IMM (0xff),d7 andl d0,d7 | Anything down here? jeq Lxfsf$to$zero movel d7,d3 movel IMM (INEXACT_RESULT),d7 | Set Inexact flag movel IMM (0x100),d2 | Increment constant movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lxfsf$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Lxfsf$to$minus jlt Lxfsf$to$zero | jra Lxfsf$to$plus Lxfsf$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lxfsf$zero | Yep, truncate jra Lxfsf$inc | Nope, increment Lxfsf$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lxfsf$zero | Nope, truncate jra Lxfsf$inc | Yep, increment Lxfsf$to$nearest: cmpl IMM (0x80),d3 | Which way is nearest? jlt Lxfsf$to$zero jne Lxfsf$inc andl d0,d2 | Half way; round to low bit zero Lxfsf$inc: addl d2,d0 | Increment low bits jcc Lxfsf$to$zero | Need to renormalize? bset IMM (31),d0 | Yep, we're power of 2 then! ' addl IMM (1),d4 Lxfsf$to$zero: movel a0@(TRAPE),d5 | Get trap enables addl IMM (SF_RBIAS-XF_BIAS),d4 | Adjust exponent jle Lxfsf$underflow | Underflow? cmpl IMM (SF_MAX_EXP-1),d4 | Overflow? jge Lxfsf$overflow | We're in range here. The only trap we can have is inexact. ' Lxfsf$check: lsrl IMM (8),d0 | Realign significand to float swap d4 | Put the exponent in place lsll IMM (7),d4 addl d4,d0 | and add one to exponent | Check for trapped exceptions. First we put the number back together. addl a1,d0 | Add in sign tstl d7 | Any exceptions? jne $_exception_handler | Yep, go handle Lxfsf$ret: moveml sp@(4),d2-d7 | We're done; ' unlk a6 | return rts Lxfsf$underflow: jlt 1f | Really small (& trap enabled)? movel d0,d2 | |value| < 2^Emin? jpl 1f | Yep, set Underflow addl d2,d2 orl d1,d2 | |value| == 2^Emin? jne Lxfsf$check | Nope 1: btst IMM (1),d5 | Underflow trap enabled? jne 2f | Yep, return rounded DOUBLE btst IMM (0),d7 | Result Inexact? jeq Lxfsf$check | Nope, won't set Underflow ' addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Lxfsf$check 2: addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Lxfsf$dtrap Lxfsf$overflow: addl IMM (OVERFLOW),d7 | Set the Overflow flag btst IMM (2),d5 | Overflow trap enabled? jne Lxfsf$dtrap | Yep, return rounded DOUBLE movel IMM (INEXACT_RESULT+OVERFLOW),d7 | Inexact too movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lxfsf$nearest | Yep, return infinity subl IMM (ROUND_TO_PLUS),d6 jgt Lxfsf$minus jlt Lxfsf$zero | jra Lxfsf$plus Lxfsf$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lxfsf$nearest | Nope, go to infinity jra Lxfsf$zero | Yep, just go big Lxfsf$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lxfsf$nearest | Yep, go to infinity jra Lxfsf$zero | Nope, just go big Lxfsf$nearest: movel IMM (0),d0 | Build an infinity value movel IMM (SF_MAX_EXP),d4 jra Lxfsf$check Lxfsf$zero: movel IMM (-1),d0 | Build max value movel IMM (SF_MAX_EXP-2),d4 jra Lxfsf$check | We've got an Underflow or Overflow trap. Pass a rounded LONG DOUBLE ' | to the trap handler. It had better convert it to a FLOAT as that's ' | what the program wants, not the first word of a LONG DOUBLE! | (Come back some day & fix for case when rounding causes overflow | of long double format!!!!!) Lxfsf$dtrap: movel IMM (0xffffff00),d1 | Strip off low bits andl d0,d1 | & put in correct registers clrl d2 subl IMM (SF_RBIAS-XF_BIAS),d4 | Restore the exponent swap d4 movel d4,d0 addl a1,d0 | Combine with sign jra $_exception_handler | Significand is either zero or unnormalized Lxfsf$mbz: movel d1,d2 orl d0,d2 | Zero? jeq Lxfsf$sign btst IMM (1),a0@(TRAPE+3) | Underflow trap enabled? jeq 4f | Nope, just summarize | Shift left until normalized or at the denormal boundary tstl d4 jeq Lxfsf$under 1: subl IMM (1),d4 jeq 3f addl d1,d1 addxl d0,d0 jpl 1b jra Lxfsf$0 3: addl d1,d1 addxl d0,d0 | We have an underflow. If it's not too far out, and the underflow ' | trap isn't enabled, we may return a FLOAT denormal ' Lxfsf$under: btst IMM (1),a0@(TRAPE+3) | Underflow trap enabled? jne Lxfdf$1 | Yep, go round whatever | we've got & trap ' addl IMM (SF_RBIAS-XF_BIAS),d4 | Adjust the bias cmpl IMM (-DF_SIGNIF_DIG-1),d4 | Anything going to be left? jlt 4f | Nope, just summarize tstl d4 | Denormal yet? jeq 2f 1: lsrl IMM (1),d0 subxl d2,d2 negl d2 orl d2,d0 | Retain sticky bit addl IMM (1),d4 jlt 1b 2: subl IMM (SF_RBIAS-XF_BIAS),d4 | Adjust the bias jra Lxfsf$1 | Go round the denormal 4: movel IMM (1),d1 movel IMM (0),d0 movel IMM (XF_BIAS-SF_RBIAS),d4 jra Lxfsf$1 | +/-INFINITY or NaN Lxfsf$nf: bclr IMM (31),d0 | Dump the integer bit movel d0,d2 orl d1,d2 | +/-INFINITY? jne Lxfsf$nan | Nope, go handle NaN movel IMM (SF_INFINITY),d0 | Make a FLOAT +/-INFINITY Lxfsf$sign: addl a1,d0 | Add in sign jra Lxfsf$ret Lxfsf$nan: lsrl IMM (8),d0 | Preserve pattern orl IMM (SF_INFINITY),d0 bset IMM (22),d0 | Signaling NaN? jne Lxfsf$sign movel IMM (INVALID_OPERATION),d7 jra $_exception_handler |============================================================================= | __floatsixf |============================================================================= | long double __floatsixf (int) SYM (__floatsixf): movel IMM (0),d0 | Assume plus movel sp@(4),d1 | Get the operand jeq Lsixf$zer | Go return zero jpl 1f | Correct assumption? negl d1 | Nope, take abs movel IMM (0x8000),d0 | & show minus 1: | Now we have a non-zero number we need to normalize. First we shift | four places at a time to get there faster. lea 0x10000000,a1 addl IMM (XF_BIAS+31),d0 cmpl a1,d1 jcc 2f | Already in range? 1: subl IMM (4),d0 | Adjust exponent lsll IMM (4),d1 | and shift significand cmpl a1,d1 jcs 1b 2: | We're now within just a few places of being normalized. Do one ' | shift at a time tstl d1 | There yet? jmi 2f | Yep, go pack up 1: subl IMM (1),d0 | Adjust exponent lsll IMM (1),d1 | and shift significand jpl 1b 2: | The significand is normalized in d1, the exponent is in d0. swap d0 | Put the exponent in its place Lsixf$ret: movel d0,a0@ movel d1,a0@(4) clrl a0@(8) rts Lsixf$zer: movel d1,d0 jra Lsixf$ret |============================================================================= | __fixxfsi |============================================================================= | int __fixxfsi (double) | Note that gcc wants round-to-zero mode here SYM (__fixxfsi): link a6,IMM (-24) moveml d2-d7,sp@ pea FIXXFSI+OP1_LDOUBLE+DST_INT lea SYM (_fpCCR),a0 movel a6@(8),d4 | Get operand clrw d4 | Zap unused bits movel a6@(12),d0 movel a6@(16),d1 movel d4,a1 | Isolate sign bclr IMM (31),d4 subl d4,a1 swap d4 cmpl IMM (XF_MAX_EXP),d4 | special case? jeq Lxfsi$inop | Yep, go handle tstl d0 jpl Lxfsi$mbz | Zero or unnormalized? Lxfsi$1: subl IMM (XF_BIAS-1),d4 | Less than 1/2? jlt Lxfsi$small | Yep, go show small movel IMM (32),d5 cmpl d5,d4 | >= 2**32? jgt Lxfsi$max | Yep, illegal subl d4,d5 | Find how far to shift jeq Lxfsi$round | None at all? movel d1,d3 lsrl d5,d1 | Align lower bits lsll d4,d3 | & set sticky bit jeq 1f orl IMM (1),d1 1: movel d0,d3 lsll d4,d3 orl d3,d1 | Carry over from upper lsrl d5,d0 | Align upper bits | Round & check for overflow Lxfsi$round: movel d1,d7 | Anything down here? jeq Lxfsi$to$zero movel IMM (INEXACT_RESULT),d7 | Set Inexact flag #ifdef FIX_ROUND_MODE movel IMM (1),d2 | Increment constant movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lxfsi$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Lxfsi$to$minus jlt Lxfsi$to$zero | jra Lxfsi$to$plus Lxfsi$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lxfsi$to$zero | Yep, truncate jra Lxfsi$inc | Nope, increment Lxfsi$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lxfsi$to$zero | Nope, truncate jra Lxfsi$inc | Yep, increment Lxfsi$to$nearest: addl d1,d1 | Which way is nearest? jcc Lxfsi$to$zero jne Lxfsi$inc andl d0,d2 Lxfsi$inc: addl d2,d0 jcs Lxfsi$max | Overflow! #endif /* FIX_ROUND_MODE */ Lxfsi$to$zero: tstl d0 | Safe? jmi Lxfsi$mimnax | Maybe minimum? #ifndef __mcf5200__ cmpa #0,a1 | Minus? #else tstl a1 | Minus? #endif jpl Lxfsi$eck negl d0 | Return minus Lxfsi$eck: tstl d7 | Exact? jne $_exception_handler | Nope, check for traps Lxfsi$ret: moveml sp@(4),d2-d7 | Restore regs unlk a6 | and return rts Lxfsi$mimnax: #ifndef __mcf5200__ cmpa #0,a1 | Minus? #else tstl a1 | Minus? #endif jpl Lxfsi$max | Nope, Overflow addl d0,d0 | Min value? jne Lxfsi$max | Nope, Overflow movel a1,d0 | Yep, get it back jra Lxfsi$eck | and check for exceptions | Have zero or an unnormalized number Lxfsi$mbz: movel d0,d2 orl d1,d2 | Zero? jeq Lxfsi$ret | Yep, return it | Shift left until normalized or at the denormal boundary tstl d4 jeq Lxfsi$small 1: subl IMM (1),d4 jeq Lxfsi$small addl d1,d1 addxl d0,d0 jpl 1b jra Lxfsi$1 | Summarize and round for non-zero |values| < 1/2 Lxfsi$small: movel IMM (0),d0 | Show very small movel IMM (1),d1 jra Lxfsi$round | Number won't fit in int, Invalid Operation Exception ' Lxfsi$max: movel a1,d0 | Get sign jmi Lxfsi$inv | If minus, it's min ' movel IMM (0x7fffffff),d0 | Get max for plus jra Lxfsi$inv Lxfsi$inop: movel IMM (0),d0 | Nothing good to return Lxfsi$inv: movel IMM (INVALID_OPERATION),d7 jra $_exception_handler #endif /* L_extended */ #ifdef L_cvt | Entry points: .globl SYM (__floatsidf) .globl SYM (__floatsisf) .globl SYM (__extendsfdf2) .globl SYM (__truncdfsf2) .globl SYM (__fixdfsi) .globl SYM (__fixsfsi) .text .even |============================================================================= |============================================================================= | conversion routines |============================================================================= |============================================================================= |============================================================================= | __floatsidf |============================================================================= | double __floatsidf (int) SYM (__floatsidf): subl a1,a1 | Assume plus movel sp@(4),d0 | Get the operand jeq Lsidf$zer | Go return zero jpl 1f | Correct assumption? negl d0 | Nope, take abs lea 0x80000000,a1 | & show minus 1: | Now we have a non-zero number we need to normalize. First we shift | four places at a time to get there faster. lea 0x10000000,a0 movel IMM (DF_RBIAS+31),d1 cmpl a0,d0 jcc 2f | Already in range? 1: subl IMM (4),d1 | Adjust exponent lsll IMM (4),d0 | and shift significand cmpl a0,d0 jcs 1b 2: | We're now within just a few places of being normalized. Do one ' | shift at a time tstl d0 | There yet? jmi 2f | Yep, go pack up 1: subl IMM (1),d1 | Adjust exponent lsll IMM (1),d0 | and shift significand jpl 1b 2: | The significand is normalized in d0, the exponent is in d1, and the | sign in a1. We just need to shuffle them around and combine them. lsll IMM (4),d1 | Put the exponent in its place swap d1 addl d1,a1 | and combine with the sign movel d0,d1 lsrl IMM (8),d0 lsrl IMM (3),d0 swap d1 clrw d1 lsll IMM (5),d1 addl a1,d0 rts Lsidf$zer: movel d0,d1 rts |============================================================================= | __floatsisf |============================================================================= | float __floatsisf (int) SYM (__floatsisf): link a6,IMM (-24) moveml d2-d7,sp@ pea FLOATSISF+OP1_INT+DST_FLOAT subl a1,a1 | Assume plus movel a6@(8),d0 | Get the operand jeq Lsisf$ret | Go return zero jpl 1f | Correct assumption? negl d0 | Nope, take abs lea 0x80000000,a1 | & show minus 1: | Now we have a non-zero number we need to normalize. First we shift | four places at a time to get there faster. lea 0x10000000,a0 movel IMM (SF_RBIAS+31),d4 cmpl a0,d0 jcc 2f | Already in range? 1: subl IMM (4),d4 | Adjust exponent lsll IMM (4),d0 | and shift significand cmpl a0,d0 jcs 1b 2: | We're now within just a few places of being normalized. Do one ' | shift at a time tstl d0 | There yet? jmi 2f | Yep, go pack up 1: subl IMM (1),d4 | Adjust exponent lsll IMM (1),d0 | and shift significand jpl 1b 2: | The significand is normalized in d0, the exponent is in d4, and the | sign in a1. Round to fit in a float significand lea SYM (_fpCCR),a0 movel IMM (0xff),d7 andl d0,d7 | Anything down here? jeq Lsisf$to$zero movel d7,d3 movel IMM (INEXACT_RESULT),d7 | Set Inexact flag movel IMM (0x100),d2 | Increment constant movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Lsisf$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Lsisf$to$minus jlt Lsisf$to$zero | jra Lsisf$to$plus Lsisf$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Lsisf$to$zero | Yep, truncate jra Lsisf$inc | Nope, increment Lsisf$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Lsisf$to$zero | Nope, truncate jra Lsisf$inc | Yep, increment Lsisf$to$nearest: cmpl IMM (0x80),d3 | Which way is nearest? jcs Lsisf$to$zero jne Lsisf$inc andl d0,d2 Lsisf$inc: addl d2,d0 jcc Lsisf$to$zero bset IMM (31),d0 addl IMM (1),d4 | The only trap we can have is inexact. Lsisf$to$zero: lsrl IMM (8),d0 swap d4 | Put the exponent in place lsll IMM (7),d4 addl d4,d0 | and add one to exponent | Check for trapped exceptions. First we put the number back together. addl a1,d0 | Add in sign tstl d7 | Any exceptions? jne $_exception_handler | Yep, go handle Lsisf$ret: moveml sp@(4),d2-d7 | We're done ' unlk a6 rts |============================================================================= | __extendsfdf2 |============================================================================= | double __extendsfdf2 (float) SYM (__extendsfdf2): movel sp@(4),d0 | Get sign & exponent movel d0,a0 | Save sign away bclr IMM (31),d0 | & seperate from eponent subl d0,a0 movel IMM (SF_INFINITY),d1 andl d0,d1 | Zero or denormal? jeq 3f cmpl IMM (SF_INFINITY),d1 | +/-INFINITY or NaN? jeq 6f subl d1,d0 | Seperate from significand 1: asrl IMM (3),d1 | Put exponent in place addl IMM ((DF_BIAS-SF_BIAS)*0x100000),d1 | and adjust 8: addl d1,a0 | Save away with sign movel IMM (7),d1 andl d0,d1 | Put low 3 bits in low word swap d1 | and shift into place lsll IMM (8),d1 lsll IMM (5),d1 lsrl IMM (3),d0 | Shift upper bits into place 7: addl a0,d0 | and combine with sign/exponent rts | Here we have a denormalized float. Normalize it & correct the exponent 3: tstl d0 | Zero? jeq 7b lea 0x800000,a1 | We leave the hidden bit set, and decrement the exponent once too many | times. They compensate for each other, and it takes fewer instructions | this way. 4: subl a1,d1 | Adjust exponent lsll IMM (1),d0 | and shift significand cmpl a1,d0 | Normalized yet? jcs 4b jra 1b | We've got a +/-INFINITY or a NaN ' 6: subl d1,d0 | Seperate from significand jne 1f | NaN? movel IMM (DF_INFINITY),d0 | Nope, return double +/-INFINITY movel IMM (0),d1 jra 7b |We've a NaN here. Is it a signaling one? ' 1: movel IMM (DF_INFINITY),d1 bset IMM (22),d0 jne 8b | Nope, just pass it on jsr 8b | Fix up to QNaN link a6,IMM (-24) moveml d2-d7,sp@ pea EXTENDSFDF+OP1_FLOAT+DST_DOUBLE movel IMM (INVALID_OPERATION),d7 jra $_exception_handler |============================================================================= | __truncdfsf2 |============================================================================= | float __truncdfsf2 (double) SYM (__truncdfsf2): link a6,IMM (-24) moveml d2-d7,sp@ pea TRUNCDFSF+OP1_DOUBLE+DST_FLOAT lea SYM (_fpCCR),a0 movel a6@(8),d0 | Get operand movel a6@(12),d1 movel d0,a1 | Isolate sign bclr IMM (31),d0 subl d0,a1 movel IMM (0xfff00000),d3 | Exponent mask movel d0,d4 andl d3,d4 addl d3,d4 | Adjust exponent subl d4,d0 | Fix up significand asrl IMM (4),d4 | Put exponent in place swap d4 cmpl IMM (0x7fe),d4 | special case? jcc Ldfsf$spec | Yep, go handle | We've got a normalized number here. Check for underflow/overflow. ' cmpl IMM (DF_RBIAS-SF_RBIAS),d4 jle Ldfsf$under | Have underflow, go handle | Round & check for overflow Ldfsf$1: movel IMM(0x1fffffff),d7 andl d1,d7 | Anything down here? jeq Ldfsf$to$zero movel d7,d3 movel IMM (INEXACT_RESULT),d7 | Set Inexact flag movel IMM (0x20000000),d2 | Increment constant movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Ldfsf$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Ldfsf$to$minus jlt Ldfsf$to$zero | jra Ldfsf$to$plus Ldfsf$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Ldfsf$zero | Yep, truncate jra Ldfsf$inc | Nope, increment Ldfsf$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Ldfsf$zero | Nope, truncate jra Ldfsf$inc | Yep, increment Ldfsf$to$nearest: cmpl IMM (0x10000000),d3 | Which way is nearest? jcs Ldfsf$to$zero jne Ldfsf$inc andl d1,d2 Ldfsf$inc: addl d2,d1 movel IMM (0),d2 addxl d2,d0 btst IMM (21),d0 jeq Ldfsf$to$zero lsrl IMM (1),d0 addl IMM (1),d4 Ldfsf$to$zero: movel a0@(TRAPE),d5 | Get trap enables addl IMM (SF_RBIAS-DF_RBIAS),d4 | Adjust exponent jle Ldfsf$underflow | Underflow? cmpl IMM (SF_MAX_EXP-1),d4 | Overflow? jge Ldfsf$overflow | Nope, go check for Underflow | We're in range here. The only trap we can have is inexact. ' Ldfsf$check: addl d1,d1 | Realign significand addxl d0,d0 addl d1,d1 addxl d0,d0 addl d1,d1 addxl d0,d0 swap d4 | Put the exponent in place lsll IMM (7),d4 addl d4,d0 | and add one to exponent | Check for trapped exceptions. First we put the number back together. addl a1,d0 | Add in sign tstl d7 | Any exceptions? jne $_exception_handler | Yep, go handle Ldfsf$ret: moveml sp@(4),d2-d7 | We're done ' unlk a6 rts Ldfsf$underflow: jlt 2f | Really small (& trap enabled)? cmpl IMM (0x100000),d0 | |value| <= 2^Emin? jgt Ldfsf$check | Nope, return jlt 1f | Yep, set Underflow andl IMM (0xe0000000),d1 | |value| == 2^Emin? jne Ldfsf$check | Nope 1: btst IMM (1),d5 | Underflow trap enabled? jne 2f | Yep, return rounded DOUBLE btst IMM (0),d7 | Result Inexact? jeq Ldfsf$check | Nope, won't set Underflow ' addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Ldfsf$check 2: addl IMM (UNDERFLOW),d7 | Set Underflow flag jra Ldfsf$dtrap Ldfsf$overflow: addl IMM (OVERFLOW),d7 | Set the Overflow flag btst IMM (2),d5 | Overflow trap enabled? jne Ldfsf$dtrap | Yep, return rounded DOUBLE movel IMM (INEXACT_RESULT+OVERFLOW),d7 | Inexact too movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Ldfsf$nearest | Yep, return infinity subl IMM (ROUND_TO_PLUS),d6 jhi Ldfsf$minus jlt Ldfsf$zero | jra Ldfsf$plus Ldfsf$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Ldfsf$nearest | Nope, go to infinity jra Ldfsf$zero | Yep, just go big Ldfsf$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Ldfsf$nearest | Yep, go to infinity jra Ldfsf$zero | Nope, just go big Ldfsf$nearest: movel IMM (0),d0 | Build an infinity value movel d0,d1 movel IMM (SF_MAX_EXP),d4 jra Ldfsf$check Ldfsf$zero: movel IMM (0xffffff),d0 | Build max value movel IMM (-1),d1 movel IMM (SF_MAX_EXP-2),d4 jra Ldfsf$check | We've got an Underflow or Overflow trap. Pass a rounded DOUBLE ' | to the trap handler. Ldfsf$dtrap: andl IMM (0xe0000000),d1 | Strip off low bits subl IMM (SF_RBIAS-DF_RBIAS),d4 | Restore the exponent swap d4 lsll IMM (4),d4 addl d4,d0 | Combine with significand addl a1,d0 | Add in sign movel IMM (TRUNCDFSF+OP1_DOUBLE+DST_DOUBLE),sp@ jra $_exception_handler | Zero, denormal, +/- INFINITY or NaN Ldfsf$spec: jeq Ldfsf$nf | Go handle big stuff movel IMM (0),d4 | Fix up the exponent bclr IMM (20),d0 | and significand movel d1,d2 orl d0,d2 | Have a zero? jeq Ldfsf$sign | Yep, return it | We have an underflow. If it's not too far out, and the underflow ' | trap isn't enabled, we may return a FLOAT denormal ' Ldfsf$under: btst IMM (1),a0@(TRAPE+3) | Underflow trap enabled? jne Ldfsf$1 | Yep, go round whatever | we've got & trap ' addl IMM (SF_RBIAS-DF_RBIAS),d4 | Adjust the bias cmpl IMM (-SF_SIGNIF_DIG-1),d4 | Anything going to be left? jlt 4f | Nope, just summarize movel IMM (31),d7 tstl d4 | Denormal yet? jeq 2f 1: lsrl IMM (1),d1 subxl d2,d2 negl d2 orl d2,d1 lsrl IMM (1),d0 subxl d2,d2 lsll d7,d2 orl d2,d1 addl IMM (1),d4 jlt 1b 2: subl IMM (SF_RBIAS-DF_RBIAS),d4 | Adjust the bias jra Ldfsf$1 | Go round the denormal 4: movel IMM (1),d1 movel IMM (0),d0 movel IMM (DF_RBIAS-SF_RBIAS),d4 jra Ldfsf$1 | +/-INFINITY or NaN Ldfsf$nf: bclr IMM (20),d0 | Zap the hidden bit movel d1,d2 orl d0,d2 | +/-INFINITY? jne Ldfsf$nan | Nope, go handle NaN movel IMM (SF_INFINITY),d0 | Make a FLOAT +/-INFINITY Ldfsf$sign: addl a1,d0 | Add in sign jra Ldfsf$ret Ldfsf$nan: addl d1,d1 | Preserve pattern addxl d0,d0 addl d1,d1 addxl d0,d0 addl d1,d1 addxl d0,d0 addl IMM (SF_INFINITY),d0 bset IMM (22),d0 | Signaling NaN? jne Ldfsf$sign movel IMM (INVALID_OPERATION),d7 movel IMM (TRUNCDFSF+OP1_DOUBLE+DST_DOUBLE),sp@ jra $_exception_handler |============================================================================= | __fixdfsi |============================================================================= | int __fixdfsi (double) | Note that gcc wants round-to-zero mode here SYM (__fixdfsi): link a6,IMM (-24) moveml d2-d7,sp@ pea FIXDFSI+OP1_DOUBLE+DST_INT lea SYM (_fpCCR),a0 movel a6@(8),d0 | Get operand movel a6@(12),d1 movel d0,a1 | Isolate sign bclr IMM (31),d0 subl d0,a1 movel IMM (0xfff00000),d3 | Exponent mask movel d0,d4 andl d3,d4 addl d3,d4 | Adjust exponent subl d4,d0 | Fix up significand asrl IMM (4),d4 | Right justify the exponent swap d4 cmpl IMM (0x7fe),d4 | special case? jcc Ldfsi$spec | Yep, go handle subl IMM (DF_RBIAS-1),d4 | Less than 1/2? jlt Lgfsi$small | Yep, go show small cmpl IMM (32),d4 | >= 2**32? jgt Lgfsi$max | Yep, illegal subl IMM (21),d4 | Find how far to shift jeq Lgfsi$round | None at all? jlt Ldfsi$right | Go right-shift lsll d4,d0 | Align upper bits movel IMM (32),d5 subl d4,d5 movel d1,d3 lsrl d5,d3 orl d3,d0 | Carry over from lower bits lsll d4,d1 | Align lower bits jra Lgfsi$round Ldfsi$right: negl d4 | Get right-shift count movel IMM (32),d5 subl d4,d5 | and left-shift count movel d1,d3 lsrl d4,d1 | Align lower bits lsll d5,d3 | Anything down here? jeq 1f bset IMM (1),d1 | Set sticky bit 1: movel d0,d3 lsll d5,d3 | Carry over from upper bits orl d3,d1 lsrl d4,d0 | Align upper bits | Round & check for overflow Lgfsi$round: movel d1,d7 | Anything down here? jeq Ldfsi$to$zero movel IMM (INEXACT_RESULT),d7 | Set Inexact flag #ifdef FIX_ROUND_MODE movel IMM (1),d2 | Increment constant movel a0@(ROUND),d6 | Get rounding mode in d6 jeq Ldfsi$to$nearest subl IMM (ROUND_TO_PLUS),d6 jhi Ldfsi$to$minus jlt Ldfsi$to$zero | jra Ldfsi$to$plus Ldfsi$to$plus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jmi Ldfsi$to$zero | Yep, truncate jra Ldfsi$inc | Nope, increment Ldfsi$to$minus: #ifndef __mcf5200__ cmpa #0,a1 | Minus result? #else tstl a1 | Minus result? #endif jpl Ldfsi$to$zero | Nope, truncate jra Ldfsi$inc | Yep, increment Ldfsi$to$nearest: addl d1,d1 | Which way is nearest? jcc Ldfsi$to$zero jne Ldfsi$inc andl d0,d2 Ldfsi$inc: addl d2,d0 jcs Lgfsi$max | Overflow! #endif /* FIX_ROUND_MODE */ Ldfsi$to$zero: tstl d0 | Safe? jmi Ldfsi$mimnax | Maybe minimum? #ifndef __mcf5200__ cmpa #0,a1 | Minus? #else tstl a1 | Minus? #endif jpl Ldfsi$eck negl d0 | Return minus Ldfsi$eck: tstl d7 | Exact? jne $_exception_handler | Nope, check for traps Lgfsi$ret: moveml sp@(4),d2-d7 | Restore regs unlk a6 | and return rts Ldfsi$mimnax: #ifndef __mcf5200__ cmpa #0,a1 | Minus? #else tstl a1 | Minus? #endif jpl Lgfsi$max | Nope, Overflow addl d0,d0 | Min value? jne Lgfsi$max | Nope, Overflow movel a1,d0 | Yep, get it back jra Ldfsi$eck | and check for exceptions Ldfsi$spec: jeq Lgfsi$inop | +/-INFINITY or NaN? bclr IMM (20),d0 | Remove false hidden bit orl d0,d1 | Zero? jeq Lgfsi$ret | Yep, return zero | Summarize and round for non-zero |values| < 1/2 Lgfsi$small: movel IMM (0),d0 | Show very small movel IMM (1),d1 jra Lgfsi$round | Number won't fit in int, Invalid Operation Exception ' Lgfsi$max: movel a1,d0 | Get sign jmi Lgfsi$inv | If minus, it's min ' movel IMM (0x7fffffff),d0 | Get max for plus jra Lgfsi$inv Lgfsi$inop: movel IMM (0),d0 | Nothing good to return Lgfsi$inv: movel IMM (INVALID_OPERATION),d7 jra $_exception_handler |============================================================================= | __fixsfsi |============================================================================= | int __fixsfsi (float) | Note that gcc wants round-to-zero mode here SYM (__fixsfsi): link a6,IMM (-24) moveml d2-d7,sp@ pea FIXSFSI+OP1_FLOAT+DST_INT lea SYM (_fpCCR),a0 movel a6@(8),d0 | Get operand movel IMM (0),d1 movel d0,a1 | Isolate sign bclr IMM (31),d0 subl d0,a1 movel IMM (0xff800000),d3 | Exponent mask movel d0,d4 andl d3,d4 addl d3,d4 | Adjust exponent subl d4,d0 | Fix up significand asrl IMM (7),d4 | Right justify the exponent swap d4 cmpl IMM (0xfe),d4 | special case? jcc Lsfsi$spec | Yep, go handle subl IMM (SF_RBIAS-1),d4 | Less than 1/2? jlt Lgfsi$small | Yep, go show small cmpl IMM (32),d4 | >= 2**32? jgt Lgfsi$max | Yep, illegal subl IMM (24),d4 | Find how far to shift jeq Lgfsi$round | None at all? jlt Lsfsi$right | Go right-shift lsll d4,d0 | Align upper bits jra Lgfsi$round Lsfsi$right: negl d4 | Get right-shift count movel IMM (32),d5 subl d4,d5 | and left-shift count movel d0,d1 lsll d5,d1 | Carry over from upper bits lsrl d4,d0 | Align upper bits jra Lgfsi$round Lsfsi$ret: moveml sp@(4),d2-d7 | Restore regs unlk a6 | and return rts Lsfsi$spec: jeq Lgfsi$inop | +/-INFINITY or NaN? andl IMM (0x7fffff),d0 | Remove false hidden bit jeq Lgfsi$ret | Zero? jra Lgfsi$small | Go handle denormals #endif /* L_cvt */