lb1sf68.asm

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

	bclr	IMM (31),d2	|
	movel	d2,d4		| check for zero
	orl	d1,d4		|
	beq	2f		| if zero (either sign) return +zero
	cmpl	IMM (0x7ff00000),d2 | compare to +INFINITY
	blt	1f		| if finite, return
	bhi	Ld$inop		| if larger (fraction not zero) is NaN
	tstl	d1		| if d2 == 0x7ff00000 check d1
	bne	Ld$inop		|
	movel	d0,d7		| else get sign and return INFINITY
	andl	IMM (0x80000000),d7
	bra	Ld$infty		
1:	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
2:	bclr	IMM (31),d0
	bra	1b

|=============================================================================
|                              __cmpdf2
|=============================================================================

GREATER =  1
LESS    = -1
EQUAL   =  0

| int __cmpdf2(double, double);
SYM (__cmpdf2):
#ifndef __mcf5200__
	link	a6,IMM (0)
	moveml	d2-d7,sp@- 	| save registers
#else
	link	a6,IMM (-24)
	moveml	d2-d7,sp@
#endif
	movew	IMM (COMPARE),d5
	movel	a6@(8),d0	| get first operand
	movel	a6@(12),d1	|
	movel	a6@(16),d2	| get second operand
	movel	a6@(20),d3	|
| First check if a and/or b are (+/-) zero and in that case clear
| the sign bit.
	movel	d0,d6		| copy signs into d6 (a) and d7(b)
	bclr	IMM (31),d0	| and clear signs in d0 and d2
	movel	d2,d7		|
	bclr	IMM (31),d2	|
	cmpl	IMM (0x7fff0000),d0 | check for a == NaN
	bhi	Ld$inop		| if d0 > 0x7ff00000, a is NaN
	beq	Lcmpdf$a$nf	| if equal can be INFINITY, so check d1
	movel	d0,d4		| copy into d4 to test for zero
	orl	d1,d4		|
	beq	Lcmpdf$a$0	|
Lcmpdf$0:
	cmpl	IMM (0x7fff0000),d2 | check for b == NaN
	bhi	Ld$inop		| if d2 > 0x7ff00000, b is NaN
	beq	Lcmpdf$b$nf	| if equal can be INFINITY, so check d3
	movel	d2,d4		|
	orl	d3,d4		|
	beq	Lcmpdf$b$0	|
Lcmpdf$1:
| Check the signs
	eorl	d6,d7
	bpl	1f
| If the signs are not equal check if a >= 0
	tstl	d6
	bpl	Lcmpdf$a$gt$b	| if (a >= 0 && b < 0) => a > b
	bmi	Lcmpdf$b$gt$a	| if (a < 0 && b >= 0) => a < b
1:
| If the signs are equal check for < 0
	tstl	d6
	bpl	1f
| If both are negative exchange them
#ifndef __mcf5200__
	exg	d0,d2
	exg	d1,d3
#else
	movel	d0,d7
	movel	d2,d0
	movel	d7,d2
	movel	d1,d7
	movel	d3,d1
	movel	d7,d3
#endif
1:
| Now that they are positive we just compare them as longs (does this also
| work for denormalized numbers?).
	cmpl	d0,d2
	bhi	Lcmpdf$b$gt$a	| |b| > |a|
	bne	Lcmpdf$a$gt$b	| |b| < |a|
| If we got here d0 == d2, so we compare d1 and d3.
	cmpl	d1,d3
	bhi	Lcmpdf$b$gt$a	| |b| > |a|
	bne	Lcmpdf$a$gt$b	| |b| < |a|
| If we got here a == b.
	movel	IMM (EQUAL),d0
#ifndef __mcf5200__
	moveml	sp@+,d2-d7 	| put back the registers
#else
	moveml	sp@,d2-d7
	| XXX if frame pointer is ever removed, stack pointer must
	| be adjusted here.
#endif
	unlk	a6
	rts
Lcmpdf$a$gt$b:
	movel	IMM (GREATER),d0
#ifndef __mcf5200__
	moveml	sp@+,d2-d7 	| put back the registers
#else
	moveml	sp@,d2-d7
	| XXX if frame pointer is ever removed, stack pointer must
	| be adjusted here.
#endif
	unlk	a6
	rts
Lcmpdf$b$gt$a:
	movel	IMM (LESS),d0
#ifndef __mcf5200__
	moveml	sp@+,d2-d7 	| put back the registers
#else
	moveml	sp@,d2-d7
	| XXX if frame pointer is ever removed, stack pointer must
	| be adjusted here.
#endif
	unlk	a6
	rts

Lcmpdf$a$0:	
	bclr	IMM (31),d6
	bra	Lcmpdf$0
Lcmpdf$b$0:
	bclr	IMM (31),d7
	bra	Lcmpdf$1

Lcmpdf$a$nf:
	tstl	d1
	bne	Ld$inop
	bra	Lcmpdf$0

Lcmpdf$b$nf:
	tstl	d3
	bne	Ld$inop
	bra	Lcmpdf$1

|=============================================================================
|                           rounding routines
|=============================================================================

| The rounding routines expect the number to be normalized in registers
| d0-d1-d2-d3, with the exponent in register d4. They assume that the 
| exponent is larger or equal to 1. They return a properly normalized number
| if possible, and a denormalized number otherwise. The exponent is returned
| in d4.

Lround$to$nearest:
| We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
| Here we assume that the exponent is not too small (this should be checked
| before entering the rounding routine), but the number could be denormalized.

| Check for denormalized numbers:
1:	btst	IMM (DBL_MANT_DIG-32),d0
	bne	2f		| if set the number is normalized
| Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent 
| is one (remember that a denormalized number corresponds to an 
| exponent of -D_BIAS+1).
#ifndef __mcf5200__
	cmpw	IMM (1),d4	| remember that the exponent is at least one
#else
	cmpl	IMM (1),d4	| remember that the exponent is at least one
#endif
 	beq	2f		| an exponent of one means denormalized
	addl	d3,d3		| else shift and adjust the exponent
	addxl	d2,d2		|
	addxl	d1,d1		|
	addxl	d0,d0		|
#ifndef __mcf5200__
	dbra	d4,1b		|
#else
	subql	IMM (1), d4
	bpl	1b
#endif
2:
| Now round: we do it as follows: after the shifting we can write the
| fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
| If delta < 1, do nothing. If delta > 1, add 1 to f. 
| If delta == 1, we make sure the rounded number will be even (odd?) 
| (after shifting).
	btst	IMM (0),d1	| is delta < 1?
	beq	2f		| if so, do not do anything
	orl	d2,d3		| is delta == 1?
	bne	1f		| if so round to even
	movel	d1,d3		| 
	andl	IMM (2),d3	| bit 1 is the last significant bit
	movel	IMM (0),d2	|
	addl	d3,d1		|
	addxl	d2,d0		|
	bra	2f		| 
1:	movel	IMM (1),d3	| else add 1 
	movel	IMM (0),d2	|
	addl	d3,d1		|
	addxl	d2,d0
| Shift right once (because we used bit #DBL_MANT_DIG-32!).
2:
#ifndef __mcf5200__
	lsrl	IMM (1),d0
	roxrl	IMM (1),d1		
#else
	lsrl	IMM (1),d1
	btst	IMM (0),d0
	beq	10f
	bset	IMM (31),d1
10:	lsrl	IMM (1),d0
#endif

| Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
| 'fraction overflow' ...).
	btst	IMM (DBL_MANT_DIG-32),d0	
	beq	1f
#ifndef __mcf5200__
	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
1:
| If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we 
| have to put the exponent to zero and return a denormalized number.
	btst	IMM (DBL_MANT_DIG-32-1),d0
	beq	1f
	jmp	a0@
1:	movel	IMM (0),d4
	jmp	a0@

Lround$to$zero:
Lround$to$plus:
Lround$to$minus:
	jmp	a0@
#endif /* L_double */

#ifdef  L_float

	.globl	SYM (_fpCCR)
	.globl  $_exception_handler

QUIET_NaN    = 0xffffffff
SIGNL_NaN    = 0x7f800001
INFINITY     = 0x7f800000

F_MAX_EXP      = 0xff
F_BIAS         = 126
FLT_MAX_EXP    = F_MAX_EXP - F_BIAS
FLT_MIN_EXP    = 1 - F_BIAS
FLT_MANT_DIG   = 24

INEXACT_RESULT 		= 0x0001
UNDERFLOW 		= 0x0002
OVERFLOW 		= 0x0004
DIVIDE_BY_ZERO 		= 0x0008
INVALID_OPERATION 	= 0x0010

SINGLE_FLOAT = 1

NOOP         = 0
ADD          = 1
MULTIPLY     = 2
DIVIDE       = 3
NEGATE       = 4
COMPARE      = 5
EXTENDSFDF   = 6
TRUNCDFSF    = 7

UNKNOWN           = -1
ROUND_TO_NEAREST  = 0 | round result to nearest representable value
ROUND_TO_ZERO     = 1 | round result towards zero
ROUND_TO_PLUS     = 2 | round result towards plus infinity
ROUND_TO_MINUS    = 3 | round result towards minus infinity

| Entry points:

	.globl SYM (__addsf3)
	.globl SYM (__subsf3)
	.globl SYM (__mulsf3)
	.globl SYM (__divsf3)
	.globl SYM (__negsf2)
	.globl SYM (__cmpsf2)

| These are common routines to return and signal exceptions.	

	.text
	.even

Lf$den:
| Return and signal a denormalized number
	orl	d7,d0
	movew	IMM (INEXACT_RESULT+UNDERFLOW),d7
	moveq	IMM (SINGLE_FLOAT),d6
	jmp	$_exception_handler

Lf$infty:
Lf$overflow:
| Return a properly signed INFINITY and set the exception flags 
	movel	IMM (INFINITY),d0
	orl	d7,d0
	movew	IMM (INEXACT_RESULT+OVERFLOW),d7
	moveq	IMM (SINGLE_FLOAT),d6
	jmp	$_exception_handler

Lf$underflow:
| Return 0 and set the exception flags 
	movel	IMM (0),d0
	movew	IMM (INEXACT_RESULT+UNDERFLOW),d7
	moveq	IMM (SINGLE_FLOAT),d6
	jmp	$_exception_handler

Lf$inop:
| Return a quiet NaN and set the exception flags
	movel	IMM (QUIET_NaN),d0
	movew	IMM (INEXACT_RESULT+INVALID_OPERATION),d7
	moveq	IMM (SINGLE_FLOAT),d6
	jmp	$_exception_handler

Lf$div$0:
| Return a properly signed INFINITY and set the exception flags
	movel	IMM (INFINITY),d0
	orl	d7,d0
	movew	IMM (INEXACT_RESULT+DIVIDE_BY_ZERO),d7
	moveq	IMM (SINGLE_FLOAT),d6
	jmp	$_exception_handler

|=============================================================================
|=============================================================================
|                         single precision routines
|=============================================================================
|=============================================================================

| A single precision floating point number (float) has the format:
|
| struct _float {
|  unsigned int sign      : 1;  /* sign bit */ 
|  unsigned int exponent  : 8;  /* exponent, shifted by 126 */
|  unsigned int fraction  : 23; /* fraction */
| } float;
| 
| Thus sizeof(float) = 4 (32 bits). 
|
| All the routines are callable from C programs, and return the result 
| in the single register d0. They also preserve all registers except 
| d0-d1 and a0-a1.

|=============================================================================
|                              __subsf3
|=============================================================================

| float __subsf3(float, float);
SYM (__subsf3):
	bchg	IMM (31),sp@(8)	| change sign of second operand
				| and fall through
|=============================================================================
|                              __addsf3
|=============================================================================

| float __addsf3(float, float);
SYM (__addsf3):
#ifndef __mcf5200__
	link	a6,IMM (0)	| everything will be done in registers
	moveml	d2-d7,sp@-	| save all data registers but d0-d1
#else
	link	a6,IMM (-24)
	moveml	d2-d7,sp@
#endif
	movel	a6@(8),d0	| get first operand
	movel	a6@(12),d1	| get second operand
	movel	d0,d6		| get d0's sign bit '
	addl	d0,d0		| check and clear sign bit of a
	beq	Laddsf$b	| if zero return second operand
	movel	d1,d7		| save b's sign bit '
	addl	d1,d1		| get rid of sign bit
	beq	Laddsf$a	| if zero return first operand

	movel	d6,a0		| save signs in address registers
	movel	d7,a1		| so we can use d6 and d7

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

	movel	IMM (0x00ffffff),d4	| mask to get fraction
	movel	IMM (0x01000000),d5	| mask to put hidden bit back

	movel	d0,d6		| save a to get exponent
	andl	d4,d0		| get fraction in d0
	notl 	d4		| make d4 into a mask for the exponent
	andl	d4,d6		| get exponent in d6
	beq	Laddsf$a$den	| branch if a is denormalized
	cmpl	d4,d6		| check for INFINITY or NaN
	beq	Laddsf$nf
	swap	d6		| put exponent into first word
	orl	d5,d0		| and put hidden bit back
Laddsf$1:
| Now we have a's exponent in d6 (second byte) and the mantissa in d0. '
	movel	d1,d7		| get exponent in d7
	andl	d4,d7		| 
	beq	Laddsf$b$den	| branch if b is denormalized
	cmpl	d4,d7		| check for INFINITY or NaN
	beq	Laddsf$nf
	swap	d7		| put exponent into first word
	notl 	d4		| make d4 into a mask for the fraction
	andl	d4,d1		| get fraction in d1
	orl	d5,d1		| and put hidden bit back
Laddsf$2:
| Now we have b's exponent in d7 (second byte) and the mantissa in d1. '

| Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we 
| shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
| bit).

	movel	d1,d2		| move b to d2, since we want to use
				| two registers to do the sum
	movel	IMM (0),d1	| and clear the new ones
	movel	d1,d3		|

| Here we shift the numbers in registers d0 and d1 so the exponents are the
| same, and put the largest exponent in d6. Note that we are using two
| registers for each number (see the discussion by D. Knuth in "Seminumerical 
| Algorithms").
#ifndef __mcf5200__
	cmpw	d6,d7		| compare exponents
#else
	cmpl	d6,d7		| compare exponents
#endif
	beq	Laddsf$3	| if equal don't shift '
	bhi	5f		| branch if second exponent largest
1:
	subl	d6,d7		| keep the largest exponent
	negl	d7
#ifndef __mcf5200__
	lsrw	IMM (8),d7	| put difference in lower byte
#else
	lsrl	IMM (8),d7	| put difference in lower byte
#endif
| if difference is too large we don't shift (actually, we can just exit) '
#ifndef __mcf5200__
	cmpw	IMM (FLT_MANT_DIG+2),d7		
#else
	cmpl	IMM (FLT_MANT_DIG+2),d7		
#endif
	bge	Laddsf$b$small
#ifndef __mcf5200__
	cmpw	IMM (16),d7	| if difference >= 16 swap
#else
	cmpl	IMM (16),d7	| if difference >= 16 swap
#endif
	bge	4f
2:
#ifndef __mcf5200__
	subw	IMM (1),d7
#else
	subql	IMM (1), d7
#endif
3:
#ifndef __mcf5200__
	lsrl	IMM (1),d2	| shift right second operand
	roxrl	IMM (1),d3
	dbra	d7,3b
#else
	lsrl	IMM (1),d3
	btst	IMM (0),d2
	beq	10f
	bset	IMM (31),d3
10:	lsrl	IMM (1),d2
	subql	IMM (1), d7
	bpl	3b
#endif
	bra	Laddsf$3
4:
	movew	d2,d3
	swap	d3
	movew	d3,d2
	swap	d2
#ifndef __mcf5200__
	subw	IMM (16),d7
#else
	subl	IMM (16),d7
#endif
	bne	2b		| if still more bits, go back to normal case
	bra	Laddsf$3
5:
#ifndef __mcf5200__
	exg	d6,d7		| exchange the exponents
#else
	eorl	d6,d7
	eorl	d7,d6
	eorl	d6,d7
#endif
	subl	d6,d7		| keep the largest exponent
	negl	d7		|
#ifndef __mcf5200__
	lsrw	IMM (8),d7	| put difference in lower byte
#else
	lsrl	IMM (8),d7	| put difference in lower byte
#endif
| if difference is too large we don't shift (and exit!) '
#ifndef __mcf5200__
	cmpw	IMM (FLT_MANT_DIG+2),d7		
#else
	cmpl	IMM (FLT_MANT_DIG+2),d7		
#endif
	bge	Laddsf$a$small
#ifndef __mcf5200__
	cmpw	IMM (16),d7	| if difference >= 16 swap
#else
	cmpl	IMM (16),d7	| if difference >= 16 swap
#endif
	bge	8f
6:
#ifndef