#APP

 | take floating point to integer and fractional pieces
 | 
 | C interface
 |  double modf( double value, double *iptr)
 |  returns fractional part of value
 |	   in *iptr returns the integral part
 |	   such that (*iptr + fractional) == value
 | 
 |-----------------------------------------------------------------------------
 | ported to 68000 by Kai-Uwe Bloem, 12/89
 |  #1  original author: Peter S. Housel 9/21/88,01/17/89,03/19/89,5/24/89
 |  #2  replaced shifts by swap if possible for speed increase	-kub-, 01/90
 |  #3  ported to gcc ++jrb 03/90
 |-----------------------------------------------------------------------------

BIAS8	=	0x3ff - 1

	.text; .even
	.globl _modf
_modf:
	lea	sp@(4),a0	| a0 -> float argument
	movel	sp@(12),a1	| a1 -> ipart result
	moveml	d2-d7,sp@-	| save d2-d7

	movew	a0@,d0		| extract value.exp
	movew	d0,d2		| extract value.sign
	lsrw	#4,d0
	andw	#0x7ff,d0	| kill sign bit

	cmpw	#BIAS8,d0
	bge	1f		| fabs(value) >= 1.0

	clrl	a1@		| store zero as the integer part
	clrl	a1@(4)
	moveml	a0@,d0-d1	| return entire value as fractional part
	jra	L0

1:
	cmpw	#BIAS8+53,d0	| all integer, with no fractional part ?
	blt	2f		| no, mixed

	movel	a0@,a1@		| store entire value as the integer part
	movel	a0@(4),a1@(4)
	moveq	#0,d0		| return zero as fractional part
	moveq	#0,d1
L0:
	moveml	sp@+,d2-d7	| restore saved d2-d7
	rts
2:
|	moveml	d3-d7,sp@-	| save some registers (d2 is pushed already)
	moveml	a0@,d4-d5	| get value

	andl	#0x0fffff,d4	| remove exponent from value.mantissa
	orl	#0x100000,d4	| restore implied leading "1"

	moveq	#0,d6		| zero fractional part
	moveq	#0,d7
3:
	cmpw	#BIAS8+37,d0	| fast shift, 16 bits ?
	bgt	4f
	movew	d6,d7		| shift down 16 bits
	movew	d5,d6
	movew	d4,d5
	clrw	d4
	swap	d7
	swap	d6
	swap	d5
	swap	d4
	addw	#16,d0
	bra	3b
4:
	cmpw	#BIAS8+53,d0	| done ?
	bge	5f
	lsrl	#1,d4		| shift integer part
	roxrl	#1,d5

	roxrl	#1,d6		| shift high bit into fractional part
	roxrl	#1,d7

	addw	#1,d0		| increment ipart exponent
	bra	4b		| keep shifting
5:
	| normalize ipart (all values are in correct reggies)
 	| save  a1, save d2-d7 that norm_df will pop
	movel	a1,sp@-
	pea	L1		| set up return address
	moveml	d2-d7,sp@-	| norm_df will pop this
	clrw	d1
	jmp	norm_df		| go do it
L1:				| norm_df will rts to here
	movel	sp@+,a1		| pop saved a1
	moveml	d0-d1,a1@	| store result into ipart
	
	| norm fractional part
	movel	d6,d4		| get frac into d4/d5
	movel	d7,d5
	clrw	d1		| rounding = 0
	movew	#BIAS8-11,d0	| set frac part exponent, sign already in d2
	jmp	norm_df		| norm_df will pop d2/d7 we save before
				| it will return to our caller via rts
				| with result in d0-d1
