lb1sf68.asm

gcc-you can use this code to learn something about gcc, and inquire further into linux,

	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 __mcf5200__
	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
	lea	SYM (_fpCCR),a0
	movew	IMM (0),a0@
#ifndef __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	lsrw	IMM (4),d5	|
#else
	lsrl	IMM (4),d5	|
#endif
Ldivdf$2:			|
#ifndef __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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		| 
	lea	SYM (_fpCCR),a0	| clear exception flags
	movew	IMM (0),a0@	|
#ifndef __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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 __mcf5200__
	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
	beq	2f		| if not loop back
	bra	1b              |
	bra	Ld$underflow	| safety check, shouldn't execute '
2:	orl	d6,d2		| this is a trick so we don't lose  '
	orl	d7,d3		| the bits which were flushed right
	movel	a0,d7		| get back sign bit into d7
| Now call the rounding routine (which takes care of denormalized numbers):
	lea	Lround$0,a0	| to return from rounding routine
	lea	SYM (_fpCCR),a1	| check the rounding mode
#ifdef __mcf5200__
	clrl	d6
#endif
	movew	a1@(6),d6	| rounding mode in d6
	beq	Lround$to$nearest
#ifndef __mcf5200__
	cmpw	IMM (ROUND_TO_PLUS),d6
#else
	cmpl	IMM (ROUND_TO_PLUS),d6
#endif
	bhi	Lround$to$minus
	blt	Lround$to$zero
	bra	Lround$to$plus
Lround$0:
| Here we have a correctly rounded result (either normalized or denormalized).

| Here we should have either a normalized number or a denormalized one, and
| the exponent is necessarily larger or equal to 1 (so we don't have to  '
| check again for underflow!). We have to check for overflow or for a 
| denormalized number (which also signals underflow).
| Check for overflow (i.e., exponent >= 0x7ff).
#ifndef __mcf5200__
	cmpw	IMM (0x07ff),d4
#else
	cmpl	IMM (0x07ff),d4
#endif
	bge	Ld$overflow
| Now check for a denormalized number (exponent==0):
	movew	d4,d4
	beq	Ld$den
1:
| Put back the exponents and sign and return.
#ifndef __mcf5200__
	lslw	IMM (4),d4	| exponent back to fourth byte
#else
	lsll	IMM (4),d4	| exponent back to fourth byte
#endif
	bclr	IMM (DBL_MANT_DIG-32-1),d0
	swap	d0		| and put back exponent
#ifndef __mcf5200__
	orw	d4,d0		| 
#else
	orl	d4,d0		| 
#endif
	swap	d0		|
	orl	d7,d0		| and sign also

	lea	SYM (_fpCCR),a0
	movew	IMM (0),a0@
#ifndef __mcf5200__
	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

|=============================================================================
|                              __negdf2
|=============================================================================

| double __negdf2(double, double);
SYM (__negdf2):
#ifndef __mcf5200__
	link	a6,IMM (0)
	moveml	d2-d7,sp@-
#else
	link	a6,IMM (-24)
	moveml	d2-d7,sp@
#endif
	movew	IMM (NEGATE),d5
	movel	a6@(8),d0	| get number to negate in d0-d1
	movel	a6@(12),d1	|
	bchg	IMM (31),d0	| negate
	movel	d0,d2		| make a positive copy (for the tests)