lb1sf68.asm

GUN开源阻止下的编译器GCC

	.text
	.even

| These are common routines to return and signal exceptions.	

Ld$den:
| Return and signal a denormalized number
	orl	d7,d0
	movew	IMM (UNDERFLOW),d7
	orw	IMM (INEXACT_RESULT),d7
	movew	IMM (DOUBLE_FLOAT),d6
	jmp	$_exception_handler

Ld$infty:
Ld$overflow:
| Return a properly signed INFINITY and set the exception flags 
	movel	IMM (0x7ff00000),d0
	movel	IMM (0),d1
	orl	d7,d0
	movew	IMM (OVERFLOW),d7
	orw	IMM (INEXACT_RESULT),d7
	movew	IMM (DOUBLE_FLOAT),d6
	jmp	$_exception_handler

Ld$underflow:
| Return 0 and set the exception flags 
	movel	IMM (0),d0
	movel	d0,d1
	movew	IMM (UNDERFLOW),d7
	orw	IMM (INEXACT_RESULT),d7
	movew	IMM (DOUBLE_FLOAT),d6
	jmp	$_exception_handler

Ld$inop:
| Return a quiet NaN and set the exception flags
	movel	IMM (QUIET_NaN),d0
	movel	d0,d1
	movew	IMM (INVALID_OPERATION),d7
	orw	IMM (INEXACT_RESULT),d7
	movew	IMM (DOUBLE_FLOAT),d6
	jmp	$_exception_handler

Ld$div$0:
| Return a properly signed INFINITY and set the exception flags
	movel	IMM (0x7ff00000),d0
	movel	IMM (0),d1
	orl	d7,d0
	movew	IMM (DIVIDE_BY_ZERO),d7
	orw	IMM (INEXACT_RESULT),d7
	movew	IMM (DOUBLE_FLOAT),d6
	jmp	$_exception_handler

|=============================================================================
|=============================================================================
|                         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, shifted by 126 */
|  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):
	bchg	IMM (31),sp@(12) | change sign of second operand
				| and fall through, so we always add
|=============================================================================
|                              __adddf3
|=============================================================================

| double __adddf3(double, double);
SYM (__adddf3):
	link	a6,IMM (0)	| everything will be done in registers
	moveml	d2-d7,sp@-	| save all data registers and a2 (but d0-d1)
	movel	a6@(8),d0	| get first operand
	movel	a6@(12),d1	| 
	movel	a6@(16),d2	| get second operand
	movel	a6@(20),d3	| 

	movel	d0,d7		| get d0's sign bit in d7 '
	addl	d1,d1		| check and clear sign bit of a, and gain one
	addxl	d0,d0		| bit of extra precision
	beq	Ladddf$b	| if zero return second operand

	movel	d2,d6		| save sign in d6 
	addl	d3,d3		| get rid of sign bit and gain one bit of
	addxl	d2,d2		| extra precision
	beq	Ladddf$a	| if zero return first operand

	andl	IMM (0x80000000),d7 | isolate a's sign bit '
        swap	d6		| and also b's sign bit '
	andw	IMM (0x8000),d6	|
	orw	d6,d7		| and combine them into d7, so that a's sign '
				| bit is in the high word and b's is in the '
				| low word, so d6 is free to be used
	movel	d7,a0		| now save d7 into a0, so d7 is free to
                		| be used also

| Get the exponents and check for denormalized and/or infinity.

	movel	IMM (0x001fffff),d6 | mask for the fraction
	movel	IMM (0x00200000),d7 | mask to put hidden bit back

	movel	d0,d4		| 
	andl	d6,d0		| get fraction in d0
	notl	d6		| make d6 into mask for the exponent
	andl	d6,d4		| get exponent in d4
	beq	Ladddf$a$den	| branch if a is denormalized
	cmpl	d6,d4		| check for INFINITY or NaN
	beq	Ladddf$nf       | 
	orl	d7,d0		| and put hidden bit back
Ladddf$1:
	swap	d4		| shift right exponent so that it starts
	lsrw	IMM (5),d4	| in bit 0 and not bit 20
| Now we have a's exponent in d4 and fraction in d0-d1 '
	movel	d2,d5		| save b to get exponent
	andl	d6,d5		| get exponent in d5
	beq	Ladddf$b$den	| branch if b is denormalized
	cmpl	d6,d5		| check for INFINITY or NaN
	beq	Ladddf$nf
	notl	d6		| make d6 into mask for the fraction again
	andl	d6,d2		| and get fraction in d2
	orl	d7,d2		| and put hidden bit back
Ladddf$2:
	swap	d5		| shift right exponent so that it starts
	lsrw	IMM (5),d5	| in bit 0 and not bit 20

| Now we have b's exponent in d5 and fraction in d2-d3. '

| The situation now is as follows: the signs are combined in a0, the 
| numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
| and d5 (b). To do the rounding correctly we need to keep all the
| bits until the end, so we need to use d0-d1-d2-d3 for the first number
| and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
| exponents in a2-a3.

	moveml	a2-a3,sp@-	| save the address registers

	movel	d4,a2		| save the exponents
	movel	d5,a3		| 

	movel	IMM (0),d7	| and move the numbers around
	movel	d7,d6		|
	movel	d3,d5		|
	movel	d2,d4		|
	movel	d7,d3		|
	movel	d7,d2		|

| Here we shift the numbers until the exponents are the same, and put 
| the largest exponent in a2.
	exg	d4,a2		| get exponents back
	exg	d5,a3		|
	cmpw	d4,d5		| compare the exponents
	beq	Ladddf$3	| if equal don't shift '
	bhi	9f		| branch if second exponent is higher

| Here we have a's exponent larger than b's, so we have to shift b. We do 
| this by using as counter d2:
1:	movew	d4,d2		| move largest exponent to d2
	subw	d5,d2		| and subtract second exponent
	exg	d4,a2		| get back the longs we saved
	exg	d5,a3		|
| if difference is too large we don't shift (actually, we can just exit) '
	cmpw	IMM (DBL_MANT_DIG+2),d2
	bge	Ladddf$b$small
	cmpw	IMM (32),d2	| if difference >= 32, shift by longs
	bge	5f
2:	cmpw	IMM (16),d2	| if difference >= 16, shift by words	
	bge	6f
	bra	3f		| enter dbra loop

4:	lsrl	IMM (1),d4
	roxrl	IMM (1),d5
	roxrl	IMM (1),d6
	roxrl	IMM (1),d7
3:	dbra	d2,4b
	movel	IMM (0),d2
	movel	d2,d3	
	bra	Ladddf$4
5:
	movel	d6,d7
	movel	d5,d6
	movel	d4,d5
	movel	IMM (0),d4
	subw	IMM (32),d2
	bra	2b
6:
	movew	d6,d7
	swap	d7
	movew	d5,d6
	swap	d6
	movew	d4,d5
	swap	d5
	movew	IMM (0),d4
	swap	d4
	subw	IMM (16),d2
	bra	3b
	
9:	exg	d4,d5
	movew	d4,d6
	subw	d5,d6		| keep d5 (largest exponent) in d4
	exg	d4,a2
	exg	d5,a3
| if difference is too large we don't shift (actually, we can just exit) '
	cmpw	IMM (DBL_MANT_DIG+2),d6
	bge	Ladddf$a$small
	cmpw	IMM (32),d6	| if difference >= 32, shift by longs
	bge	5f
2:	cmpw	IMM (16),d6	| if difference >= 16, shift by words	
	bge	6f
	bra	3f		| enter dbra loop

4:	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
3:	dbra	d6,4b
	movel	IMM (0),d7
	movel	d7,d6
	bra	Ladddf$4
5:
	movel	d2,d3
	movel	d1,d2
	movel	d0,d1
	movel	IMM (0),d0
	subw	IMM (32),d6
	bra	2b
6:
	movew	d2,d3
	swap	d3
	movew	d1,d2
	swap	d2
	movew	d0,d1
	swap	d1
	movew	IMM (0),d0
	swap	d0
	subw	IMM (16),d6
	bra	3b
Ladddf$3:
	exg	d4,a2	
	exg	d5,a3
Ladddf$4:	
| Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
| the signs in a4.

| Here we have to decide whether to add or subtract the numbers:
	exg	d7,a0		| get the signs 
	exg	d6,a3		| a3 is free to be used
	movel	d7,d6		|
	movew	IMM (0),d7	| get a's sign in d7 '
	swap	d6              |
	movew	IMM (0),d6	| and b's sign in d6 '
	eorl	d7,d6		| compare the signs
	bmi	Lsubdf$0	| if the signs are different we have 
				| to subtract
	exg	d7,a0		| else we add the numbers
	exg	d6,a3		|
	addl	d7,d3		|
	addxl	d6,d2		|
	addxl	d5,d1		| 
	addxl	d4,d0           |

	movel	a2,d4		| return exponent to d4
	movel	a0,d7		| 
	andl	IMM (0x80000000),d7 | d7 now has the sign

	moveml	sp@+,a2-a3	

| Before rounding normalize so bit #DBL_MANT_DIG 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:
	btst	IMM (DBL_MANT_DIG+1),d0	
	beq	1f
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
	addw	IMM (1),d4
1:
	lea	Ladddf$5,a0	| to return from rounding routine
	lea	SYM (_fpCCR),a1	| check the rounding mode
	movew	a1@(6),d6	| rounding mode in d6
	beq	Lround$to$nearest
	cmpw	IMM (ROUND_TO_PLUS),d6
	bhi	Lround$to$minus
	blt	Lround$to$zero
	bra	Lround$to$plus
Ladddf$5:
| Put back the exponent and check for overflow
	cmpw	IMM (0x7ff),d4	| is the exponent big?
	bge	1f
	bclr	IMM (DBL_MANT_DIG-1),d0
	lslw	IMM (4),d4	| put exponent back into position
	swap	d0		| 
	orw	d4,d0		|
	swap	d0		|
	bra	Ladddf$ret
1:
	movew	IMM (ADD),d5
	bra	Ld$overflow

Lsubdf$0:
| Here we do the subtraction.
	exg	d7,a0		| put sign back in a0
	exg	d6,a3		|
	subl	d7,d3		|
	subxl	d6,d2		|
	subxl	d5,d1		|
	subxl	d4,d0		|
	beq	Ladddf$ret$1	| if zero just exit
	bpl	1f		| if positive skip the following
	exg	d7,a0		|
	bchg	IMM (31),d7	| change sign bit in d7
	exg	d7,a0		|
	negl	d3		|
	negxl	d2		|
	negxl	d1              | and negate result
	negxl	d0              |
1:	
	movel	a2,d4		| return exponent to d4
	movel	a0,d7
	andl	IMM (0x80000000),d7 | isolate sign bit
	moveml	sp@+,a2-a3	|

| Before rounding normalize so bit #DBL_MANT_DIG 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:
	btst	IMM (DBL_MANT_DIG+1),d0	
	beq	1f
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
	addw	IMM (1),d4
1:
	lea	Lsubdf$1,a0	| to return from rounding routine
	lea	SYM (_fpCCR),a1	| check the rounding mode
	movew	a1@(6),d6	| rounding mode in d6
	beq	Lround$to$nearest
	cmpw	IMM (ROUND_TO_PLUS),d6
	bhi	Lround$to$minus
	blt	Lround$to$zero
	bra	Lround$to$plus
Lsubdf$1:
| Put back the exponent and sign (we don't have overflow). '
	bclr	IMM (DBL_MANT_DIG-1),d0	
	lslw	IMM (4),d4	| put exponent back into position
	swap	d0		| 
	orw	d4,d0		|
	swap	d0		|
	bra	Ladddf$ret

| If one of the numbers was too small (difference of exponents >= 
| DBL_MANT_DIG+1) we return the other (and now we don't have to '
| check for finiteness or zero).
Ladddf$a$small:
	moveml	sp@+,a2-a3	
	movel	a6@(16),d0
	movel	a6@(20),d1
	lea	SYM (_fpCCR),a0
	movew	IMM (0),a0@
	moveml	sp@+,d2-d7	| restore data registers
	unlk	a6		| and return
	rts

Ladddf$b$small:
	moveml	sp@+,a2-a3	
	movel	a6@(8),d0
	movel	a6@(12),d1
	lea	SYM (_fpCCR),a0
	movew	IMM (0),a0@
	moveml	sp@+,d2-d7	| restore data registers
	unlk	a6		| and return
	rts

Ladddf$a$den:
	movel	d7,d4		| d7 contains 0x00200000
	bra	Ladddf$1

Ladddf$b$den:
	movel	d7,d5           | d7 contains 0x00200000
	notl	d6
	bra	Ladddf$2

Ladddf$b:
| Return b (if a is zero)
	movel	d2,d0
	movel	d3,d1
	bra	1f
Ladddf$a:
	movel	a6@(8),d0
	movel	a6@(12),d1
1:
	movew	IMM (ADD),d5
| Check for NaN and +/-INFINITY.
	movel	d0,d7         		|
	andl	IMM (0x80000000),d7	|
	bclr	IMM (31),d0		|
	cmpl	IMM (0x7ff00000),d0	|
	bge	2f			|
	movel	d0,d0           	| check for zero, since we don't  '
	bne	Ladddf$ret		| want to return -0 by mistake
	bclr	IMM (31),d7		|
	bra	Ladddf$ret		|
2:
	andl	IMM (0x000fffff),d0	| check for NaN (nonzero fraction)
	orl	d1,d0			|
	bne	Ld$inop         	|
	bra	Ld$infty		|
	
Ladddf$ret$1:
	moveml	sp@+,a2-a3	| restore regs and exit

Ladddf$ret:
| Normal exit.
	lea	SYM (_fpCCR),a0
	movew	IMM (0),a0@
	orl	d7,d0		| put sign bit back
	moveml	sp@+,d2-d7
	unlk	a6
	rts

Ladddf$ret$den:
| Return a denormalized number.
	lsrl	IMM (1),d0	| shift right once more
	roxrl	IMM (1),d1	|
	bra	Ladddf$ret

Ladddf$nf:
	movew	IMM (ADD),d5
| This could be faster but it is not worth the effort, since it is not
| executed very often. We sacrifice speed for clarity here.
	movel	a6@(8),d0	| get the numbers back (remember that we
	movel	a6@(12),d1	| did some processing already)
	movel	a6@(16),d2	| 
	movel	a6@(20),d3	| 
	movel	IMM (0x7ff00000),d4 | useful constant (INFINITY)
	movel	d0,d7		| save sign bits
	movel	d2,d6		| 
	bclr	IMM (31),d0	| clear sign bits
	bclr	IMM (31),d2	| 
| We know that one of them is either NaN of +/-INFINITY