| @(#)fract-mc68020.s 9.2 88/01/19	Copyright (c) 1985 by Sun Microsystems, Inc.

| fract.s -- a set of routines for dealing with fixed point numbers.
|	A "fract" is representes as a 32 bit integer.
|	The bottom 16 bits are the fraction, the top 16 are the integer

| Originally written sometime in 1984 by Vaughan Pratt.
| Hacked over, clean up, commented, and made to deal with overflow
| by James Gosling, February, 1985

	.text
	.globl	_ceilingfr, _roundfr, _floorfr, _floatfr, _fractf
	.globl	_frmul, _frdiv, _fridivu, _frsdivu, _frsqrt
|	.globl	_getdouble, _getfract
	.globl	_fract_overflows	| Count of the number of overflows
					| that have occurred.
	.globl	_vfrmul, _vfrdiv, _vfradd, _vfrsub
|	.globl	ieeeused


_ceilingfr:
	addl	#0x7fff,sp@(4)
_roundfr:
	addl	#0x8000,sp@(4)
_floorfr:
	movl	sp@(4),d0
|	clrw	d0
	swap	d0
	extl	d0
	rts

| Compute A*B, A,B fracts
|	k = 2**16
|	A = Ah+Al/k		(Ah and Al are the two 16 bit halves)
|	B = Bh+Bl/k
| A*B	= (Ah+Al/k)*(Bh+Bl/k)
|	= Al*Bl/(k**2) + Ah*Bh + Bh*Al + Ah*Bl
_frmul:				| A*B	A at sp@(4)m B at sp@(8)
|#ifdef MC68020
|	movl	sp@(4),d1
|	mulsl	sp@(8),d0:d1
|	swap	d0
|	swap	d1
|	movw	d1,d0
|#else
	subl	a0,a0
	tstl	sp@(4)		| A := abs(A)
	jge	L1
	negl	sp@(4)
	addql	#1,a0
L1:
	tstl	sp@(8)		| B := abs(B)
	jge	L2
	negl	sp@(8)
	subql	#1,a0		| a0 == 0 iff A&B have the same sign.
L2:
	movw	sp@(6),d0
	mulu	sp@(10),d0	| Al*Bl
	clrw	d0
	swap	d0		| top 16 of Al*Bl to bottom 16 of d0
	movw	sp@(4),d1
	mulu	sp@(8),d1	| Ah*Bh
	swap	d1
	clrw	d1		| bottom 16 of Ah*Bh to top 16 of d1
	addl	d1,d0
	movw	sp@(6),d1
	mulu	sp@(8),d1	| Al*Bh
	addl	d1,d0
	movw	sp@(4),d1
	mulu	sp@(10),d1	| Ah*Bl
	addl	d1,d0
	cmpl	#0,a0
	beq	RTS
	negl	d0
|#endif
	rts

_floatfr:
	subl	a0,a0
	clrl	d1
	movl	sp@(4),d0
	jeq	RTS
	jpl	L4
	addl	#0x800,a0
	negl	d0
	jvc	L4
	movl	#0xc0e00000,d0
	rts
L4:
	subql	#1,d1
	lsll	#1,d0
	jpl	L4
L5:
	lsll	#1,d0
	lsrl	#8,d0
	lsrl	#4,d0
	addl	#0x40e,d1
	addl	a0,d1
	rorl	#8,d1
	rorl	#4,d1
	addl	d1,d0
	clrl	d1
RTS:
	rts
_fractf:
	subl	a0,a0
	movl	sp@(4),d0
	jpl	L7
	addql	#1,a0
	andl	#0x7fffffff,d0
L7:
	jeq	L10
	andl	#0xfffff,d0
	orl	#0x100000,d0
	movl	sp@(8),d1
	andl	#0xffe00000,d1
	orl	d1,d0
	roll	#8,d0
	roll	#3,d0
	movw	sp@(4),d1
	andw	#0x7ff0,d1
	lsrw	#4,d1
	subw	#0x40e,d1
	jlt	L8
	movl	#0x7fffffff,d0
	movl	a0,d1
	jeq	RTS
	notl	d0
	rts
L8:
	negl	d1
	cmpw	#32,d1
	jlt	L9
	clrl	d0
L9:
	lsrl	d1,d0
L10:
	movl	a0,d1
	jeq	RTS1
	negl	d0
	rts

_vfrdiv:			| THIS IS A LIE

_frdiv:
	clrl	sp@-
	tstl	sp@(8)
	jge	L12
	negl	sp@(8)
	addql	#1,sp@
L12:
	tstl	sp@(12)
	jge	L13
	negl	sp@(12)
	subql	#1,sp@
L13:
	movl	sp@(12),sp@-
	movl	sp@(12),sp@-
	jsr	_fridivu	| floor(a/b)
	addql	#8,sp
	movl	d0,a1
	movl	sp@(12),sp@-
	movl	d0,sp@-
	jsr	_frmul		| floor(a/b)*b
	addql	#8,sp
	subl	sp@(8),d0	| floor(a/b)*b - a
	negl	d0		| a - floor(a/b)*b
	jge	L14		| floor(a/b) might be 1 too big
	addl	sp@(12),d0	| so patch by adding b
	subl	#0x10000,a1	| and decrement floor(a/b)
L14:
	movl	sp@(12),sp@-
	movl	d0,sp@-
	jsr	_frsdivu	| (a - floor(a/b)*b)/b
	addql	#8,sp
	addl	a1,d0
	tstl	sp@+
	beq	L14a
	negl	d0
L14a:
	rts
| fract small divide unsigned: a/b assuming a < b, both unsigned
_frsdivu:
	movl	sp@(4),d0	| a	unsigned
	movl	sp@(8),d1	| b	unsigned, > a
	jra	L16
L15:
	rorl	#1,d0		| this works because 0 <= a < b
	lsrl	#1,d1
L16:
	cmpl	#0x10000,d1	| get b down to 16 bits
	jge	L15
	cmpw	d1,d0
	jcs	L16a
	movl	#0x10000,d0
	rts
L16a:	swap	d0
	divu	d1,d0
	swap	d0
	clrw	d0		| clear remainder
	swap	d0
RTS1:
	rts
| fract integer divide unsigned:  floor(a/b), a,b unsigned fracts
_fridivu:
	movl	sp@(4),d0	| a	unsigned
	movl	sp@(8),d1	| b	unsigned
	jra	L18
L17:
	lsrl	#1,d0
	lsrl	#1,d1
L18:
	cmpl	#0x10000,d1
	jge	L17
	divu	d1,d0
	swap	d0
	clrw	d0
	rts
_frsqrt:
	movl	sp@(4),d0
	cmpl	#1,d0
	jhi	L19		| 0 <= d <= 1/65536 -> special case
	lsll	#8,d0
	rts
L19:
	exg	a0,d2
	exg	a1,d3
	movl	d4,sp@-
	cmpl	#0xFFFE0000,d0	| special-case high values
	jhi	ffff		| FFFE0001 <= d <= FFFFFFFF -> sqrt(d) = FFFF
	jeq	fffe		| d = FFFE0000  ->  sqrt(d) = FFFE
	movl	d0,d1		| d
	movl	d1,d2		| x
half:	lsrl	#1,d2		| Take initial approximation to
	lsrl	#2,d1		|   be input shifted
	jne	half		|   right halfway
	addqw	#1,d2		| fixes the 40000000 <= d < FFFE0000 problem
clr:	clrl	d3		| place to put quotient
	jra	entry		| enter loop
more:	addl	d3,d2		| x+d/x -> x
	lsrl	#1,d2		| x+d/x)/2
entry:	movl	d0,d1		| Make copy of d
	divu	d2,d1		| d/x
	movw	d1,d3		| extract quotient
	cmpl	d3,d2		| compare x with d/x
	jgt	more		| if x > d/x then reiterate
	movl	d2,d4		| move x to form
	addql	#2,d4		|   x+2
	cmpl	d4,d3		| compare d/x with x+2
	jlt	done		| if d/x < x+2 then done
	jgt	more		| if d/x > x+2 then reiterate
	eorl	d3,d1		| extract remainder
	jne	more		| if remainder ne 0, d/x > x+2 so reiterate
done:	movw	d2,d1		| Final stage: get low order byte of sqrt
	mulu	d2,d1
	subl	d1,d0		| d - x^2 (<= 2x)
	lsll	#7,d0		| <= 256x
	divu	d2,d0		| floor(128(d - x^2)/x) <= 256
	lsll	#8,d2
	extl	d0
	addl	d0,d2
	movl	d2,d0
sqexit:
	exg	a0,d2
	exg	a1,d3
	movl	sp@+,d4
	rts
fffe:	movl	#0xFFFE,d0
	jra	done
ffff:	movl	#0xFFFF,d0
	jra	done
|_getfract:
|	link	a6,#-4
|	pea	a6@(-4)
|	pea	L25
|	movl	a6@(8),sp@-
|	jbsr	_fscanf
|	addqw	#8,sp
|	movl	a6@(-4),d0
|	jsr	fvdoublei
|	movl	d1,sp@-
|	movl	d0,sp@-
|	jbsr	_fractf
|	addqw	#8,sp
|	unlk	a6
|	rts
|_getdouble:
|	pea	sp@(-8)
|	pea	L26
|	movl	sp@(12),sp@-
|	jbsr	_fscanf
|	lea	sp@(12),sp
|	movl	sp@(-8),d0
|	movl	sp@(-4),d1
|	rts

| Routines for doing fract arithmetic that check for overflow.
| In general, the variable frack_overflows is incremented each
| time an operation overflows

| Compute A*B, A,B fracts.
|  Overflow is carefully watched.  On overflow, the result is 0.
|	k = 2**16
|	A = Ah+Al/k		(Ah and Al are the two 16 bit halves)
|	B = Bh+Bl/k
| A*B	= (Ah+Al/k)*(Bh+Bl/k)
|	= Al*Bl/(k**2) + Ah*Bh + Bh*Al + Ah*Bl
_vfrmul:			| A*B	A at sp@(4)m B at sp@(8)
	subl	a0,a0
	tstl	sp@(4)		| A := abs(A)
	jge	L1v
	negl	sp@(4)
	addql	#1,a0
L1v:
	tstl	sp@(8)		| B := abs(B)
	jge	L2v
	negl	sp@(8)
	subql	#1,a0		| a0 == 0 iff A&B have the same sign.
L2v:
	movw	sp@(6),d0
	mulu	sp@(10),d0	| Al*Bl
	clrw	d0
	swap	d0		| top 16 of Al*Bl to bottom 16 of d0
	movw	sp@(4),d1
	mulu	sp@(8),d1	| Ah*Bh
	swap	d1
	jmi	overflow
	tstw	d1
	jne	overflow
	clrw	d1		| bottom 16 of Ah*Bh to top 16 of d1
	addl	d1,d0
	jvs	overflow
	movw	sp@(6),d1
	mulu	sp@(8),d1	| Al*Bh
	addl	d1,d0
	jvs	overflow
	movw	sp@(4),d1
	mulu	sp@(10),d1	| Ah*Bl
	addl	d1,d0
	jvs	overflow
	cmpl	#0,a0
	beq	RTS
	negl	d0
	jvs	overflow
	rts

overflow:
	addqw	#1,_fract_overflows
	clrl	d0
	rts

_vfradd:
	movl	sp@(8), d0
	addl	sp@(4), d0
	jvs	overflow
	rts

_vfrsub:
	movl	sp@(4),d0
	subl	sp@(8),d0
	jvs	overflow
	rts

	.data1
L25:
	.ascii	"%f\0"	
L26:
	.ascii	"%F\0"	

_fract_overflows:
	.word	0