dtof.a68
;
;	FLOATING POINT ERROR VALUES
;
ERANGE		equ	13
EDOM		equ	14
;
;	convert double to float
;
;	argument in d0/d1	result in d0
;

		dseg
		xref	_errno

		cseg
		public	.dtof
.dtof:	movem.l	d2/d3/d4,-(sp)
		move.l	d0,d2			;save a copy for exponent extraction
		beq		exit			;no work if zero
		smi		d4				;get sign of result
		swap	d2
		lsr.w	#4,d2
		and.w	#$7ff,d2		;extract exponent
		sub.w	#1023,d2		;unbias
		and.l	#$fffff,d0		;extract mantissa
		or.l	#$100000,d0		;turn hidden bit back on
		move.l	#2,d3
dagain:
		lsl.l	#1,d1			;shift mantissa into place
		roxl.l	#1,d0
		dbra	d3,dagain
		cmp.l	#$80000000,d1	;check for round
		beq		rstar
		bcs		nornd
round:
		add.l	#1,d0			;round up
		cmp.l	#$00ffffff,d0	;check for carry
		bls		nornd
		add.w	#1,d2			;bump exponent
		lsr.l	#1,d0			;and shift result back into place
		bra		nornd
rstar:
		or.w	#1,d0			;rstar case, force low bit to 1
nornd:
		add.w	#127,d2			;add in float bias
		bmi		underflow		;can't represent double in float
		cmp.w	#255,d2
		bgt		overflow		;ditto
		lsl.w	#7,d2			;back in place
		tst.b	d4				;need to set sign?
		beq		ndsign			;no
		or.w	#$8000,d2		;yes, set it
ndsign:
		swap	d2				;get exponent in high word
		clr.w	d2				;clear out
		and.l	#$7fffff,d0		;get rid of hidden bit
		or.l	d2,d0			; put back in number
		bra		exit
underflow:
		IF INT32
		move.l	#ERANGE,_errno	;under flow
		ELSE
		move.w	#ERANGE,_errno	;under flow
		ENDC
		move.l	#$00800000,d0	;set to smallest number
dosign:
		tst.b	d4
		beq		exit
		or.l	#$80000000,d0	;set sign
exit:
		movem.l	(sp)+,d2/d3/d4
		rts
;
overflow:
		IF INT32
		move.l	#ERANGE,_errno	;over flow
		ELSE
		move.w	#ERANGE,_errno	;over flow
		ENDC
		move.l	#$7fffffff,d0	;set to largest number
		tst.b	d4
		beq		1$
		or.l	#$80000000,d0	;set sign
1$
		ori		#2,ccr			;set overflow condition
		movem.l	(sp)+,d2/d3/d4
		rts
;
dtox.a68
; Copyright (C) 1989 by Manx Software Systems, Inc.
;:ts=8
;
; Convert double to extended
;
; arg in d0/d1, result in .p0
;

	dseg
	xref	.p0
	cseg


	public	.dtox
.dtox
	movem.l	d2/d3,-(sp)		;preserve registers
	move.l	d0,d2			;save sign and mantissa
	swap	d0			;get exponent for double
	and.w	#$7ff0,d0		;check if zero exponent
	beq	zero			;yes, so result is zero
	lsr.w	#4,d0			;extract exponent
	add.w	#(16383-1023),d0	;add in differences of biases
	tst.l	d2			;check if negative
	bpl	positive		;no, don't set sign bit
	or.w	#$8000,d0		;yes, set sign bit
positive:
	swap	d0			;put in high word for 96 bit format
	move.w	#0,d0			;clear low word
	move.w	#11,d3
	lsl.l	d3,d2			;extract mantissa
	or.l	#$80000000,d2		;set normalized bit
	rol.l	d3,d1			;shift low long word
	move.w	d1,d3			;save low word
	and.w	#$7ff,d3		;isolate bits
	or.w	d3,d2			;fill in 11 bits
	and.w	#$f800,d1		;isolate high bits
	exg	d1,d2			;put in the right order
exit:
	movem.l	d0-d2,.p0		;place results into extended
	movem.l	(sp)+,d2/d3		;restore registers
	rts

zero:
	move.l	#0,d0
	move.l	#0,d1
	move.l	#0,d2
	bra	exit

fabs.a68
; Copyright (C) 1985 by Manx Software Systems, Inc.
; :ts=8

	public	.Pabs
	public	_fabs
_fabs
	movem.l	4(sp),d0/d1
	jmp	.Pabs

floats.a68
; Copyright (C) 1986 by Manx Software Systems, Inc.
; :ts=8

	public	.Ffix
.Ffix:
	move.l	d1,-(sp)
	jsr	.ftod#
	jsr	.Pfix#
	movem.l	(sp)+,d1
	rts

	public	.Fflt
.Fflt:
	move.l	d1,-(sp)
	jsr	.Pflt#
	jsr	.dtof#
	movem.l	(sp)+,d1
	rts

	public	.Fneg
.Fneg:
	move.l	d1,-(sp)
	jsr	.ftod#
	jsr	.Pneg#
	jsr	.dtof#
	movem.l	(sp)+,d1
	rts

	public	.Ftst
.Ftst:
	movem.l	d0/d1,-(sp)
	jsr	.ftod#
	jsr	.Ptst#
	movem.l	(sp)+,d0/d1
	rts

	public	.Fadd
.Fadd:
	movem.l	d2/d3,-(sp)
	move.l	d1,d2
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Padd#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fcmp
.Fcmp:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Pcmp#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fdiv
.Fdiv:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Pdiv#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fmul
.Fmul:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Pmul#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fsub
.Fsub:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Psub#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

flt68.a68
; Copyright 1987 Manx Software Systems, Inc
;
;	FLOATING POINT ERRORS
;
ERANGE		equ	13
EDOM		equ	14

FAST_SHIFT	equ	1
NGUARD		equ	11
LOWMANT		equ	1<<NGUARD
ROUND		equ	LOWMANT>>1
STICKY		equ	1
;
INTOVER equ	2
;
		dseg
		xref	_errno
		cseg
;
;	fix a floating point argument
;
;	floating argument in d0/d1, integer result in d0
;
		public	.Pfix
.Pfix:	movem.l	d1/d2,-(sp)		; save registers
		move.l	d0,d2			; copy high order into d2
		smi.b	d1				; save sign in low 8 bits of d1
		swap	d2				; get exponent into low word
		lsr.w	#4,d2			; move into place
		and.w	#$7ff,d2		; extract exponent
		sub.w	#1023,d2		; unbias exponent
		bmi		zeroresult		; zero or negative exponent is zero
		sub.w	#31,d2			; exponent greater than 32 is overflow
		bgt		overflow
		neg.w	d2
		lsl.l	#8,d0			; get mantissa over to left
		lsl.l	#3,d0
		or.l	#$80000000,d0	; turn hidden bit back on
		ext.w	d1				; extend sign flag to a word value
		swap	d1				; get high order word of low order longword
;
		lsr.w	#5,d1
		and.w	#$7ff,d1		; get eleven bits;
		or.w	d1,d0			; insert into shifted mantissa
		lsr.l	d2,d0			; shift by remaining exponent
		move.l	d1,d1			; negative input?
		bmi		signed			; bif yes
tstfix:	tst.l	d0				; set condition code (probably a waste of time)
notest:	movem.l	(sp)+,d1/d2
		rts
signed:	neg.l	d0				; negate result
		bra		notest
zeroresult:
		move.l	#0,d0
		bra		notest
overflow:
		move.l	#$7fffffff,d0	; positive overflow
		move.l	d1,d1
		bpl		over1
		neg.l	d0				; negative overflow
over1:	ori		#INTOVER,ccr 	;set overflow bit
		bra		notest
;
;	float a fixed point argument
;
;	argument in d0, result in d0/d1
;
		public .Pflt
.Pflt:
		tst.l	d0				; value is zero ?
		beq		iszero			; make float zero
		movem.l	d4/d5/d6/d7,-(sp) ; (note: movem.l doesn't change cc)
		smi.b	d6				; determine sign (note: smi doesn't change cc)
		bgt		notneg			; test for negative value
		neg.l	d0				; get absolute value
notneg:
		move.l	#0,d1			; zero out LL of result
		move.l	#31,d4			; set exponent to end of HL
		jsr		donormnc		; go normalize without carry in
		movem.l	(sp)+,d4/d5/d6/d7
		rts
iszero:
		move.l	#0,d1			; answer is zero
		rts
;
;	compare two floating point values
;
;	arg1 in d0,d1  arg2 in d2,d3  answer in d0	trashes d2
;	-1 arg1 < arg2
;	 0 arg1 = arg2
;   +1 arg1 > arg2
;
		public .Pcmp
		IF 0
.Pcmp:	movem.l	d4/d5/d6/d7,-(sp)

		move.l	d2,d5
		smi.b	d7			; set sign of second
		swap	d5
		lsr.w	#4,d5		; extract exponent of second
		and.w	#$7ff,d5

		move.l	d0,d4
		smi.b	d6			; set sign of first
		swap	d4
		lsr.w	#4,d4		; extract exponent of first
		and.w	#$7ff,d4

		beq		firzero		; special code for first == 0
		tst.w	d5
		beq		seczero		; special code for second == 0
		cmp.b	d6,d7
		bne		difsign		; special code for differing signs
		move.l	#-1,d7		; assume first is less
		cmp.w	d5,d4		; compare exponents
		beq		cmpmant		; special code for exponents equal
		blt		cmpdone
		move.l	#1,d7		; second is less
		bra		cmpdone		; all except for sign check
firzero:
		move.b	d7,d6		; preserve sign of second
		move.l	#-1,d7		; assume first is less
		tst.w	d5			; check for other also zero
		bne		cmpdone		; all done except for sign check
		move.l	#0,d0
		bra		cmpexitz
seczero:
		move.l	#1,d7		; assume first is greater
		bra		cmpdone		; all done except for sign check
difsign:
		move.l	#-1,d0		; assume first is less
		tst.b	d6			; is first negative?
		bne		cmpexitz	; yes, all done
		move.l	#1,d0		; no, second is less because it's negative
		bra		cmpexitz	; all done
cmpmant:
		and.l	#$000fffff,d0	; extract mantissa of first
		;or.l	#$00100000,d0	; restore hidden bit
		and.l	#$000fffff,d2	; extract mantissa of second
		;or.l	#$00100000,d2	; restore hidden bit
		move.l	#-1,d7			; assume first is less
		cmp.l	d2,d0
		blt		cmpdone			; all done except for sign check
		beq		moretst			; go check low order
		move.l	#1,d7			; second is less
		bra		cmpdone			; go do sign check
moretst:
		cmp.l	d3,d1			; check low order
		bne		oktodo			; second is less
		move.l	#0,d0			; values are equal
		bra		cmpexitz
oktodo:
		bls		cmpdone			; first is less, go do sign check
		move.l	#1,d7			; second is less
cmpdone:
		tst.b	d6				; is sign set?
		beq		cmpexit			; no, don't have to flip 
		neg.l	d7				; flip
cmpexit:
		move.l	d7,d0			; set condition code
cmpexitz:
		movem.l	(sp)+,d4/d5/d6/d7
		rts
		ELSE

;	compare two floating point values
;
;	arg1 in d0,d1  arg2 in d2,d3  answer in d0
;	-1 arg1 < arg2
;	 0 arg1 = arg2
;   +1 arg1 > arg2
;
.Pcmp:
		movem.l	d1/d2/d3/d4,-(sp)
;		dc.w	%0100101011111100
		move.l	d0,d4			; test if arg1 is zero
		swap	d4
		and.w	#$7ff0,d4		; by checking if exponent is zero
		bne		cmpnz1			; bif not
		move.l	#0,d0			; insure a clean zero for arg1
		move.l	#0,d1
cmpnz1:
		move.l	d2,d4			; repeat test for arg2
		swap	d4
		and.w	#$7ff0,d4
		bne		cmpnz2
		move.l	#0,d2
		move.l	#0,d3
cmpnz2:
		move.l	d2,d2			; is second negative?
		bmi		cmpmi2			; bif so
		move.l	d0,d0			; 2nd is pos; is first negative?
		bmi		aLTb			; bif so
		cmp.l	d2,d0			; signed test of 1st longword of each arg
atstb:
		blt		aLTb
		bgt		aGTb
		cmp.l	d3,d1			; use an unsigned test for 2nd longword
		bcs		aLTb
		bhi		aGTb
aEQb:
		move.l	#0,d0			; arg1 == arg2
		movem.l	(sp)+,d1/d2/d3/d4
		rts
aLTb: 	move.l	#-1,d0
		movem.l	(sp)+,d1/d2/d3/d4
		rts
aGTb:	move.l	#1,d0
		movem.l	(sp)+,d1/d2/d3/d4
		rts
cmpmi2:							; 2nd is negative.
		move.l	d0,d0			; is 1st zero or positive?
		bpl		aGTb			; bif yes.
		exg		d1,d3
		cmp.l	d0,d2			; a bigger value means a smaller number
		bra		atstb
		ENDC
;
;	test a floating point number for zero
;
;	argument in d0,d1	result in d0	trashes d1
;
		public .Ptst
.Ptst:	move.l	d0,d1			; copy high order long word into d1
		swap	d1
		and.w	#$7ff0,d1		; test for zero exponent
		beq		tstzero			; argument is zero
		move.l	d0,d0			; test sign bit
		bpl		tstpos			; value is positive
		move.l	#-1,d0			; value is negative
		rts
tstpos:
		move.l	#1,d0
		rts

;
;	make a floating pooint value positive
;
;	argument and result are both in d0,d1
;
		public	.Pabs
.Pabs:	move.l	d0,-(sp)
		swap	d0
		and.w	#$7ff0,d0		; extract exponent
		beq		abszero			; zero exponent is zero
		move.l	(sp)+,d0		
		and.l	#$7fffffff,d0	; force sign bit off
		rts
abszero:
negzero:
		add.l	#4,sp			; pop stack
retzero:
		move.l	#0,d1
tstzero:
		move.l	#0,d0
		rts
;
;	negate a floating point argument
;
;	argument and result in d0,d1
;
		public	.Pneg
.Pneg:	move.l	d0,-(sp)
		swap	d0
		and.w	#$7ff0,d0		; extract exponent
		beq		negzero			; zero exponent is zero
		move.l	(sp)+,d0		; restore d0
		eor.l	#$80000000,d0	; flip sign bit
		rts
;
;	subtract/add two floating point numbers
;
;	arg1 in d0,d1i	arg2 in d2,d3	result in d0,d1	trashes d2,d3

addret2:							; return second operand
		move.l	d3,d1
		move.l	d2,d0
		swap	d2					; return true zero if second
		and.w	#$7ff0,d2			; operand exponent is zero.
		bne		addret1
divzero:							; return 0.0 = 0.0/x
		move.l	#0,d1
		move.l	#0,d0
;		FALLTHROUGH

addret1:							; return first operand
		movem.l	(sp)+,d4/d5/d6/d7/a0
		move.l	d0,d0				; (which is known not to be a funny zero)
		rts							; return with condition code set

		public	.Psub
.Psub:
		eor.l	#$80000000,d2		; flip sign bit of second operand
;		bra		.Padd				; go add values

		public	.Padd
.Padd:
		movem.l	d4/d5/d6/d7/a0,-(sp)
;		dc.w	%0100101011111100
		move.l	#0,d6
		move.l	#0,d7

		move.l	d0,d4		; break up first operand mantissa in d0,d1
		smi.b	d6			; sign in d6
		swap	d4
		and.w	#$7ff0,d4	; exponent in low word of d4.
		beq		addret2		; zero exponent is zero, result is second
		lsr.w	#4,d4		; low d4 := exponent of first operand

		move.l	d2,d5		; break up second operand mantissa in d2,d3
		smi.b	d7			; sign in d7
		swap	d5
		and.w	#$7ff0,d5	; exponent in low word of d5.
		beq		addret1		; an exponent of zero == 0.0, result is first
		lsr.w	#4,d5		; low d5 := exponent of second operand

		sub.w	#1023,d4	; adjust exponent of first
		sub.w	#1023,d5	; adjust exponent of second
		and.l	#$0fffff,d0	; extract mantissa of first
		or.l	#$100000,d0	; turn hidden bit back on

		move.l	d6,a0		; save d6
		IF FAST_SHIFT==0
		move.w	#(NGUARD-1),d6
l1:
		lsl.l	#1,d1		; now adjust mantissa to get guard digits
		roxl.l	#1,d0
		dbra	d6,l1
;
		and.l	#$0fffff,d2	; extract mantissa of second
		or.l	#$100000,d2	; turn hidden bit back on

		move.w	#(NGUARD-1),d6
l2:
		lsl.l	#1,d3		; adjust secondary for guard digits
		roxl.l	#1,d2
		dbra	d6,l2
;
		ELSE				; Fast left shift by 11
		swap	d5			; save second operand exponent in high d5
		move.w	#NGUARD,d5	; load the shift count for lsl.l instruction
		lsl.l	d5,d0		; adjust first mantissa to get guard bits
		move.l	d1,d6
		swap	d6
		lsr.w	#(16-NGUARD),d6
		or.w	d6,d0
		lsl.l	d5,d1		; have just shifted d0,d1 left by 11 bits
;
		and.l	#$0fffff,d2	; extract mantissa of second
		or.l	#$100000,d2	; turn hidden bit back on
;
		lsl.l	d5,d2		; adjust second mantissa to get guard bits
		move.l	d3,d6
		swap	d6
		lsr.w	#(16-NGUARD),d6
		or.w	d6,d2
		lsl.l	d5,d3		; have just shifted d2,d3 left by 11 bits
		swap	d5			; restore exponent
		ENDC
;
		move.l	a0,d6		; restore d6
;
		cmp.w	d5,d4		; find larger exponent
		beq		equdone		; don't need to equalize
		bgt		noswap		; first is larger, don"t swap
		exg.l	d0,d2		; swap for equalize
		exg.l	d1,d3
		exg.l	d4,d5
		exg.l	d6,d7
noswap:
		move.l	d4,a0		; save final exponent
		sub.w	d5,d4		; get equalization shift count
;		cmp.w	#53,d4		; going to shift it all away? *jd 9 Oct 86
		cmp.w	#55,d4
		bmi		doshft		; no, do the shift
		move.l	#0,d2		; yes, zero the secondary
		move.l	#0,d3
		bra		equdone0
doshft:
		IF FAST_SHIFT==1
;		dc.w	%0100101011111100
		cmp.w	#5,d4
		bmi		rotadd		; shift counts <=4 can be done faster with rotate.
		cmp.w	#32,d4
		bmi		shflt32		; fallthrough if right shifting more than 31 bits
		sub.w	#32,d4
		IFND STICKY
		move.l	d2,d3
		lsr.l	d4,d3
		move.l	#0,d2
		bra		equdone0
		ELSE
		move.l	d3,d3		; if any of the bits being lost are 1
		beq		sh32lowz
		move.l	d2,d3
		lsr.l	d4,d3
sh32sky:
		or.w	#STICKY,d3	; then or in STICKY bit.
		move.l	#0,d2
		bra		equdone0
sh32lowz:
		move.l	d2,d3		; 32 bit right shift of d2
		lsr.l	d4,d3		; d3 >> shift
		move.l	d3,d5
		lsl.l	d4,d5		; (d3 >> shift) << shift
		cmp.l	d2,d5		; if ((d3 >> shift) << shift) != d2,
		bne		sh32sky		; then bits will be lost during right shift of d3
		move.l	#0,d2
		bra		equdone0
		ENDC
shflt32:					; shift count < 32
		IFD STICKY
		;cmp.w	#NGUARD,d4	; if shifting <= number of guard bits,
		;ble		shflegu		; then no possibility of losing any bits.
		move.w	d6,-(sp)	; free up a register
		move.l	d3,d5
		lsr.l	d4,d3
		move.l	d3,d6
		lsl.l	d4,d6
		cmp.l	d5,d6
		beq		shfnone
		or.w	#STICKY,d3
shfnone:
		move.w	(sp)+,d6
		;bra		shflsr
shflegu:
		ENDC
		;lsr.l	d4,d3
shflsr:
		move.l	d2,d5
		lsr.l	d4,d2		d2 >>= count
		neg.w	d4
		add.w	#32,d4
		lsl.l	d4,d5		;d5 = first << (32 - count)
		or.l	d5,d3
		bra		equdone0
		ENDC
rotadd:						; Note that since shifting < 5 bits right, that
		sub.w	#1,d4		; no nonzero bits will drop off the right end and
shfttop:					; so no need to handle sticky bit.
		lsr.l	#1,d2		; equalize the secondary
		roxr.l	#1,d3
		dbra	d4,shfttop
equdone0:
		move.l	a0,d4		; restore exponent
equdone:
		cmp.w	d7,d6
		bgt		sub_2_1		; first is signed, sub first from second
		blt		sub_1_2		; second is signed, sub second from first
		add.l	d3,d1		; same sign, just add
		addx.l	d2,d0
		bcc		normit		; check for carry out on addition
		roxr.l	#1,d0		; carry, restore carried bit to result
		roxr.l	#1,d1
		add.w	#1,d4		; bump exponent
		bra		normit
sub_2_1:
		exg.l	d0,d2		; swap values
		exg.l	d1,d3
		move.l	d7,d6
sub_1_2:
		sub.l	d3,d1		; do subtraction
		subx.l	d2,d0
		bcc		normit		; did sign change (generated a borrow)?
		negx.b	d6			; yes, flip sign (-1 -> 0; 0 -> -1)
		neg.l	d1			; sign changed, subtract result from zero
		negx.l	d0
normit:
		jsr		donormnc	; go normalize without carry in
		movem.l	(sp)+,d4/d5/d6/d7/a0
		rts

multzero:
		move.l	#0,d0
		move.l	#0,d1
		movem.l	(sp)+,d4/d5/d6/d7/a0/a1/a2
		rts
;
;	multiply two floating point values
;
;	arg1 in d0,d1  arg2 in d2,d3 result in d0,d1
;
		public	.Pmul
.Pmul:	movem.l	d4/d5/d6/d7/a0/a1/a2,-(sp)
;		dc.w	%0100101011111100

		move.l	d0,d4		; break up first operand mantissa in d0,d1
		smi.b	d6			; sign in d6
		swap	d4			; exponent in d4
		and.w	#$7ff0,d4
		beq		multzero	; zero exponent is zero, result is zero
;
		move.l	d2,d5		; break up second operand mantissa in d2,d3
		smi.b	d7			; sign in d7
		swap	d5			; exponent in d5
		and.w	#$7ff0,d5
		beq		multzero	; exponent of zero is zero, result is zero
;
		and.l	#$0fffff,d0	; extract mantissa of first
		or.l	#$100000,d0	; turn hidden bit back on

		;lsr.w	#4,d4
		;lsr.w	#4,d5
		;sub.w	#1023,d4	; adjust exponent of first
		;sub.w	#1023,d5	; adjust exponent of second
		;add.w	d5,d4		; add exponents
		;add.w	#1,d4		; adjust for result of multiply
		add.w	d5,d4		; add exponents
		lsr.w	#4,d4
		sub.w	#(1023+1023-1),d4 ; unbias result exponent
		move.l	d4,a0		; save result exponent

		eor.b	d7,d6		; figure out sign of result
		move.l	d6,a1		; save sign

		IF FAST_SHIFT==0
		move.w	#(NGUARD-1),d6
m1:
		lsl.l	#1,d1		; now adjust mantissa to get guard digits
		roxl.l	#1,d0
		dbra	d6,m1
;
		and.l	#$0fffff,d2	; extract mantissa of second
		or.l	#$100000,d2	; turn hidden bit back on
		move.w	#(NGUARD-1),d6
m2:
		lsl.l	#1,d3		; adjust secondary for guard digits
		roxl.l	#1,d2
		dbra	d6,m2
		ELSE
		move.l	#NGUARD,d5	; load the shift count for lsl.l instruction
;
		lsl.l	d5,d0		; adjust first mantissa to get guard bits
		move.l	d1,d6
		swap	d6
		lsr.w	#(16-NGUARD),d6
		or.w	d6,d0
		lsl.l	d5,d1		; have just shifted d0,d1 left by 11 bits
;
		and.l	#$0fffff,d2	; extract mantissa of second
		or.l	#$100000,d2	; turn hidden bit back on
;
		lsl.l	d5,d2
		move.l	d3,d6
		swap	d6
		lsr.w	#(16-NGUARD),d6
		or.w	d6,d2
		lsl.l	d5,d3		; have just shifted d2,d3 left by 11 bits
		ENDC
;
		move.l	#0,d7		; make a zero for carry addition and donorm
							; because 68000 architecture is brain-damaged
		move.l	d2,a2		; save HL of multiplier for later
		move.l	d3,d2		; get LL of mmultiplier into multiplier reg
		move.l	#0,d3		; clear temporary accumulator
		move.l	#0,d4		; clear result registers	
		move.l	#0,d5
		move.l	#0,d6
		tst.w	d2			; check for zero LLLW multiplier
		beq		skip1		; don't bother, it's zero
		move.w	d0,d3		; check HLLW multiplicand
		beq		skip11		; don't bother it's zero
		mulu.w	d2,d3
		swap	d3			; get high order word of result
		move.w	d3,d6		; save in LLHW of final result 
skip11:
		move.l	d0,d3		; check HLHW multiplicand
		swap	d3
		tst.w	d3
		beq		skip1		; don't bother it's zero
		mulu.w	d2,d3
		add.l	d3,d6		; save in LL of final result
skip1:
		swap	d2			; check for zero LLHW multiplier
		tst.w	d2
		beq		skip2		; don't bother, it's zero
		move.l	d1,d3		; get LLHW of multiplicand
		swap	d3
		tst.w	d3			; check for zero
		beq		skip2x		; don't bother it's zero
		mulu.w	d2,d3
		clr.w	d3			; zero LW of result
		swap	d3			; get HW of result
		add.l	d3,d6		; add to LL of final result
		addx.w	d7,d5		; carry into HLLW of final result
skip2x:
		move.w	d0,d3		; check for zero HLLW
		beq		skip21		; don't bother, it's zero
		mulu.w	d2,d3
		add.l	d3,d6		; add to LL of final result
		addx.w	d7,d5		; carry into MLLW of final result
skip21:
		move.l	d0,d3		; get HLHW of multiplicand
		swap	d3

		mulu.w	d2,d3
		swap	d6
		add.w	d3,d6		; add to LLHW of final result
		swap	d6			; put LLHW back in place
		clr.w	d3			; zero low part of result
		swap	d3			; get high part of result
		addx.l	d3,d5		; add with carry into ML of final result
skip2:
		move.l	a2,d2		; get back HL of multiplier
		tst.w	d2			; check for zero HLLW of multiplier
		beq		skip3		; don't bother, it's zero
		move.w	d1,d3		; get LLLW of multiplicand
		beq		skip3x
		mulu.w 	d2,d3
		clr.w	d3
		swap	d3			; get HW of result
		add.l	d3,d6		; add into LL of final result
		addx.l	d7,d5		; carry into ML of final result
skip3x:
		move.l	d1,d3
		swap	d3
		tst.w	d3			; check for zero
		beq		skip31		; don't bother, it's zero
		mulu.w	d2,d3
		add.l	d3,d6		; add to LL of final result
		addx.l	d7,d5		; carry into ML of final result
skip31:
		move.w	d0,d3		; get HLLW of multiplicand
		beq		skip32		; don't bother, it's zero
		mulu.w	d2,d3
		swap	d6
		add.w	d3,d6		; add LW of result to LLHW of final result
		swap	d6			; put LLHW back into place
		clr.w	d3			; get rid of LW of result
		swap	d3			; get HW of result
		addx.l	d3,d5		; carry into ML of final result
skip32:
		move.l	d0,d3		; get HLHW of multiplicand
		swap	d3
							; don't need usual test for zero, can't happen
		mulu.w	d2,d3
		add.l	d3,d5		; add to ML of final result
		addx.l	d7,d4		; carry into HL of final result
skip3:
		swap	d2			; get HLHW of multiplier
							; don't need usual test for zero, can't happen
		move.w	d1,d3		; get LLLW of multiplicand
		beq		skip41
		mulu.w	d2,d3
		add.l	d3,d6		; add to LL of final result
		addx.l	d7,d5		; carry into ML of final result
		addx.l	d7,d4		; carry into HL of final result
skip41:
		swap	d1			; get LLHW of multiplicand
		move.w	d1,d3
		beq		skip42		; if zero, skip
		mulu.w	d2,d3
		swap	d6			; get LLHW of final result
		add.w	d3,d6		; add LW of result into LLHW of final result
		swap	d6			; get LLHW of final result back in place
		clr.w	d3			; get rid of LW of result
		swap	d3			; get high word of result
		;add.l	d3,d5		; add into ML of final result
		addx.l	d3,d5		; add into ML of final result + carry from prev add
		addx.l	d7,d4		; carry into HL of final result
skip42:
		move.w	d0,d3		; get HLLW of multiplicand
		beq		skip43		; don't bother, it's zero
		mulu.w	d2,d3
		add.l	d3,d5		; add to ML of final result
		addx.l	d7,d4		; carry into HL of final result
skip43:
		swap	d0			; get HLHW of multiplicand
		move.w	d0,d3
							; don't need check for zero, can't happen
		mulu.w	d2,d3
		swap	d5			; get MLHW of final result
		add.w	d3,d5		; add LW of result to it
		swap	d5			; put MLHW of final result back into place
		clr.w	d3			; get rid of LW of result
		swap	d3			; get HW of result
		addx.l	d3,d4		; add to HL of final result
		cmp.l	#$ffff,d4	; check for carry into HLHW of result
		bls		skip4
		add.w	#1,a1		; adjust expoonent
		lsr.l	#1,d4		; shift mantissa over one
		roxr.l	#1,d5
		roxr.l	#1,d6
skip4:
		cmp.w	#$8000,d6	; check for round
		beq		rstar		; this is the rstar special case
		bls		noround
mulround:
		swap	d6			; for round, get LLHW of result
		add.w	#1,d6		; round it up
		swap	d6			; put it back in place
		addx.l	d7,d5		; propagate carry
		addx.l	d7,d4		; ditto
		cmp.l	#$ffff,d4	; check for overflow out of HLLW
		bls		noround
		add.w	#1,a0		; bump exponent
		lsr.l	#1,d4		; shift it down
		roxr.l	#1,d5
		roxr.l	#1,d6
		bra 	noround
rstar:
		;or.l	#10000,d6	; force low order bit to one
		or.l	#$10000,d6	; force low order bit to one
; IEEE standard says that round to nearest is the default rounding mode; the
; representable value nearest to the infinitely precise result shall be
; delivered;  if the two nearest representable values are equally near, the one
; with its least significant bit zero shall be delivered.
		;move.l	d6,d1
		;and.l	#$10000,d1
		;bne		mulround
noround:
		move.w	d5,d6		; construct result
		swap	d6			;  d4  |  d5  |  d6  |
		move.l	d6,d1		;ov |   result   |rnd|
		move.w	d4,d5
		swap	d5
		move.l	d5,d0
		move.l	a0,d4		; move expoonent back into d4
		move.l	a1,d6		; sign into d6
		jsr	donorm			; normalize	
		movem.l	(sp)+,d4/d5/d6/d7/a0/a1/a2
		rts
;
;	divide two floating point values
;
;	arg1 in d0,d1 arg2 in d2,d3 result in d0,d1
;	trashes d2,d3
;
		public	.Pdiv
.Pdiv	movem.l	d4/d5/d6/d7/a0,-(sp)

		move.l	d0,d4		; break up first operand mantissa in d0,d1
		smi.b	d6			; sign in d6
		swap	d4			; exponent in d4
		and.w	#$7ff0,d4
		beq		divzero		; zero exponent is zero, result is zero
;
		move.l	d2,d5		; break up second operand mantissa in d2,d3
		smi.b	d7			; sign in d7
		swap	d5			; exponent in d5
		and.w	#$7ff0,d5
		beq		divbyzero	; exponent of zero is zero, result is fault
;
		and.l	#$0fffff,d0	; extract mantissa of first
		or.l	#$100000,d0	; turn hidden bit back on

		;lsr.w	#4,d4
		;lsr.w	#4,d5
		;sub.w	#1023,d4	; adjust exponent of first
		;sub.w	#1023,d5	; adjust exponent of second
		;sub.w	d5,d4		; subtract exponents
		;sub.w	#1,d4		; adjust for result of divide
		sub.w	d5,d4		; subtract exponents
		asr.w	#4,d4
		sub.w	#1,d4		; adjust for result of divide
		move.w	d4,a0		; sign extend result exponent and save it
		eor.b	d7,d6		; figure out sign of result
		swap	d6			; save sign in high half of d6

		IF FAST_SHIFT==0
		move.w	#(NGUARD-1),d7
n1:
		lsl.l	#1,d1		; now adjust mantissa to get guard digits
		roxl.l	#1,d0
		dbra	d7,n1
;
		and.l	#$0fffff,d2	; extract mantissa of second
		or.l	#$100000,d2	; turn hidden bit back on
		move.w	#(NGUARD-1),d7
n2:
		lsl.l	#1,d3		; adjust secondary for guard digits
		roxl.l	#1,d2
		dbra	d7,n2
		ELSE
		move.l	#NGUARD,d5	; load the shift count for lsl.l instruction
;
		lsl.l	d5,d0		; adjust first mantissa to get guard bits
		move.l	d1,d7
		swap	d7
		lsr.w	#(16-NGUARD),d7
		or.w	d7,d0
		lsl.l	d5,d1		; have just shifted d0,d1 left by 11 bits
;
		and.l	#$0fffff,d2	; extract mantissa of second
		or.l	#$100000,d2	; turn hidden bit back on
;
		lsl.l	d5,d2
		move.l	d3,d7
		swap	d7
		lsr.w	#(16-NGUARD),d7
		or.w	d7,d2
		lsl.l	d5,d3		; have just shifted d2,d3 left by 11 bits
		ENDC
;
;		check for easy divide case
;
		tst.w	d2			; check for zero HLLW divisor
		bne		harddiv
		tst.l	d3			; check for zero LL divisor
		bne		harddiv
;
		swap	d2			; get HLHW into place as divisor
		swap	d0
		cmp.w	d2,d0		; make sure divide won't overflow
		beq		okdiv
		bls		ediv
okdiv:
		swap	d0			; get dividend back in place
		add.w	#1,a0		; bump exponent
		lsr.l	#1,d0		; make dividend smaller
		roxr.l	#1,d1
		swap	d0			; make nop of next instruction
ediv:
		swap	d0			; get dividend back in place
		divu.w	d2,d0
		move.w	d0,d4		; save result as HLHW of final result
		swap	d4			
		swap	d1			; get next word of dividend
		move.w	d1,d0		; move it in with remainder of last calc
		divu.w	d2,d0
		move.w	d0,d4		; save result as HLLW of final result
		swap	d1			; get next word of dividend
		move.w	d1,d0
		divu.w	d2,d0
		move.w	d0,d5		; save result as LLHW of result
		swap	d5
		clr.w	d0			; get rid of quotient
		divu.w	d2,d0
		move.w	d0,d5		; save result as LLLW of final result
		move.l	d4,d0		; get into place for normalization
		move.l	d5,d1
		bra 	divnorm		; go normalize result
harddiv:
		cmp.l	d0,d2
		bne		do_div
		cmp.l	d1,d3
		bne		do_div
		move.l	#$80000000,d0	; mantissas same, answer is one with exponent
		add.w	#1,a0
		move.l	#0,d1
		bra		divnorm		; go normalize and exit
do_div:
		bhi		noadjust	;		*jd 9 Oct 86
		add.w	#1,a0		; bump exponent
		lsr.l	#1,d0		; shift dividend
		roxr.l	#1,d1
noadjust:
		lsr.l	#1,d2
		roxr.l	#1,d3
		lsr.l	#1,d0
		roxr.l	#1,d1
		move.w	#2,d6		; outer loop count (don't use move.l since high
							; half of d6 contains result sign flag)
bdiv_loop:
		move.l	#0,d4		; clear result
		move.l	#31,d7		; inner loop count
div_loop:
		lsl.l	#1,d4		; shift result
		lsl.l	#1,d1		; shift dividend
		roxl.l	#1,d0
		sub.l	d3,d1		; subtract divisor from dividend
		subx.l	d2,d0
		bmi		zero_bit
one_bit:
		;add.l	#1,d4		; set bit in quotient
		add.w	#1,d4		; set bit in quotient
		dbra	d7,div_loop
		sub.w	#1,d6
		beq		divend
		move.l	d4,d5
		bra		bdiv_loop
;
bzero_loop:
		move.l	d4,d5
		move.l	#0,d4
		move.l	#31,d7		; inner loop count
zero_loop:
		lsl.l	#1,d4		; shift zero quotient bit in
		lsl.l	#1,d1		; shift dividend
		roxl.l	#1,d0
		add.l	d3,d1		; add divisor back in
		addx.l	d2,d0
		bpl		one_bit
zero_bit:
		dbra	d7,zero_loop
		sub.w	#1,d6
		bne		bzero_loop
;
divend:
		move.l	d5,d0		; move result into d0,d1
		move.l	d4,d1
divnorm:
		move.l	a0,d4		; get exponent back
		swap	d6			; get sign back
		jsr		donormnc	; go normalize without carry in
divret:
		movem.l	(sp)+,d4/d5/d6/d7/a0
		rts
divbyzero:
		IF INT32
		move.l	#EDOM,_errno
		ELSE
		move.w	#EDOM,_errno
		ENDC
		jsr		setinf
		bra		divret

normzero:
		move.l	#0,d1				; d0 is already known to contain 0
		rts

normover:
		IF INT32
		move.l	#ERANGE,_errno
		ELSE
		move.w	#ERANGE,_errno
		ENDC
setinf:	move.l	#$ffffffff,d1
		move.l	#$7fffffff,d0
		tst.b	d6
		beq		setover
		or.l	#$80000000,d0
setover:
		ori		#INTOVER,ccr
		rts
;
;		normalize a floating point value and pack it
;
;		mantissa in d0,d1, exponent in d4, sign in d6, carry in d7
;		result in d0,d1; trashes d4,d5,d7
;		condition code will be set correctly.

donormnc:
		move.l	#0,d7		; make a zero for donorm
donorm:
		tst.l	d0			; test for HL is zero
		bne		normnext	; is not, go on
		exg.l	d0,d1		; move up LL
		sub.w	#32,d4		; adjust exponent
		tst.l	d0			; is it zero too
		beq		normzero	; yes, result is zero
normnext:
		bmi		nonorm		; is sign bit on (i.e. number normalized)?
normnn:
		sub.w	#1,d4		; no, adjust exponent
		lsl.l	#1,d1		; and shift number
		roxl.l	#1,d0
		;bpl		normnext	; sign on now? no, keep looping
		;beq		normnext
		bge		normnn		; loop until sign bit is set
nonorm:						; normalized, now check for round
		move.w	d1,d5		; get low order word
		and.w	#$7ff,d5	; extract eleven bits
		cmp.w	#ROUND,d5	; compare to .5
		beq		nrstar		; special rstar case
		bls		normfin		; round down if less
normrnd:
		add.l	#LOWMANT,d1	; round up
		addx.l	d7,d0		; and carry into HL
		bcc		normfin		; check for overflow
		roxr.l	#1,d0		; and re-normalize if necessary
		roxr.l	#1,d1
		add.w	#1,d4		; bump exponent
		bra		normfin
nrstar:
		;or.l	#LOWMANT,d1	; rstar round
; IEEE standard says that round to nearest is the default rounding mode; the
; representable value nearest to the infinitely precise result shall be
; delivered;  if the two nearest representable values are equally near, the one
; with its least significant bit zero shall be delivered.
		move.w	d1,d5
		and.w	#LOWMANT,d5	; if last bit is on,
		bne		normrnd		; then round up
							; else round down
normfin:
		IF FAST_SHIFT==0
		move.l	#(NGUARD-1),d5
packloop:
		lsr.l	#1,d0		; put mantissa back in place
		roxr.l	#1,d1
		dbra	d5,packloop
		ELSE
		move.l	#NGUARD,d5
		move.w	d0,d7
		lsl.w	#(16-NGUARD),d7
		lsr.l	d5,d0
		lsr.l	d5,d1
		swap	d1
		or.w	d7,d1
		swap	d1
		ENDC
		and.l	#$0fffff,d0	; get rid of hidden bit
		add.w	#1023,d4	; add bias back into exponent
		bmi		normunder	; check for underflow case
		cmp.w	#2047,d4	; check for overflow case
		bgt		normover
		tst.b	d6			; check for sign
		beq		normnsign	
		or.w	#$800,d4	; turn sign bit on
normnsign:
		lsl.w	#4,d4
		swap	d0
		or.w	d4,d0		; or in the exponent
		swap	d0
		rts
normunder:
		move.l	#0,d1
		IF INT32
		move.l	#ERANGE,_errno
		ELSE
		move.w	#ERANGE,_errno
		ENDC
		move.l	#$00100000,d0	; smallest number
		tst.b	d6
		beq		exit
		or.l	#$80000000,d0
exit:
		rts
		end
frexp.a68
;
;	FLOATING POINT ERROR VALUES
;
UNDER_FLOW	equ	1
OVER_FLOW	equ	2
DIV_BY_ZERO	equ	3
;
;	frexp(d, &i)
;
;	returns 1/2 <= |x| < 1
;	such that d = x *2^i
;
		public	_frexp
_frexp:	move.l	d2,-(sp)
		move.l	16(sp),a0		;get address for int
		move.l	8(sp),d0		; get double value
		move.l	12(sp),d1
		move.l	d0,d2
		swap	d2
		and.w	#$7ff0,d2
		bne		notzero
		move.l	#0,d2
		beq		done
notzero:
		and.l	#$800fffff,d0	;get rid of old exponent
		lsr.w	#4,d2
		sub.w	#1022,d2
		or.l	#$3fe00000,d0	;change exponent to -1
done:
		if INT32
		ext.l	d2
		move.l	d2,(a0)
		else
		move.w	d2,(a0)
		endc
		move.l	(sp)+,d2
		rts
ftoa.c
/* Copyright (C) 1984 by Manx Software Systems, Inc. */

static double round[] = { 10, 1,
	5e-1, 5e-2, 5e-3, 5e-4, 5e-5, 5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
	5e-11, 5e-12, 5e-13, 5e-14, 5e-15, 5e-16 };

void
ftoa(double number, register char *buffer, int maxwidth, int mode)
{
	register int i;
	register int exp, digit, decpos, ndig;
	register int flag;

	if ((*(short *)&number&0x7ff0) == 0x7ff0) {
		exp = (*(short *)&number&0x8000) ? '-' : '+';
		while (maxwidth--)
			*buffer++ = exp;
		*buffer = 0;
		return;
	}

	/*
	 *	Convert the number to a number between 1 and 10 or 0. Calculate the
	 *	exponent at the same time.
	 */

	exp = 0;
	if (number < 0.0) {
		number = -number;
		*buffer++ = '-';
	}
	if (number > 0.0) {
		while (number < round[1]) {
			number *= round[0];
			--exp;
		}
		while (number >= round[0]) {
			number /= round[0];
			++exp;
		}
	}


	/*
	 *	The number of digits to display is equal to the number of digits
	 *	following the decimal point (maxwidth) + digits to the left (exp) +
	 *	a guaranteed digit since the number 10 > x >= 1. For E format, the
	 *	exp digits are not counted, while G format works differently.
	 */

	flag = mode & 3;
	if (flag == 2) {						/* 'g' format */
		if (maxwidth == 0)
			maxwidth = 1;
		if (exp < -4 || exp >= maxwidth)
			flag = -1;						/* switch to 'e' format */
		ndig = maxwidth;
	} else if (flag == 1)					/* 'f' format */
		ndig = maxwidth + exp + 1;
	else									/* 'e' format */
		ndig = maxwidth + 1;

	/*
	 *	Round the number up using 5e-x where x is one more than the
	 *	precision of the number.
	 */

	if (ndig > 0) {
		if ((number += round[(ndig>16?16:ndig)+1]) >= round[0]) {
			number = round[1];
			++exp;
			if (flag > 0)
				++ndig;
		}
	}

	if (flag > 0) {
		if (exp < 0) {
			*buffer++ = '0';
			*buffer++ = '.';
			i = -exp - 1;
			if (ndig <= 0)
				i = maxwidth;
			while (i--)
				*buffer++ = '0';
			decpos = 0;
		} else {
			decpos = exp+1;
		}
	} else {
		decpos = 1;
	}

	if (ndig > 0) {
		for (i = 0 ; ; ++i) {
			if (i < 16) {
				digit = (int)number;
				*buffer++ = digit+'0';
				number = (number - digit) * round[0];
			} else
				*buffer++ = '0';
			if (--ndig == 0)
				break;
			if (decpos && --decpos == 0)
				*buffer++ = '.';
		}
	}

	if (decpos && (mode&0x20))
		*buffer++ = '.';
	if (flag <= 0) {
		*buffer++ = (mode&0x10)?'E':'e';
		if (exp < 0) {
			exp = -exp;
			*buffer++ = '-';
		} else
			*buffer++ = '+';
		*buffer++ = exp/100 + '0';
		exp %= 100;
		*buffer++ = exp/10 + '0';
		*buffer++ = exp%10 + '0';
	}
	if (decpos == 0 && (flag == 2 || flag == -1) && (mode&0x20) == 0) {
		while (*--buffer == '0')
			;
		if (*buffer != '.')
			++buffer;
	}
	*buffer = 0;
}

ftod.a68
;
;	convert float to double
;
;	argument in d0 result in d0/d1
;
		public	.ftod
.ftod:	movem.l	d3/d4,-(sp)
		tst.l	d0				;save sign
		smi		d3
		clr.l	d1				;clear out LL of result
		and.l	#$7fffffff,d0	;extract exponent and mantissa
		move.l	#2,d4
fagain:
		lsr.l	#1,d0			;shift exponent and mantissa into position
		roxr.l	#1,d1
		dbra	d4,fagain
		move.l	d0,d4
		or.l	d1,d4			;check for 0
		beq		1$
		add.l	#$38000000,d0	;adjust bias for double (1023-127) << 20
1$:
		tst.b	d3				; have to set sign?
		beq		fexit			; no
		or.l	#$80000000,d0	;yes, set it
fexit:
		movem.l	(sp)+,d3/d4
		rts
ftox.a68
; Copyright (C) 1989 by Manx Software Systems, Inc.
;:ts=8
;
; Convert float to extended
;
; argument in d0, result in .p0
;

	dseg
	xref	.p0

	cseg

	public	.ftox
.ftox
	movem.l	d1/d2,-(sp)		;preserve registers
	move.l	d0,d1			;save for mantissa work
	smi.b	d2			;save sign
	swap	d0			;get high word
	and.w	#$7f80,d0		;check if zero exponent
	beq	zero			;yes, so result is zero
	lsr.w	#7,d0			;extract exponent
	add.w	#(16383-127),d0		;add in differences of biases
	tst.b	d2			;check if negative
	beq	positive		;no, don't set sign bit
	or.w	#$8000,d0		;yes, set sign bit
positive:
	swap	d0			;put in high word for 96 bit format
	move.w	#0,d0			;clear low word
	lsl.l	#8,d1			;extract mantissa
	or.l	#$80000000,d1		;set normalized bit
exit:
	move.l	#0,d2			;low long word always zero
	movem.l	d0-d2,.p0		;place results into extended
	movem.l	(sp)+,d1/d2		;restore registers
	rts

zero:
	move.l	#0,d0
	move.l	#0,d1
	bra	exit

ldexp.a68
;
;	FLOATING POINT ERROR VALUES
;
UNDER_FLOW	equ	1
OVER_FLOW	equ	2
DIV_BY_ZERO	equ	3
;
ERANGE	equ	13
EDOM	equ	14
;
		xdef	_errno
;
;		ldexp(d, i)
;
;		return x = d * 2^i
;
		public	_ldexp
_ldexp	move.l	d2,-(sp)
		movem.l	8(sp),d0/d1		;get the number
		move.l	d0,d2			;get exponent
		and.l	#$800fffff,d0	;get old exponent out
		swap	d2
		and.l	#$7ff0,d2
		lsr.w	#4,d2
		if INT32
		add.w	18(sp),d2
		else
		add.w	16(sp),d2
		endc
		ble		ldunder
		cmp.l	#2047,d2
		bgt		ldover
		lsl.w	#4,d2
		swap	d2
		or.l	d2,d0			;put new exponent in
		move.l	(sp)+,d2
		rts
ldunder:
		move.w	#UNDER_FLOW,_flterr
		if INT32
		move.l	#ERANGE,_errno
		else
		move.w	#ERANGE,_errno
		endc
;		or.l	#$00100000,d0
		move.l	#$c3a99999,d0	;-HUGE_VAL
		move.l	#$8d2c7175,d1
		move.l	(sp)+,d2
		rts
ldover:
		move.w	#OVER_FLOW,_flterr
		if INT32
		move.l	#ERANGE,_errno
		else
		move.w	#ERANGE,_errno
		endc
;		or.l	#$7ff00000,d0
		move.l	#$43a99999,d0
		move.l	#$8d2c7175,d1
		move.l	(sp)+,d2
		rts
;
		global	_flterr,2
makefile
CC=qcc
AS=as
CFLAGS=
AFLAGS=
O=oss

.c.$O:
	$(CC) -o $@ $(CFLAGS) $*.c
.a68.$O:
	$(AS) -o $@ $(AFLAGS) $*.a68

ALIB=\
	dtof.$O dtox.$O fabs.$O floats.$O flt68.$O \
	frexp.$O ftod.$O ftox.$O ldexp.$O modf.$O \
	regs.$O x80to96.$O x96to80.$O xflt68.$O xtod.$O \
	xtof.$O

mlib:	$(ALIB)

modf.a68
;
;	FLOATING POINT ERROR VALUES
;
UNDER_FLOW	equ	1
OVER_FLOW	equ	2
DIV_BY_ZERO	equ	3
;
;		modf(d, dptr)
;
;		return fractional part of d and
;		stores integral part in	*dptr
;
		public	.Pfix
		public	.Pflt
		public	.Psub
		public	_modf
_modf:
		movem.l	d2/d3,-(sp)
		move.l	12(sp),d0		;pick up double
		move.l	16(sp),d1
		move.l	d0,d2			;copy number
		move.l	d1,d3			;
		jsr		.Pfix			;fix the number
;		bvc		fits			;fits in long
		cmp.l	#$7fffffff,d0	;positive overflow
		beq		over1
		cmp.l	#$80000001,d0	;negative overflow
		bne		fits
over1:
		move.l	d2,d0			;overflow
		move.l	d3,d1			;so return original
		bra		over2
fits:
		jsr		.Pflt			;float it again
over2:
		move.l	20(sp),a0
		move.l	d0,(a0)+		;store integral part
		move.l	d1,(a0)
		exg.l	d0,d2
		exg.l	d1,d3
		jsr		.Psub			;get rid of integral part
		movem.l	(sp)+,d2/d3
		rts
regs.a68
	global	.p0,12
	global	.p1,12

x80to96.a68
; Copyright (C) 1989 by Manx Software Systems, Inc.
;:ts=8
;
; Convert 80 bit extended to 96 bit extended: arg1 = arg2
; Declared as: void _X80to96(extended96 *arg1, extended *arg2)
; Return address = (sp)
; arg1 = 4(sp)
; arg2 = 8(sp)
;
	public	_X80to96
_X80to96
	movem.l	a0/a1,-(sp)	;preserve a0/a1
	move.l	12(sp),a0	;get address of extended96, arg1
	move.l	16(sp),a1	;get address of extended, arg2
	move.w	(a1)+,(a0)+	;copy 16 bit exponent
	move.w	#0,(a0)+	;put in 16 bit zero for padding
	move.l	(a1)+,(a0)+	;copy first 32 bits of mantissa
	move.l	(a1)+,(a0)+	;copy last 32 bits of mantissa
	movem.l	(sp)+,a0/a1	;restore a0/a1
	rts
;
x96to80.a68
; Copyright (C) 1989 by Manx Software Systems, Inc.
;:ts=8
;
; Convert 96 bit extended to 80 bit extended: arg1 = arg2
; Declared as: void _X96to80(extended *arg1, extended96 *arg2)
; Return address = (sp)
; arg1 = 4(sp)
; arg2 = 8(sp)
;
	public	_X96to80
_X96to80
	movem.l	a0/a1,-(sp)	;preserve a0/a1
	move.l	12(sp),a0	;get address of extended, arg1
	move.l	16(sp),a1	;get address of extended96, arg2
	move.w	(a1)+,(a0)+	;copy 16 bit exponent
	add.w	#2,a1		;skip zero bytes
	move.l	(a1)+,(a0)+	;copy first 32 bits of mantissa
	move.l	(a1),(a0)	;copy last 32 bits of mantissa
	movem.l	(sp)+,a0/a1	;restore a0/a1
	rts
;
xflt68.a68
; Copyright 1989 Manx Software Systems, Inc
;
;	FLOATING POINT ERRORS
;
ERANGE		equ	13
EDOM		equ	14

FAST_SHIFT	equ	1
STICKY		equ	1
;
INTOVER equ	2
;
		dseg
		xref	_errno
		xref	.p0,.p1


		cseg
;
;
; Convert extended to long integer
;
; argument in .p0, result in d0
;

		public	.Xfix
.Xfix:
		move.l	d1,-(sp)
		movem.l	.p0,d0-d1			; put extended in registers
		swap	d0					; get exponent in low word
		ext.l	d0					; shift sign into high word
		and.w	#$7fff,d0			; extract exponent
		sub.w	#$3fff,d0			; unbias exponent
		bmi		zeroresult			; zero or negative exponent is zero
		sub.w	#31,d0				; exponent greater than 32 is overflow
		bgt		overflow
		neg.w	d0					; change sign to positive
		lsr.l	d0,d1				; (divide) shift by remaining exponent
		move.l	d0,d0				; negative input?
		bmi		signed
exitfix:
		move.l	d1,d0				; copy in result
notest:
		movem.l	(sp)+,d1			; restore registers
		rts
signed:
		neg.l	d1					; negate result
		bra		exitfix
zeroresult:
		move.l	#0,d1
		bra		exitfix
overflow:
		move.l	#$7fffffff,d1		; positive overflow
		move.l	d0,d0				; negative input?
		bpl		over1				; bif no
		neg.l	d1					; negative overflow
		move.l	d1,d0				; copy in long result
over1:
		ori		#INTOVER,ccr		; set overflow bit
		bra		notest


;
; Convert long integer to extended
;
; argument in d0, result in .p0
;

	public	.Xflt
.Xflt
		movem.l	d1-d3,-(sp)			; save registers
		move.l	d0,d1				; put long in register
		tst.l	d1					; value is zero ?
		beq		iszero				; make extended zero
		movem.l	d6/d7,-(sp)			; (note: movem.l doesn't change cc)
		smi.b	d6					; determine sign (note: smi doesn't change cc)
		bgt		notneg				; test for negative value
		neg.l	d1					; get absoulte value
notneg:
		move.l	#0,d2				; zero out LL of result
		move.l	#0,d3				; zero out guard digits
		move.l	#31,d0				; set exponent
		jsr	donormnc				; go normalize without carry in
		movem.l	(sp)+,d6/d7
exitflt:
		movem.l	(sp)+,d1-d3		; restore registers
		rts
iszero:
		move.l	#0,.p0			; answer is zero
		move.l	#0,.p0+4
		move.l	#0,.p0+8
		bra		exitflt


;
; Extended compare: (arg1 - arg2)
;
; arg1 in .p0, arg2 in .p1
;
; returns answer in d0
;	-1; arg1 < arg2
;	 0; arg1 = arg2
;   +1; arg1 > arg2
;
		public	.Xcmp
.Xcmp:
		movem.l	d1-d5,-(sp)
		move.l	8*4+4(sp),a0
		movem.l	.p0,d0-d2		; test if arg1 is zero
		and.l	#$7fff0000,d0	; by checking if exponent is zero
		bne		cmpnz1			; bif not
		move.l	#0,d1			; insure a clean zero for arg1
		move.l	#0,d2
cmpnz1:
		movem.l	.p1,d3-d5			; test if arg2 is zero
		and.l	#$7fff0000,d3	; by checking if exponent is zero
		bne		cmpnz2			; bif not
		move.l	#0,d4			; insure a clean zero for arg2
		move.l	#0,d5
cmpnz2:
		move.l	.p1,d3			; is second negative
		bmi		cmpmi2			; bif so
		move.l	.p0,d0			; 2nd is pos; is first negative?
		bmi		aLTb			; bif so
		cmp.l	d3,d0			; signed test of 1st longword of each arg
atstb:
		blt		aLTb
		bgt		aGTb
		cmp.l	d4,d1			; use an unsigned test for 2nd longword
		bcs		aLTb
		bhi		aGTb
		cmp.l	d5,d2			; use an unsigned test for 3rd longword
		bcs		aLTb
		bhi		aGTb
aEQb:
		move.l	#0,d0			; *arg1 == *arg2
		movem.l	(sp)+,d1-d5
		rts
aLTb:
		move.l	#-1,d0			; *arg1 < *arg2
		movem.l	(sp)+,d1-d5
		rts
aGTb:
		move.l	#1,d0			; *arg1 > *arg2
		movem.l	(sp)+,d1-d5
		rts
cmpmi2:							; 2nd is negative.
		move.l	.p0,d0			; is 1st zero or positive
		bpl		aGTb			; bif yes
		exg.l	d1,d4
		exg.l	d2,d5
		cmp.l	d0,d3
		bra		atstb


;
; Extended test: (arg1 - 0L)
;
; argument in .p0, result in d0
;
	public	.Xtst
.Xtst
		move.l	.p0,d0			; test if zero
		and.l	#$7fff0000,d0	; by checking if exponent is zero
		beq		tstzero			; argument is zero
		move.w	.p0,d0
		bpl		tstpos			; value is positive
		move.l	#-1,d0			; value is negative
tstzero:
		rts
tstpos:
		move.l	#1,d0
		rts


;
; Get absolute value of extended: arg1 = |arg1|
;
; argument in .p0, result in .p0
;
	public	.Xabs
.Xabs
		movem.l	d0,-(sp)
		move.w	.p0,d0			; test if zero
		and.w	#$7fff,d0		; by checking if exponent is zero
		beq		abszero			; argument is zero
		and.w	#$7fff,.p0		; force sign bit off
		movem.l	(sp)+,d0
		rts
abszero:
negzero:
		move.l	#0,.p0
		move.l	#0,.p0+4
		move.l	#0,.p0+8
		movem.l	(sp)+,d0
		rts


;
; Negate an extended: arg1 = -arg1
;
; argument in .p0, result in .p0
;
	public	.Xneg
.Xneg
		movem.l	d0,-(sp)
		move.w	.p0,d0			; test if zero
		and.w	#$7fff,d0		; by checking if exponent is zero
		beq		negzero			; argument is zero
		eor.w	#$8000,.p0		; flip sign bit
		movem.l	(sp)+,d0
		rts


;
; subtract/add two extended floating point numbers
; 1st in .p0
; 2st in .p1
;
addret2:
		movem.l	d3-d5,.p0
		and.l	#$7fff0000,d3
		bne		addret1
divzero:
		move.l	#0,.p0
		move.l	#0,.p0+4
		move.l	#0,.p0+8
addret1:
		move.l	.p0,d0				; return with condition code set
		movem.l	(sp)+,a2/a3/d0-d7
		rts

;
; Extended subtract: .p0 = .p0 - .p1
;
		public	.Xsub
.Xsub
		eor.w	#$8000,.p1		; flip sign bit of second operand
; Fall through to .Xadd


;
; Extended add: .p0 = .p0 + .p1
;
		public	.Xadd
.Xadd
		movem.l	a2/a3/d0-d7,-(sp)
		movem.l	.p0,d0-d2			; load arg1
		movem.l	.p1,d3-d5			; load arg2
		move.l	#0,d6
		move.l	#0,d7
		move.l	d0,d0
		smi.b	d6					; sign for arg1 in d6
		swap	d0
		and.w	#$7fff,d0
		beq		addret2				; zero exponent is zero, result is second
		move.l	d3,d3
		smi.b	d7					; sign for arg2 in d7
		swap	d3
		and.w	#$7fff,d3
		beq		addret1				; zero exponent is zero, result is first
		sub.w	#16383,d0			; adjust exponent of arg1
		sub.w	#16383,d3			; adjust exponent of arg2
		cmp.w	d3,d0				; find larger exponent
		bne		1$
		move.w	#0,d3
		bra		equdone				; don't need to equalize
1$
		bgt		noswap				; first is larger, don't swap
		exg.l	d0,d3				; swap to equalize
		exg.l	d1,d4
		exg.l	d2,d5
		exg.l	d6,d7
noswap:
		move.l	d0,a2				; save final exponent
		IFD STICKY
		movem.w	d6/d7,-(sp)			; free up d6/d7
		move.l	#0,d6
		ELSE
		move.w	d7,-(sp)			; free up d7
		ENDC
		sub.w	d3,d0				; get equalization shift count
		move.w	#0,d3				; zero out guard bits
		cmp.w	#80,d0				; going to shift it all away?
		bmi		doshft				; no, do the shift
		move.l	#0,d4				; yes, zero the secondary
		move.l	#0,d5
		bra		equdone0
doshft:
		IF FAST_SHIFT==1
		cmp.w	#5,d0
		bmi		rotadd		; shift counts <5 can be done faster with rotate.
		cmp.w	#64,d0
		bmi		shflt64		; fallthrough if right shifting more than 63 bits
		sub.w	#64,d0
		IFND STICKY
		swap	d4
		move.w	d4,d3
		lsr.w	d0,d3
		move.l	#0,d4
		move.l	#0,d5
		bra		equdone0
		ELSE
		move.w	d4,d4		; if any of the bits being lost are 1
		bne		1$
		tst.l	d5
		beq		sh64lowz
1$
		swap	d4
		move.w	d4,d3
		lsr.w	d0,d3
sh64sky:
		or.w	#STICKY,d3	; then or in STICKY bit.
		move.l	#0,d4
		move.l	#0,d5
		bra		equdone0
sh64lowz:
		swap	d4
		move.w	d4,d3		; 64 bit right shift of d4
		lsr.l	d0,d3		; d3 >> shift
		move.l	d3,d7
		lsl.l	d0,d7		; (d7 >> shift) << shift
		cmp.w	d4,d7		; if ((d7 >> shift) << shift) != d4,
		bne		sh64sky		; then bits will be lost during right shift of d7
		move.l	#0,d4
		move.l	#0,d5
		bra		equdone0
		ENDC
shflt64:					; shift count < 64
		cmp.w	#48,d0
		bmi		shflt48		; fallthrough if right shifting more than 47 bits
		IFD		STICKY
		tst.l	d5			; if any of the bits being lost are 1
		beq		1$
		move.l	#1,d6
1$
		ENDC
		move.w	d4,d3
		move.l	#0,d5
		swap	d4
		move.w	d4,d5
		move.l	#0,d4
		sub.w	#48,d0
		bra		shflt16
shflt48:					; shift count < 48
		cmp.w	#32,d0
		bmi		shflt32		; fallthrough if right shifting more than 31 bits
		swap	d5
		move.w	d5,d3		; if any of the bits being lost are 1
		IFD		STICKY
		beq		1$
		move.l	#1,d6
1$
		ENDC
		move.l	d4,d5
		move.l	#0,d4
		sub.w	#32,d0
		bra		shflt16
shflt32:					; shift count < 32
		cmp.w	#16,d0
		bmi		shflt16		; fallthrough if right shifting more than 15 bits
		move.w	d5,d3		; if any of the bits being lost are 1
		move.w	d4,d5
		swap	d5
		move.w	#0,d4
		swap	d4
		sub.w	#16,d0
shflt16
		beq		chksky
		move.w	d3,d7
		lsr.w	d0,d7
		lsl.w	d0,d7
		cmp.w	d3,d7
		beq		1$
		move.l	#1,d6
1$
		cmp.w	#5,d0
		bmi		rotadd		; shift counts <5 can be done faster with rotate.
		swap	d6			; preserve sticky bit
		move.w	#16,d6
		sub.w	d0,d6
		lsr.w	d0,d3
		move.w	d5,d7
		lsl.w	d6,d7
		or.w	d7,d3
		move.w	#32,d6
		sub.w	d0,d6
		lsr.l	d0,d5
		move.l	d4,d7
		lsl.l	d6,d7
		or.l	d7,d5
		lsr.l	d0,d4
		swap	d6			; restore sticky bit
		bra		chksky
		ENDC
rotadd:
		sub.w	#1,d0
shfttop:
		lsr.l	#1,d4		; equalize the secondary
		roxr.l	#1,d5
		roxr.w	#1,d3
		dbra	d0,shfttop
		IFD STICKY
chksky
		tst.l	d6
		beq		equdone0
		or.w	#STICKY,d3
		ENDC
equdone0:
		move.l	a2,d0		; restore exponent
		IFD STICKY
		movem.w	(sp)+,d6/d7	; restore d6/d7
		ELSE
		move.w	(sp)+,d7	; restore d7
		ENDC
equdone:
		cmp.w	d7,d6
		bgt		sub_2_1		; first is signed, sub first from second
		blt		sub_1_2		; second is signed, sub second from first
		add.l	d5,d2		; same sign, just add
		addx.l	d4,d1
		bcc		normit		; check for carry out on addition
		roxr.l	#1,d1		; carry, restore carried bit to result
		roxr.l	#1,d2
		roxr.w	#1,d3
		add.w	#1,d0		; bump exponent
		bra		normit
sub_2_1:
		exg.l	d1,d4		; swap values
		exg.l	d2,d5
		move.l	d7,d6
		sub.l	d5,d2
		bra		sub_1_2a
sub_1_2:
		neg.w	d3			; do subtraction
		subx.l	d5,d2
sub_1_2a:
		subx.l	d4,d1
		bcc		normit		; did sign change (generated a borrow)?
		negx.b	d6			; yes, flip sign (-1 -> 0; 0 -> -1)
		neg.w	d3			; sign changed, subtract result from zero
		negx.l	d2
		negx.l	d1
normit:
		jsr		donormnc	; go normalize without carry in
		movem.l	(sp)+,a2/a3/d0-d7
		rts


multzero:
		move.l	#0,.p0
		move.l	#0,.p0+4
		move.l	#0,.p0+8
		movem.l	(sp)+,a1-a3/d0-d7
		rts
;
; Extended multiply: .p0 = .p0 * .p1
;
		public	.Xmul
.Xmul
		movem.l	a1-a3/d0-d7,-(sp)
		movem.l	.p0,d0-d2			; load arg1
		movem.l	.p1,d3-d5			; load arg2
		move.l	d0,d0
		smi.b	d6					; sign for arg1 in d6
		swap	d0					; exponent in d0
		and.w	#$7fff,d0
		beq		multzero			; zero exponent is zero, result is zero
		move.l	d3,d3
		smi.b	d7					; sign for arg2 in d7
		swap	d3					; exponent in d3
		and.w	#$7fff,d3
		beq		multzero			; zero exponent is zero, result is zero
		add.w	d3,d0				; add exponents
		sub.w	#(16383+16383-1),d0	; unbias result exponent
		move.w	d0,a2				; save result exponent
		eor.b	d7,d6				; figure out sign of result
		move.w	d6,a3				; save sign
; don't need to adjust mantissa since it fits exactly
		move.l	#0,d7				; make a zero for carry addition and donorm
									; because 68000 architecture is brain-damaged
		move.l	d4,a1				; save HL of multiplier for later
		move.l	d5,d4				; get LL of mmultiplier into multiplier reg
		move.l	#0,d5				; clear temporary accumulator
		move.l	#0,d0				; clear result registers	
		move.l	#0,d3
		move.l	#0,d6
		tst.w	d4					; check for zero LLLW multiplier
		beq		skip1				; don't bother, it's zero
		move.w	d1,d5				; check HLLW multiplicand
		beq		skip11				; don't bother it's zero
		mulu.w	d4,d5
		swap	d5					; get high order word of result
		move.w	d5,d6				; save in LLHW of final result 
skip11:
		move.l	d1,d5				; check HLHW multiplicand
		swap	d5
		tst.w	d5
		beq		skip1				; don't bother it's zero
		mulu.w	d4,d5
		add.l	d5,d6				; save in LL of final result
skip1:
		swap	d4					; check for zero LLHW multiplier
		tst.w	d4
		beq		skip2				; don't bother, it's zero
		move.l	d2,d5				; get LLHW of multiplicand
		swap	d5
		tst.w	d5					; check for zero
		beq		skip2x				; don't bother it's zero
		mulu.w	d4,d5
		clr.w	d5					; zero LW of result
		swap	d5					; get HW of result
		add.l	d5,d6				; add to LL of final result
		addx.w	d7,d3				; carry into HLLW of final result
skip2x:
		move.w	d1,d5				; check for zero HLLW
		beq		skip21				; don't bother, it's zero
		mulu.w	d4,d5
		add.l	d5,d6				; add to LL of final result
		addx.w	d7,d3				; carry into MLLW of final result
skip21:
		move.l	d1,d5				; get HLHW of multiplicand
		swap	d5

		mulu.w	d4,d5
		swap	d6
		add.w	d5,d6				; add to LLHW of final result
		swap	d6					; put LLHW back in place
		clr.w	d5					; zero low part of result
		swap	d5					; get high part of result
		addx.l	d5,d3				; add with carry into ML of final result
skip2:
		move.l	a1,d4				; get back HL of multiplier
		tst.w	d4					; check for zero HLLW of multiplier
		beq		skip3				; don't bother, it's zero
		move.w	d2,d5				; get LLLW of multiplicand
		beq		skip3x
		mulu.w 	d4,d5
		clr.w	d5
		swap	d5					; get HW of result
		add.l	d5,d6				; add into LL of final result
		addx.l	d7,d3				; carry into ML of final result
skip3x:
		move.l	d2,d5
		swap	d5
		tst.w	d5					; check for zero
		beq		skip31				; don't bother, it's zero
		mulu.w	d4,d5
		add.l	d5,d6				; add to LL of final result
		addx.l	d7,d3				; carry into ML of final result
skip31:
		move.w	d1,d5				; get HLLW of multiplicand
		beq		skip32				; don't bother, it's zero
		mulu.w	d4,d5
		swap	d6
		add.w	d5,d6				; add LW of result to LLHW of final result
		swap	d6					; put LLHW back into place
		clr.w	d5					; get rid of LW of result
		swap	d5					; get HW of result
		addx.l	d5,d3				; carry into ML of final result
skip32:
		move.l	d1,d5				; get HLHW of multiplicand
		swap	d5
									; don't need usual test for zero, can't happen
		mulu.w	d4,d5
		add.l	d5,d3				; add to ML of final result
		addx.l	d7,d0				; carry into HL of final result
skip3:
		swap	d4					; get HLHW of multiplier
									; don't need usual test for zero, can't happen
		move.w	d2,d5				; get LLLW of multiplicand
		beq		skip41
		mulu.w	d4,d5
		add.l	d5,d6				; add to LL of final result
		addx.l	d7,d3				; carry into ML of final result
		addx.l	d7,d0				; carry into HL of final result
skip41:
		swap	d2					; get LLHW of multiplicand
		move.w	d2,d5
		beq		skip42				; if zero, skip
		mulu.w	d4,d5
		swap	d6					; get LLHW of final result
		add.w	d5,d6				; add LW of result into LLHW of final result
		swap	d6					; get LLHW of final result back in place
		clr.w	d5					; get rid of LW of result
		swap	d5					; get high word of result
		addx.l	d5,d3				; add into ML of final result + carry from prev add
		addx.l	d7,d0				; carry into HL of final result
skip42:
		move.w	d1,d5				; get HLLW of multiplicand
		beq		skip43				; don't bother, it's zero
		mulu.w	d4,d5
		add.l	d5,d3				; add to ML of final result
		addx.l	d7,d0				; carry into HL of final result
skip43:
		swap	d1					; get HLHW of multiplicand
		move.w	d1,d5
									; don't need check for zero, can't happen
		mulu.w	d4,d5
		swap	d3					; get MLHW of final result
		add.w	d5,d3				; add LW of result to it
		swap	d3					; put MLHW of final result back into place
		clr.w	d5					; get rid of LW of result
		swap	d5					; get HW of result
		addx.l	d5,d0				; add to HL of final result
		cmp.l	#$ffff,d0			; check for carry into HLHW of result
		bls		skip4
		add.w	#1,a2				; adjust expoonent
		lsr.l	#1,d0				; shift mantissa over one
		roxr.l	#1,d3
		roxr.l	#1,d6
skip4:
		move.w	d6,d5
		move.w	d3,d6				; construct result
		swap	d6					;  d0  |  d3  |  d6  |
		move.l	d6,d2				;ov |   result   |rnd|
		move.w	d0,d3
		swap	d3
		move.l	d3,d1
		move.w	d5,d3
		move.w	a2,d0				; move expoonent back into d0
		move.w	a3,d6				; sign into d6
		jsr	donorm					; normalize	
		movem.l	(sp)+,a1-a3/d0-d7
		rts


;
; Extended division: .p0 = .p0 / .p1
;
		public	.Xdiv
.Xdiv
		movem.l	a2/a3/d0-d7,-(sp)
		movem.l	.p0,d0-d2			; load arg1
		movem.l	.p1,d3-d5			; load arg2
		move.l	d0,d0
		smi.b	d6					; sign for arg1 in d6
		swap	d0					; exponent in d0
		and.w	#$7fff,d0
		beq		divzero				; zero exponent is zero, result is zero
		move.l	d3,d3
		smi.b	d7					; sign for arg2 in d7
		swap	d3					; exponent in d3
		and.w	#$7fff,d3
		beq		divbyzero			; zero exponent is zero, result is fault
		sub.w	d3,d0				; subtract exponents
		sub.w	#1,d0				; adjust for result of divide
		move.l	d0,a2				; save result exponent
		eor.b	d7,d6				; figure out sign of result
		swap	d6					; save sign in high half of d6
; don't need to adjust mantissa since it fits exactly
;
;		check for easy divide case
;
		move.l	#0,d3		; zero out guard bits
		tst.w	d4			; check for zero HLLW divisor
		bne		harddiv
		tst.l	d5			; check for zero LL divisor
		bne		harddiv
;
		swap	d4			; get HLHW into place as divisor
		swap	d1
		cmp.w	d4,d1		; make sure divide won't overflow
		beq		okdiv
		bls		ediv
okdiv:
		swap	d1			; get dividend back in place
		add.w	#1,a2		; bump exponent
		lsr.l	#1,d1		; make dividend smaller
		roxr.l	#1,d2
		roxr.w	#1,d3
		swap	d1			; make nop of next instruction
ediv:
		swap	d1			; get dividend back in place
		divu.w	d4,d1
		move.w	d1,d0		; save result as HLHW of final result
		swap	d0			
		swap	d2			; get next word of dividend
		move.w	d2,d1		; move it in with remainder of last calc
		divu.w	d4,d1
		move.w	d1,d0		; save result as HLLW of final result
		swap	d2			; get next word of dividend
		move.w	d2,d1
		divu.w	d4,d1
		move.w	d1,d7		; save result as LLHW of result
		swap	d7
;		clr.w	d1			; get rid of quotient
		move.w	d3,d1
		divu.w	d4,d1
		move.w	d1,d7		; save result as LLLW of final result
		clr.w	d1			; get rid of quotient
		divu.w	d4,d1
		move.w	d1,d3		; save result as guard bits
		move.l	d0,d1		; get into place for normalization
		move.l	d7,d2
		bra 	divnorm		; go normalize result
harddiv:
		cmp.l	d1,d4
		bne		do_div
		cmp.l	d2,d5
		bne		do_div
		move.l	#$80000000,d1	; mantissas same, answer is one with exponent
		add.w	#1,a2
		move.l	#0,d2
		bra		divnorm		; go normalize and exit
do_div:
		bhi		noadjust	;		*jd 9 Oct 86
		add.w	#1,a2		; bump exponent
		lsr.l	#1,d1		; shift dividend
		roxr.l	#1,d2
		roxr.w	#1,d3
noadjust:
		lsr.l	#1,d4
		roxr.l	#1,d5
		lsr.l	#1,d1
		roxr.l	#1,d2
		roxr.w	#1,d3
		move.w	#2,d6		; outer loop count (don't use move.l since high
							; half of d6 contains result sign flag)
bdiv_loop:
		move.l	#0,d0		; clear result
		move.l	#31,d7		; inner loop count
div_loop:
		lsl.l	#1,d0		; shift result
		lsl.w	#1,d3		; shift dividend
		roxl.l	#1,d2
		roxl.l	#1,d1
		sub.l	d5,d2		; subtract divisor from dividend
		subx.l	d4,d1
		bmi		zero_bit
one_bit:
		;add.l	#1,d0		; set bit in quotient
		add.w	#1,d0		; set bit in quotient
		dbra	d7,div_loop
		sub.w	#1,d6
		beq		divguard
		move.l	d0,a3
		bra		bdiv_loop
;
bzero_loop:
		move.l	d0,a3
		move.l	#0,d0
		move.l	#31,d7		; inner loop count
zero_loop:
		lsl.l	#1,d0		; shift zero quotient bit in
		lsl.w	#1,d3		; shift dividend
		roxl.l	#1,d2
		roxl.l	#1,d1
		add.l	d5,d2		; add divisor back in
		addx.l	d4,d1
		bpl		one_bit
zero_bit:
		dbra	d7,zero_loop
		sub.w	#1,d6
		bne		bzero_loop
;
divguard:
; calculate guard digits
		move.w	#0,d3		; clear result
		move.l	#15,d7		; inner loop count
g_div_loop:
		lsl.l	#1,d3		; shift result
		lsl.l	#1,d2		; shift dividend
		roxl.l	#1,d1
		sub.l	d5,d2		; subtract divisor from dividend
		subx.l	d4,d1
		bmi		g_zero_bit
g_one_bit:
		add.w	#1,d3		; set bit in quotient
		dbra	d7,g_div_loop
		bra		divend
;
g_zero_loop:
		lsl.l	#1,d3		; shift zero quotient bit in
		lsl.l	#1,d2		; shift dividend
		roxl.l	#1,d1
		add.l	d5,d2		; add divisor back in
		addx.l	d4,d1
		bpl		g_one_bit
g_zero_bit:
		dbra	d7,g_zero_loop
;
divend:
		move.l	a3,d1		; move result into d0,d1
		move.l	d0,d2
divnorm:
		move.l	a2,d0		; get exponent back
		swap	d6			; get sign back
		jsr		donormnc	; go normalize without carry in
divret:
		movem.l	(sp)+,a2/a3/d0-d7
		rts
divbyzero:
		IF INT32
		move.l	#EDOM,_errno
		ELSE
		move.w	#EDOM,_errno
		ENDC
		jsr		setinf
		bra		divret


normzero:
		move.l	#0,d0				; d1,d2 are already known to contain 0
		bra		exitnorm

normover:
		IF INT32
		move.l	#ERANGE,_errno
		ELSE
		move.w	#ERANGE,_errno
		ENDC
setinf:	move.l	#$ffffffff,d2
		move.l	#$ffffffff,d1
		move.l	#$7fff0000,d0
		tst.b	d6
		beq		setover
		or.l	#$80000000,d0
setover:
		ori		#INTOVER,ccr
		bra		exitnorm
;
;		normalize an extended floating point value and pack it
;
;		mantissa in d1,d2, exponent in d0, sign in d6, carry in d7
;		guard bits in low word of d3
;		store result in .p0
;		trashes d0-d2,d7
;		condition code will be set correctly.

donormnc:
		move.l	#0,d7		; make a zero for donorm
donorm:
		tst.l	d1			; test for HL is zero
		bne		normnext	; is not, go on
		exg.l	d1,d2		; move up LL
		exg.l	d3,d2		; move up guard bits
		swap	d2
		move.w	#0,d2
		sub.w	#32,d0		; adjust exponent for zero in HL
		tst.l	d1			; is LL zero too
		bne		normnext	; is not, go on
		exg.l	d1,d2		; check guard bits
		sub.w	#32,d0		; adjust exponent for zero in LL
		tst.l	d1			; are guard bits zero too
		beq		normzero	; yes, result is zero
normnext:
		bmi		nonorm		; is sign bit on (i.e. number normalized)?
normnn:
		sub.w	#1,d0		; no, adjust exponent
		lsl.w	#1,d3		; and shift number
		roxl.l	#1,d2
		roxl.l	#1,d1
		;bpl		normnext	; sign on now? no, keep looping
		;beq		normnext
		bge		normnn		; loop until sign bit is set
nonorm:						; normalized, now check for round
		cmp.w	#$8000,d3			; check for round
		beq		rstar				; this is the rstar special case
		bls		noround
normrnd:
		add.l	#1,d2
		addx.l	d7,d1
		bcc		noround		; check for overflow
		roxr.l	#1,d1		; and re-normalize if necessary
		roxr.l	#1,d2
		add.w	#1,d0		; bump exponent
		bra 	noround
rstar:
;		or.b	#1,d2				; force low order bit to one
; IEEE standard says that round to nearest is the default rounding mode; the
; representable value nearest to the infinitely precise result shall be
; delivered;  if the two nearest representable values are equally near, the one
; with its least significant bit zero shall be delivered.
		btst	#0,d2
		bne		normrnd
noround:
; Don't have to readjust mantissa by NGUARD digits since exact fit
;		or.l	#$80000000,d1	;normalized bit should already be set
		add.w	#16383,d0	; add bias back into exponent
		bmi		normunder	; check for underflow case
		cmp.w	#32767,d0	; check for overflow case
		bgt		normover
		tst.b	d6			; check for sign
		beq		normnsign	
		or.w	#$8000,d0	; turn sign bit on
normnsign:
		and.l	#$ffff,d0
		swap	d0
exitnorm:
		movem.l	d0-d2,.p0	; store exponent & mantissa
		rts
normunder:
		move.l	#$80000000,d1	;set normalized bit
		move.l	#0,d2
		IF INT32
		move.l	#ERANGE,_errno
		ELSE
		move.w	#ERANGE,_errno
		ENDC
		move.l	#$00010000,d0	; smallest number
		tst.b	d6
		beq		exitnorm
		or.l	#$80000000,d0
		bra		exitnorm
		end
xtod.a68
; Copyright (C) 1989 by Manx Software Systems, Inc.
;:ts=8
;
;	FLOATING POINT ERROR VALUES
;
ERANGE	equ	13
EDOM	equ	14
;

	dseg
	xref	_errno
	xref	.p0

	cseg
;
; Convert extended to double
;
; argument in .p0, result in d0/d1
;

	public	.xtod
.xtod
	movem.l	d2-d5,-(sp)
	movem.l	.p0+4,d0-d1		;get extended mantissa
	move.w	.p0,d2			;get extended exponent
	smi	d3			;get sign of result
	and.w	#$7fff,d2		;extract exponent
	beq	zero			;no work if zero
	sub.w	#16383,d2		;unbias
	move.w	d1,d5			;get low order word
	and.w	#$7ff,d5		;extract 11 bits
	cmp.w	#$400,d1		;check for round
	beq	rstar
	bcs	nornd
round:
	add.l	#$800,d1		;round up
	bcc	1$
	add.l	#1,d0
1$
	bcc	nornd			;check for overflow
	roxr.l	#1,d0			;and re-normalize if necessary
	roxr.l	#1,d1
	add.w	#1,d2			;bump exponent
	bra	nornd
rstar:
;	or.l	#$800,d1		;rstar case, force low bit to 1
; IEEE standard way of handling rstar case
	move.w	d1,d5
	and.w	#$800,d5		;if last bit is on
	bne	round			;then round up
nornd:					;else round down
	move.l	#11,d5			;pack mantissa for double
	move.w	d0,d4
	lsl.w	#5,d4
	lsr.l	d5,d0
	lsr.l	d5,d1
	swap	d1
	or.w	d4,d1
	swap	d1
	and.l	#$fffff,d0		;get rid of hidden bit
	add.w	#1023,d2		;add in double bias
	bmi	underflow		;check for underflow
	cmp.w	#2047,d2		;check for overflow
	bgt	overflow
	tst.b	d3			;need to set sign?
	beq	nosign			;no
	or.w	#$800,d2		;yes, set it
nosign:
	lsl.w	#4,d2			;back in place
	swap	d0
	or.w	d2,d0
	swap	d0
	bra	exit
underflow:
	move.l	#0,d1
	IF INT32
	move.l	#ERANGE,_errno		;under flow
	ELSE
	move.w	#ERANGE,_errno		;under flow
	ENDC
	move.l	#$00100000,d0		;set to smallest number
dosign:
	tst.b	d3
	beq	exit
	or.l	#$80000000,d0		;set sign
exit:
	movem.l	(sp)+,d2-d5
	rts
overflow:
	IF INT32
	move.l	#ERANGE,_errno		;over flow
	ELSE
	move.w	#ERANGE,_errno		;over flow
	ENDC
	move.l	#$ffffffff,d1		;set to largest number
	move.l	#$7fffffff,d0
	tst.b	d3
	beq	exit
	or.l	#$80000000,d0		;set sign
	ori	#2,ccr			;set overflow condition
	bra	exit
zero:
	move.l	#0,d0
	move.l	#0,d1
	bra	exit
;
xtof.a68
; Copyright (C) 1989 by Manx Software Systems, Inc.
;:ts=8
;
;	FLOATING POINT ERROR VALUES
;
ERANGE	equ	13
EDOM	equ	14
;

	dseg
	xref	_errno
	xref	.p0

	cseg
;
; Convert extended to float
;
; argument in .p0, result in d0
;

	public	.xtof
.xtof
	movem.l	d1-d4,-(sp)
	movem.l	.p0+4,d0-d1		;get extended mantissa
	move.w	.p0,d2			;get extended exponent
	smi	d4			;get sign of result
	and.w	#$7fff,d2		;extract exponent
	beq	zero			;no work if zero
	sub.w	#16383,d2		;unbias
	lsr.l	#8,d1
	move.b	d0,d3
	lsl.w	#8,d3
	swap	d3
	move.w	#0,d3
	or.l	d3,d1
	lsr.l	#8,d0
	cmp.l	#$80000000,d1		;check for round
	beq	rstar
	bcs	nornd
round:
	add.l	#1,d0			;round up
	cmp.l	#$00ffffff,d0		;check for carry
	bls	nornd
	add.w	#1,d2			;bump exponent
	lsr.l	#1,d0			;and shift result back into place
	bra	nornd
rstar:
;	or.w	#1,d0			;rstar case, force low bit to 1
; IEEE standard way of handling rstar case
	btst	#0,d0			;if last bit is on
	bne	round
nornd:
	add.w	#127,d2			;add in float bias
	bmi	underflow		;can't represent extended in float
	cmp.w	#255,d2
	bgt	overflow		;ditto
	lsl.w	#7,d2			;back in place
	tst.b	d4			;need to set sign?
	beq	ndsign			;no
	or.w	#$8000,d2		;yes, set it
ndsign:
	swap	d2			;get exponent in high word
	clr.w	d2			;clear out
	and.l	#$7fffff,d0		;get rid of hidden bit
	or.l	d2,d0			; put back in number
	bra	exit
underflow:
	IF INT32
	move.l	#ERANGE,_errno		;under flow
	ELSE
	move.w	#ERANGE,_errno		;under flow
	ENDC
	move.l	#$00800000,d0		;set to smallest number
dosign:
	tst.b	d4
	beq	exit
	or.l	#$80000000,d0		;set sign
exit:
	movem.l	(sp)+,d1-d4
	rts
overflow:
	IF INT32
	move.l	#ERANGE,_errno		;over flow
	ELSE
	move.w	#ERANGE,_errno		;over flow
	ENDC
	move.l	#$7fffffff,d0		;set to largest number
	tst.b	d4
	beq	1$
	or.l	#$80000000,d0		;set sign
1$
	ori	#2,ccr			;set overflow condition
	movem.l	(sp)+,d1-d4
	rts
zero:
	move.l	#0,d0
	bra	exit
;