lb1sf68.asm

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

	.globl SYM (__negdf2)
	.globl SYM (__cmpdf2)

	.text
	.even

| These are common routines to return and signal exceptions.	

Ld$den:
| Return and signal a denormalized number
	orl	d7,d0
	movew	IMM (INEXACT_RESULT+UNDERFLOW),d7
	moveq	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 (INEXACT_RESULT+OVERFLOW),d7
	moveq	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 (INEXACT_RESULT+UNDERFLOW),d7
	moveq	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 (INEXACT_RESULT+INVALID_OPERATION),d7
	moveq	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 (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
	moveq	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):
#ifndef __mcf5200__
	link	a6,IMM (0)	| everything will be done in registers
	moveml	d2-d7,sp@-	| save all data registers and a2 (but d0-d1)
#else
	link	a6,IMM (-24)
	moveml	d2-d7,sp@
#endif
	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 '
#ifndef __mcf5200__
	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
#else
	andl	IMM (0x8000),d6
	orl	d6,d7
#endif
	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
#ifndef __mcf5200__
	lsrw	IMM (5),d4	| in bit 0 and not bit 20
#else
	lsrl	IMM (5),d4	| in bit 0 and not bit 20
#endif
| 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
#ifndef __mcf5200__
	lsrw	IMM (5),d5	| in bit 0 and not bit 20
#else
	lsrl	IMM (5),d5	| in bit 0 and not bit 20
#endif

| 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.

#ifndef __mcf5200__
	moveml	a2-a3,sp@-	| save the address registers
#else
	movel	a2,sp@-	
	movel	a3,sp@-	
	movel	a4,sp@-	
#endif

	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.
#ifndef __mcf5200__
	exg	d4,a2		| get exponents back
	exg	d5,a3		|
	cmpw	d4,d5		| compare the exponents
#else
	movel	d4,a4		| get exponents back
	movel	a2,d4
	movel	a4,a2
	movel	d5,a4
	movel	a3,d5
	movel	a4,a3
	cmpl	d4,d5		| compare the exponents
#endif
	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
#ifndef __mcf5200__
	subw	d5,d2		| and subtract second exponent
	exg	d4,a2		| get back the longs we saved
	exg	d5,a3		|
#else
	subl	d5,d2		| and subtract second exponent
	movel	d4,a4		| get back the longs we saved
	movel	a2,d4
	movel	a4,a2
	movel	d5,a4
	movel	a3,d5
	movel	a4,a3
#endif
| if difference is too large we don't shift (actually, we can just exit) '
#ifndef __mcf5200__
	cmpw	IMM (DBL_MANT_DIG+2),d2
#else
	cmpl	IMM (DBL_MANT_DIG+2),d2
#endif
	bge	Ladddf$b$small
#ifndef __mcf5200__
	cmpw	IMM (32),d2	| if difference >= 32, shift by longs
#else
	cmpl	IMM (32),d2	| if difference >= 32, shift by longs
#endif
	bge	5f
2:
#ifndef __mcf5200__
	cmpw	IMM (16),d2	| if difference >= 16, shift by words	
#else
	cmpl	IMM (16),d2	| if difference >= 16, shift by words	
#endif
	bge	6f
	bra	3f		| enter dbra loop

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

4:
#ifndef __mcf5200__
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1
	roxrl	IMM (1),d2
	roxrl	IMM (1),d3
#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
#endif
3:
#ifndef __mcf5200__
	dbra	d6,4b
#else
	subql	IMM (1),d6
	bpl	4b
#endif
	movel	IMM (0),d7
	movel	d7,d6
	bra	Ladddf$4
5:
	movel	d2,d3
	movel	d1,d2
	movel	d0,d1
	movel	IMM (0),d0
#ifndef __mcf5200__
	subw	IMM (32),d6
#else
	subl	IMM (32),d6
#endif
	bra	2b
6:
	movew	d2,d3
	swap	d3
	movew	d1,d2
	swap	d2
	movew	d0,d1
	swap	d1
	movew	IMM (0),d0
	swap	d0
#ifndef __mcf5200__
	subw	IMM (16),d6
#else
	subl	IMM (16),d6
#endif
	bra	3b
Ladddf$3:
#ifndef __mcf5200__
	exg	d4,a2	
	exg	d5,a3
#else
	movel	d4,a4
	movel	a2,d4
	movel	a4,a2
	movel	d5,a4
	movel	a3,d5
	movel	a4,a3
#endif
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:
#ifndef __mcf5200__
	exg	d7,a0		| get the signs 
	exg	d6,a3		| a3 is free to be used
#else
	movel	d7,a4
	movel	a0,d7
	movel	a4,a0
	movel	d6,a4
	movel	a3,d6
	movel	a4,a3
#endif
	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
#ifndef __mcf5200__
	exg	d7,a0		| else we add the numbers
	exg	d6,a3		|
#else
	movel	d7,a4
	movel	a0,d7
	movel	a4,a0
	movel	d6,a4
	movel	a3,d6
	movel	a4,a3
#endif
	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

#ifndef __mcf5200__
	moveml	sp@+,a2-a3	
#else
	movel	sp@+,a4	
	movel	sp@+,a3	
	movel	sp@+,a2	
#endif

| 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
#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:
	lea	Ladddf$5,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
Ladddf$5:
| Put back the exponent and check for overflow
#ifndef __mcf5200__
	cmpw	IMM (0x7ff),d4	| is the exponent big?
#else
	cmpl	IMM (0x7ff),d4	| is the exponent big?
#endif
	bge	1f
	bclr	IMM (DBL_MANT_DIG-1),d0
#ifndef __mcf5200__
	lslw	IMM (4),d4	| put exponent back into position
#else
	lsll	IMM (4),d4	| put exponent back into position
#endif
	swap	d0		| 
#ifndef __mcf5200__
	orw	d4,d0		|
#else
	orl	d4,d0		|
#endif
	swap	d0		|
	bra	Ladddf$ret
1:
	movew	IMM (ADD),d5
	bra	Ld$overflow