lb1sf68.asm

GUN开源阻止下的编译器GCC

| Check for NaN (if either one is NaN return NaN)
	cmpl	d4,d0		| check first a (d0)
	bhi	Ld$inop		| if d0 > 0x7ff00000 or equal and
	bne	2f
	tstl	d1		| d1 > 0, a is NaN
	bne	Ld$inop		| 
2:	cmpl	d4,d2		| check now b (d1)
	bhi	Ld$inop		| 
	bne	3f
	tstl	d3		| 
	bne	Ld$inop		| 
3:
| Now comes the check for +/-INFINITY. We know that both are (maybe not
| finite) numbers, but we have to check if both are infinite whether we
| are adding or subtracting them.
	eorl	d7,d6		| to check sign bits
	bmi	1f
	andl	IMM (0x80000000),d7 | get (common) sign bit
	bra	Ld$infty
1:
| We know one (or both) are infinite, so we test for equality between the
| two numbers (if they are equal they have to be infinite both, so we
| return NaN).
	cmpl	d2,d0		| are both infinite?
	bne	1f		| if d0 <> d2 they are not equal
	cmpl	d3,d1		| if d0 == d2 test d3 and d1
	beq	Ld$inop		| if equal return NaN
1:	
	andl	IMM (0x80000000),d7 | get a's sign bit '
	cmpl	d4,d0		| test now for infinity
	beq	Ld$infty	| if a is INFINITY return with this sign
	bchg	IMM (31),d7	| else we know b is INFINITY and has
	bra	Ld$infty	| the opposite sign

|=============================================================================
|                              __muldf3
|=============================================================================

| double __muldf3(double, double);
SYM (__muldf3):
	link	a6,IMM (0)
	moveml	d2-d7,sp@-
	movel	a6@(8),d0		| get a into d0-d1
	movel	a6@(12),d1		| 
	movel	a6@(16),d2		| and b into d2-d3
	movel	a6@(20),d3		|
	movel	d0,d7			| d7 will hold the sign of the product
	eorl	d2,d7			|
	andl	IMM (0x80000000),d7	|
	movel	d7,a0			| save sign bit into a0 
	movel	IMM (0x7ff00000),d7	| useful constant (+INFINITY)
	movel	d7,d6			| another (mask for fraction)
	notl	d6			|
	bclr	IMM (31),d0		| get rid of a's sign bit '
	movel	d0,d4			| 
	orl	d1,d4			| 
	beq	Lmuldf$a$0		| branch if a is zero
	movel	d0,d4			|
	bclr	IMM (31),d2		| get rid of b's sign bit '
	movel	d2,d5			|
	orl	d3,d5			| 
	beq	Lmuldf$b$0		| branch if b is zero
	movel	d2,d5			| 
	cmpl	d7,d0			| is a big?
	bhi	Lmuldf$inop		| if a is NaN return NaN
	beq	Lmuldf$a$nf		| we still have to check d1 and b ...
	cmpl	d7,d2			| now compare b with INFINITY
	bhi	Lmuldf$inop		| is b NaN?
	beq	Lmuldf$b$nf 		| we still have to check d3 ...
| Here we have both numbers finite and nonzero (and with no sign bit).
| Now we get the exponents into d4 and d5.
	andl	d7,d4			| isolate exponent in d4
	beq	Lmuldf$a$den		| if exponent zero, have denormalized
	andl	d6,d0			| isolate fraction
	orl	IMM (0x00100000),d0	| and put hidden bit back
	swap	d4			| I like exponents in the first byte
	lsrw	IMM (4),d4		| 
Lmuldf$1:			
	andl	d7,d5			|
	beq	Lmuldf$b$den		|
	andl	d6,d2			|
	orl	IMM (0x00100000),d2	| and put hidden bit back
	swap	d5			|
	lsrw	IMM (4),d5		|
Lmuldf$2:				|
	addw	d5,d4			| add exponents
	subw	IMM (D_BIAS+1),d4	| and subtract bias (plus one)

| 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!), which means that in the product bit 104 
| (which will correspond to bit 8 of the fourth long) is set.

| Here we have to do the product.
| To do it we have to juggle the registers back and forth, as there are not
| enough to keep everything in them. So we use the address registers to keep
| some intermediate data.

	moveml	a2-a3,sp@-	| save a2 and a3 for temporary use
	movel	IMM (0),a2	| a2 is a null register
	movel	d4,a3		| and a3 will preserve the exponent

| First, shift d2-d3 so bit 20 becomes bit 31:
	rorl	IMM (5),d2	| rotate d2 5 places right
	swap	d2		| and swap it
	rorl	IMM (5),d3	| do the same thing with d3
	swap	d3		|
	movew	d3,d6		| get the rightmost 11 bits of d3
	andw	IMM (0x07ff),d6	|
	orw	d6,d2		| and put them into d2
	andw	IMM (0xf800),d3	| clear those bits in d3

	movel	d2,d6		| move b into d6-d7
	movel	d3,d7           | move a into d4-d5
	movel	d0,d4           | and clear d0-d1-d2-d3 (to put result)
	movel	d1,d5           |
	movel	IMM (0),d3	|
	movel	d3,d2           |
	movel	d3,d1           |
	movel	d3,d0	        |

| We use a1 as counter:	
	movel	IMM (DBL_MANT_DIG-1),a1		
	exg	d7,a1

1:	exg	d7,a1		| put counter back in a1
	addl	d3,d3		| shift sum once left
	addxl	d2,d2           |
	addxl	d1,d1           |
	addxl	d0,d0           |
	addl	d7,d7		|
	addxl	d6,d6		|
	bcc	2f		| if bit clear skip the following
	exg	d7,a2		|
	addl	d5,d3		| else add a to the sum
	addxl	d4,d2		|
	addxl	d7,d1		|
	addxl	d7,d0		|
	exg	d7,a2		| 
2:	exg	d7,a1		| put counter in d7
	dbf	d7,1b		| decrement and branch

	movel	a3,d4		| restore exponent
	moveml	sp@+,a2-a3

| Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The 
| first thing to do now is to normalize it so bit 8 becomes bit 
| DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
	swap	d0
	swap	d1
	movew	d1,d0
	swap	d2
	movew	d2,d1
	swap	d3
	movew	d3,d2
	movew	IMM (0),d3
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
	
| Now round, check for over- and underflow, and exit.
	movel	a0,d7		| get sign bit back into d7
	movew	IMM (MULTIPLY),d5

	btst	IMM (DBL_MANT_DIG+1-32),d0
	beq	Lround$exit
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	addw	IMM (1),d4
	bra	Lround$exit

Lmuldf$inop:
	movew	IMM (MULTIPLY),d5
	bra	Ld$inop

Lmuldf$b$nf:
	movew	IMM (MULTIPLY),d5
	movel	a0,d7		| get sign bit back into d7
	tstl	d3		| we know d2 == 0x7ff00000, so check d3
	bne	Ld$inop		| if d3 <> 0 b is NaN
	bra	Ld$overflow	| else we have overflow (since a is finite)

Lmuldf$a$nf:
	movew	IMM (MULTIPLY),d5
	movel	a0,d7		| get sign bit back into d7
	tstl	d1		| we know d0 == 0x7ff00000, so check d1
	bne	Ld$inop		| if d1 <> 0 a is NaN
	bra	Ld$overflow	| else signal overflow

| If either number is zero return zero, unless the other is +/-INFINITY or
| NaN, in which case we return NaN.
Lmuldf$b$0:
	movew	IMM (MULTIPLY),d5
	exg	d2,d0		| put b (==0) into d0-d1
	exg	d3,d1		| and a (with sign bit cleared) into d2-d3
	bra	1f
Lmuldf$a$0:
	movel	a6@(16),d2	| put b into d2-d3 again
	movel	a6@(20),d3	|
	bclr	IMM (31),d2	| clear sign bit
1:	cmpl	IMM (0x7ff00000),d2 | check for non-finiteness
	bge	Ld$inop		| in case NaN or +/-INFINITY return NaN
	lea	SYM (_fpCCR),a0
	movew	IMM (0),a0@
	moveml	sp@+,d2-d7
	unlk	a6
	rts

| 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.
Lmuldf$a$den:
	movel	IMM (1),d4
	andl	d6,d0
1:	addl	d1,d1           | shift a left until bit 20 is set
	addxl	d0,d0		|
	subw	IMM (1),d4	| and adjust exponent
	btst	IMM (20),d0	|
	bne	Lmuldf$1        |
	bra	1b

Lmuldf$b$den:
	movel	IMM (1),d5
	andl	d6,d2
1:	addl	d3,d3		| shift b left until bit 20 is set
	addxl	d2,d2		|
	subw	IMM (1),d5	| and adjust exponent
	btst	IMM (20),d2	|
	bne	Lmuldf$2	|
	bra	1b


|=============================================================================
|                              __divdf3
|=============================================================================

| double __divdf3(double, double);
SYM (__divdf3):
	link	a6,IMM (0)
	moveml	d2-d7,sp@-
	movel	a6@(8),d0	| get a into d0-d1
	movel	a6@(12),d1	| 
	movel	a6@(16),d2	| and b into d2-d3
	movel	a6@(20),d3	|
	movel	d0,d7		| d7 will hold the sign of the result
	eorl	d2,d7		|
	andl	IMM (0x80000000),d7
	movel	d7,a0		| save sign into a0
	movel	IMM (0x7ff00000),d7 | useful constant (+INFINITY)
	movel	d7,d6		| another (mask for fraction)
	notl	d6		|
	bclr	IMM (31),d0	| get rid of a's sign bit '
	movel	d0,d4		|
	orl	d1,d4		|
	beq	Ldivdf$a$0	| branch if a is zero
	movel	d0,d4		|
	bclr	IMM (31),d2	| get rid of b's sign bit '
	movel	d2,d5		|
	orl	d3,d5		|
	beq	Ldivdf$b$0	| branch if b is zero
	movel	d2,d5
	cmpl	d7,d0		| is a big?
	bhi	Ldivdf$inop	| if a is NaN return NaN
	beq	Ldivdf$a$nf	| if d0 == 0x7ff00000 we check d1
	cmpl	d7,d2		| now compare b with INFINITY 
	bhi	Ldivdf$inop	| if b is NaN return NaN
	beq	Ldivdf$b$nf	| if d2 == 0x7ff00000 we check d3
| Here we have both numbers finite and nonzero (and with no sign bit).
| Now we get the exponents into d4 and d5 and normalize the numbers to
| ensure that the ratio of the fractions is around 1. We do this by
| making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
| set, even if they were denormalized to start with.
| Thus, the result will satisfy: 2 > result > 1/2.
	andl	d7,d4		| and isolate exponent in d4
	beq	Ldivdf$a$den	| if exponent is zero we have a denormalized
	andl	d6,d0		| and isolate fraction
	orl	IMM (0x00100000),d0 | and put hidden bit back
	swap	d4		| I like exponents in the first byte
	lsrw	IMM (4),d4	| 
Ldivdf$1:			| 
	andl	d7,d5		|
	beq	Ldivdf$b$den	|
	andl	d6,d2		|
	orl	IMM (0x00100000),d2
	swap	d5		|
	lsrw	IMM (4),d5	|
Ldivdf$2:			|
	subw	d5,d4		| subtract exponents
	addw	IMM (D_BIAS),d4	| and add bias

| We are now ready to do the division. We have prepared things in such a way
| that the ratio of the fractions will be less than 2 but greater than 1/2.
| At this point the registers in use are:
| d0-d1	hold a (first operand, bit DBL_MANT_DIG-32=0, bit 
| DBL_MANT_DIG-1-32=1)
| d2-d3	hold b (second operand, bit DBL_MANT_DIG-32=1)
| d4	holds the difference of the exponents, corrected by the bias
| a0	holds the sign of the ratio

| To do the rounding correctly we need to keep information about the
| nonsignificant bits. One way to do this would be to do the division
| using four registers; another is to use two registers (as originally
| I did), but use a sticky bit to preserve information about the 
| fractional part. Note that we can keep that info in a1, which is not
| used.
	movel	IMM (0),d6	| d6-d7 will hold the result
	movel	d6,d7		| 
	movel	IMM (0),a1	| and a1 will hold the sticky bit

	movel	IMM (DBL_MANT_DIG-32+1),d5	
	
1:	cmpl	d0,d2		| is a < b?
	bhi	3f		| if b > a skip the following
	beq	4f		| if d0==d2 check d1 and d3
2:	subl	d3,d1		| 
	subxl	d2,d0		| a <-- a - b
	bset	d5,d6		| set the corresponding bit in d6
3:	addl	d1,d1		| shift a by 1
	addxl	d0,d0		|
	dbra	d5,1b		| and branch back
	bra	5f			
4:	cmpl	d1,d3		| here d0==d2, so check d1 and d3
	bhi	3b		| if d1 > d2 skip the subtraction
	bra	2b		| else go do it
5:
| Here we have to start setting the bits in the second long.
	movel	IMM (31),d5	| again d5 is counter

1:	cmpl	d0,d2		| is a < b?
	bhi	3f		| if b > a skip the following
	beq	4f		| if d0==d2 check d1 and d3
2:	subl	d3,d1		| 
	subxl	d2,d0		| a <-- a - b
	bset	d5,d7		| set the corresponding bit in d7
3:	addl	d1,d1		| shift a by 1
	addxl	d0,d0		|
	dbra	d5,1b		| and branch back
	bra	5f			
4:	cmpl	d1,d3		| here d0==d2, so check d1 and d3
	bhi	3b		| if d1 > d2 skip the subtraction
	bra	2b		| else go do it
5:
| Now go ahead checking until we hit a one, which we store in d2.
	movel	IMM (DBL_MANT_DIG),d5
1:	cmpl	d2,d0		| is a < b?
	bhi	4f		| if b < a, exit
	beq	3f		| if d0==d2 check d1 and d3
2:	addl	d1,d1		| shift a by 1
	addxl	d0,d0		|
	dbra	d5,1b		| and branch back
	movel	IMM (0),d2	| here no sticky bit was found
	movel	d2,d3
	bra	5f			
3:	cmpl	d1,d3		| here d0==d2, so check d1 and d3
	bhi	2b		| if d1 > d2 go back
4:
| Here put the sticky bit in d2-d3 (in the position which actually corresponds
| to it; if you don't do this the algorithm loses in some cases). '
	movel	IMM (0),d2
	movel	d2,d3
	subw	IMM (DBL_MANT_DIG),d5
	addw	IMM (63),d5
	cmpw	IMM (31),d5
	bhi	2f
1:	bset	d5,d3
	bra	5f
	subw	IMM (32),d5
2:	bset	d5,d2
5:
| Finally we are finished! Move the longs in the address registers to
| their final destination:
	movel	d6,d0
	movel	d7,d1
	movel	IMM (0),d3

| Here we have finished the division, with the result in d0-d1-d2-d3, with
| 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
| If it is not, then definitely bit 21 is set. Normalize so bit 22 is
| not set:
	btst	IMM (DBL_MANT_DIG-32+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:
| Now round, check for over- and underflow, and exit.
	movel	a0,d7		| restore sign bit to d7
	movew	IMM (DIVIDE),d5
	bra	Lround$exit

Ldivdf$inop:
	movew	IMM (DIVIDE),d5
	bra	Ld$inop

Ldivdf$a$0:
| If a is zero check to see whether b is zero also. In that case return
| NaN; then check if b is NaN, and return NaN also in that case. Else
| return zero.
	movew	IMM (DIVIDE),d5
	bclr	IMM (31),d2	|
	movel	d2,d4		| 
	orl	d3,d4		| 
	beq	Ld$inop		| if b is also zero return NaN
	cmpl	IMM (0x7ff00000),d2 | check for NaN
	bhi	Ld$inop		| 
	blt	1f		|
	tstl	d3		|
	bne	Ld$inop		|
1:	movel	IMM (0),d0	| else return zero
	movel	d0,d1		| 
	lea	SYM (_fpCCR),a0	| clear exception flags
	movew	IMM (0),a0@	|
	moveml	sp@+,d2-d7	| 
	unlk	a6		| 
	rts			| 	

Ldivdf$b$0:
	movew	IMM (DIVIDE),d5
| If we got here a is not zero. Check if a is NaN; in that case return NaN,
| else return +/-INFINITY. Remember that a is in d0 with the sign bit 
| cleared already.
	movel	a0,d7		| put a's sign bit back in d7 '
	cmpl	IMM (0x7ff00000),d0 | compare d0 with INFINITY
	bhi	Ld$inop		| if larger it is NaN
	tstl	d1		| 
	bne	Ld$inop		| 
	bra	Ld$div$0	| else signal DIVIDE_BY_ZERO

Ldivdf$b$nf:
	movew	IMM (DIVIDE),d5
| If d2 == 0x7ff00000 we have to check d3.
	tstl	d3		|
	bne	Ld$inop		| if d3 <> 0, b is NaN
	bra	Ld$underflow	| else b is +/-INFINITY, so signal underflow

Ldivdf$a$nf:
	movew	IMM (DIVIDE),d5
| If d0 == 0x7ff00000 we have to check d1.
	tstl	d1		|
	bne	Ld$inop		| if d1 <> 0, a is NaN
| If a is INFINITY we have to check b
	cmpl	d7,d2		| compare b with INFINITY 
	bge	Ld$inop		| if b is NaN or INFINITY return NaN
	tstl	d3		|
	bne	Ld$inop		| 
	bra	Ld$overflow	| else return overflow