lb1sf68.asm

Mac OS X 10.4.9 for x86 Source Code gcc 实现源代码

#else
	movel	sp@+,a4
	movel	sp@+,a3
	movel	sp@+,a2
#endif

| 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
#ifndef __mcoldfire__
	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
#else
	moveq	IMM (29),d6
	lsrl	IMM (3),d3
	movel	d2,d7
	lsll	d6,d7
	orl	d7,d3
	lsrl	IMM (3),d2
	movel	d1,d7
	lsll	d6,d7
	orl	d7,d2
	lsrl	IMM (3),d1
	movel	d0,d7
	lsll	d6,d7
	orl	d7,d1
	lsrl	IMM (3),d0
#endif
	
| 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
#ifndef __mcoldfire__
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	addw	IMM (1),d4
#else
	lsrl	IMM (1),d1
	btst	IMM (0),d0
	beq	10f
	bset	IMM (31),d1
10:	lsrl	IMM (1),d0
	addl	IMM (1),d4
#endif
	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
#ifndef __mcoldfire__
	exg	d2,d0		| put b (==0) into d0-d1
	exg	d3,d1		| and a (with sign bit cleared) into d2-d3
#else
	movel	d2,d7
	movel	d0,d2
	movel	d7,d0
	movel	d3,d7
	movel	d1,d3
	movel	d7,d1
#endif
	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
	PICLEA	SYM (_fpCCR),a0
	movew	IMM (0),a0@
#ifndef __mcoldfire__
	moveml	sp@+,d2-d7
#else
	moveml	sp@,d2-d7
	| XXX if frame pointer is ever removed, stack pointer must
	| be adjusted here.
#endif
	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		|
#ifndef __mcoldfire__
	subw	IMM (1),d4	| and adjust exponent
#else
	subl	IMM (1),d4	| and adjust exponent
#endif
	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		|
#ifndef __mcoldfire__
	subw	IMM (1),d5	| and adjust exponent
#else
	subql	IMM (1),d5	| and adjust exponent
#endif
	btst	IMM (20),d2	|
	bne	Lmuldf$2	|
	bra	1b


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

| double __divdf3(double, double);
SYM (__divdf3):
#ifndef __mcoldfire__
	link	a6,IMM (0)
	moveml	d2-d7,sp@-
#else
	link	a6,IMM (-24)
	moveml	d2-d7,sp@
#endif
	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
#ifndef __mcoldfire__
	lsrw	IMM (4),d4	| 
#else
	lsrl	IMM (4),d4	| 
#endif
Ldivdf$1:			| 
	andl	d7,d5		|
	beq	Ldivdf$b$den	|
	andl	d6,d2		|
	orl	IMM (0x00100000),d2
	swap	d5		|
#ifndef __mcoldfire__
	lsrw	IMM (4),d5	|
#else
	lsrl	IMM (4),d5	|
#endif
Ldivdf$2:			|
#ifndef __mcoldfire__
	subw	d5,d4		| subtract exponents
	addw	IMM (D_BIAS),d4	| and add bias
#else
	subl	d5,d4		| subtract exponents
	addl	IMM (D_BIAS),d4	| and add bias
#endif

| 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		|
#ifndef __mcoldfire__
	dbra	d5,1b		| and branch back
#else
	subql	IMM (1), d5
	bpl	1b
#endif
	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		|
#ifndef __mcoldfire__
	dbra	d5,1b		| and branch back
#else
	subql	IMM (1), d5
	bpl	1b
#endif
	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		|
#ifndef __mcoldfire__
	dbra	d5,1b		| and branch back
#else
	subql	IMM (1), d5
	bpl	1b
#endif
	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
#ifndef __mcoldfire__
	subw	IMM (DBL_MANT_DIG),d5
	addw	IMM (63),d5
	cmpw	IMM (31),d5
#else
	subl	IMM (DBL_MANT_DIG),d5
	addl	IMM (63),d5
	cmpl	IMM (31),d5
#endif
	bhi	2f
1:	bset	d5,d3
	bra	5f
#ifndef __mcoldfire__
	subw	IMM (32),d5
#else
	subl	IMM (32),d5
#endif
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
#ifndef __mcoldfire__
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
	addw	IMM (1),d4
#else
	lsrl	IMM (1),d3
	btst	IMM (0),d2
	beq	10f
	bset	IMM (31),d3
10:	lsrl	IMM (1),d2
	btst	IMM (0),d1
	beq	11f
	bset	IMM (31),d2
11:	lsrl	IMM (1),d1
	btst	IMM (0),d0
	beq	12f
	bset	IMM (31),d1
12:	lsrl	IMM (1),d0
	addl	IMM (1),d4
#endif
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		| 
	PICLEA	SYM (_fpCCR),a0	| clear exception flags
	movew	IMM (0),a0@	|
#ifndef __mcoldfire__
	moveml	sp@+,d2-d7	| 
#else
	moveml	sp@,d2-d7	| 
	| XXX if frame pointer is ever removed, stack pointer must
	| be adjusted here.
#endif
	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

| If a number is denormalized we put an exponent of 1 but do not put the 
| bit back into the fraction.
Ldivdf$a$den:
	movel	IMM (1),d4
	andl	d6,d0
1:	addl	d1,d1		| shift a left until bit 20 is set
	addxl	d0,d0
#ifndef __mcoldfire__
	subw	IMM (1),d4	| and adjust exponent
#else
	subl	IMM (1),d4	| and adjust exponent
#endif
	btst	IMM (DBL_MANT_DIG-32-1),d0
	bne	Ldivdf$1
	bra	1b

Ldivdf$b$den:
	movel	IMM (1),d5
	andl	d6,d2
1:	addl	d3,d3		| shift b left until bit 20 is set
	addxl	d2,d2
#ifndef __mcoldfire__
	subw	IMM (1),d5	| and adjust exponent
#else
	subql	IMM (1),d5	| and adjust exponent
#endif
	btst	IMM (DBL_MANT_DIG-32-1),d2
	bne	Ldivdf$2
	bra	1b

Lround$exit:
| This is a common exit point for __muldf3 and __divdf3. When they enter
| this point the sign of the result is in d7, the result in d0-d1, normalized
| so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.

| First check for underlow in the exponent:
#ifndef __mcoldfire__
	cmpw	IMM (-DBL_MANT_DIG-1),d4		
#else
	cmpl	IMM (-DBL_MANT_DIG-1),d4		
#endif
	blt	Ld$underflow	
| It could happen that the exponent is less than 1, in which case the 
| number is denormalized. In this case we shift right and adjust the 
| exponent until it becomes 1 or the fraction is zero (in the latter case 
| we signal underflow and return zero).
	movel	d7,a0		|
	movel	IMM (0),d6	| use d6-d7 to collect bits flushed right
	movel	d6,d7		| use d6-d7 to collect bits flushed right
#ifndef __mcoldfire__
	cmpw	IMM (1),d4	| if the exponent is less than 1 we 
#else
	cmpl	IMM (1),d4	| if the exponent is less than 1 we 
#endif
	bge	2f		| have to shift right (denormalize)
1:
#ifndef __mcoldfire__
	addw	IMM (1),d4	| adjust the exponent
	lsrl	IMM (1),d0	| shift right once 
	roxrl	IMM (1),d1	|
	roxrl	IMM (1),d2	|
	roxrl	IMM (1),d3	|
	roxrl	IMM (1),d6	| 
	roxrl	IMM (1),d7	|
	cmpw	IMM (1),d4	| is the exponent 1 already?
#else
	addl	IMM (1),d4	| adjust the exponent
	lsrl	IMM (1),d7
	btst	IMM (0),d6
	beq	13f
	bset	IMM (31),d7
13:	lsrl	IMM (1),d6
	btst	IMM (0),d3
	beq	14f
	bset	IMM (31),d6
14:	lsrl	IMM (1),d3
	btst	IMM (0),d2
	beq	10f
	bset	IMM (31),d3
10:	lsrl	IMM (1),d2
	btst	IMM (0),d1
	beq	11f
	bset	IMM (31),d2
11:	lsrl	IMM (1),d1
	btst	IMM (0),d0
	beq	12f
	bset	IMM (31),d1
12:	lsrl	IMM (1),d0
	cmpl	IMM (1),d4	| is the exponent 1 already?
#endif