/* Copyright (C) 1995 Eddie C. Dost
This file is part of the HP48 C Library.

The HP48 C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public License
as published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

The HP48 C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details. */

#include <rom.h>

	xdef ___adddf3
	xdef ___subdf3
	xdef ___muldf3
	xdef ___divdf3

	xdef ___trunctfdf2
	xdef ___extenddftf2

	xdef ___fixdfsi
	xdef ___fixdfdi
	xdef ___fixdfti

	xdef ___fixunsdfsi
	xdef ___fixunsdfdi
	xdef ___fixunsdfti

	xdef ___eqdf2
	xdef ___nedf2
	xdef ___gtdf2
	xdef ___gedf2
	xdef ___ltdf2
	xdef ___ledf2

	xdef ___ckinf		; used by math library
	xdef ___packsb		; used by math library
	xdef ___pack		; used by math library

___ckinf
	;
	; Check for infinite result:
	;
	; if (XM != 0)
	;   if (P == 3)
	;     Infinite Result
	;   else
	;     Undefined Result
	;
	rteq	#1,hst		; return if XM == 0
	beq.1	#3,p,___ckinf2
	; error: Undefined Result
	rtn

___ckinf2
	; error Infinite Result
	move.1	#0,p
	move.5	#$499,c.p
	move.a	c,a
	clr.w	b
	not.m	b
	rtn

___packsb
	;
	; Round extended real in A.W/B.W to real in A.W
	;
	; CPU must be in DEC mode
	;
	; If SB is set, rounding direction is determined by
	; the last nibble in B.M (Mantissa of A.W/B.W)
	;
	move.1	#0,p
	bne.m	b,0,.t0
	clr.w	a		; Zero mantissa, return 0
	rtn
.t0
	add.x	b,b
	bcc	___pack.t3.
	clr.m	c
	inc.m	c
	bne.x	b,0,___pack.t2.
	beq	#2,hst,___pack.t1.
	bcc	___pack.t2.

___pack
	;
	; Round extended real in A.W/B.W to real in A.W
	;
	; CPU must be in DEC mode
	;
	; possible errors:
	;
	;   Positive Underflow  -->  + 0
	;   Negative Underflow  -->  - 0
	;   Overflow            -->  +/- Inf
	;
	move.1	#0,p
	bne.m	b,0,.t0
	clr.w	a		; Zero mantissa, return 0
	rtn
.t0
	add.x	b,b
	bcc	.t3		; Round down (truncate)
	clr.m	c
	inc.m	c
	bcc	.t2
.t1
	and.m	b,c
.t2
	add.m	c,b		; Round up
	bcc	.t3
	move.1	#14,p		; Overflow on mantissa
	inc.p	b
	move.1	#0,p
	inc.a	a		; Increment exponent
.t3
	move.a	a,c
	add.a	c,c
	bcs	.t5
	move.5	#$499,c
	blt.a	c,a,.t7
.t4
	move.m	b,a
	rtn
.t5
	move.5	#$99501,c
	ble.a	c,a,.t4
	beq.s	a,0,.t6
	clr.w	a		; Negative underflow
	not.s	a
	rtn
.t6
	clr.w	a		; Positive underflow
	rtn
.t7
	move.3	#$499,c.p	; Overflow
	move.x	c,a
	clr.m	a
	not.m	a
	rtn


___trunctfdf2
	;
	; Round extended real in B.W/R0.W (???) to real in C.W
	;
	move.w	r0,c
	move.w	b,a
	move.w	c,b
	setdec
	jsr	___pack
	sethex
	move.1	#7,p
	move.w	a,c
	rtn

___extenddftf2
	;
	; Extend real in B.W to extended real in A.W/C.W (???)
	;
	move.w	b,a
	jsr	SPLITA		; calls setdec
	sethex
	move.1	#7,p
	move.w	b,c
	rtn
	

___adddf3
	;
	; Floating point C.W := B.W + R0.W
	;
	exg.a   c,d1
	push			; save frame-pointer
	exg.a   c,d1
	move.w	r0,c
	move.w	b,a
	jsr	SPLITAC		; calls setdec
	clr.1	#3,hst
	jsr	ADDF
	jsr	___packsb
	sethex			; restore hexmode
	move.1	#7,p		; restore p
	pop
	move.a	c,d1		; restore frame-pointer
	move.w	a,c
	rtn

___subdf3
	;
	; Floating point C.W := B.W - R0.W
	;
	exg.a   c,d1
	push			; save frame-pointer
	exg.a   c,d1
	move.w	r0,c
	move.w	b,a
	jsr	SPLITAC		; calls setdec
	not.s	c		; change sign
	clr.1	#3,hst
	jsr	ADDF
	jsr	___packsb
	sethex			; restore hexmode
	move.1	#7,p		; restore p
	pop
	move.a	c,d1		; restore frame-pointer
	move.w	a,c
	rtn

___muldf3
	;
	; C.W := B.W * R0.W
	;
	exg.a   c,d1
	push			; save frame-pointer
	exg.a   c,d1
	move.w	r0,c
	move.w	b,a
	jsr	SPLITAC		; calls setdec
	clr.1	#3,hst
	jsr	MULTF
	jsr	___packsb
	sethex			; restore hexmode
	move.1	#7,p		; restore p
	pop
	move.a	c,d1		; restore frame-pointer
	move.w	a,c
	rtn

___divdf3
	;
	; C.W := B.W / R0.W
	;
	exg.a	c,d1
	push			; save frame-pointer
	exg.a	c,d1
	move.w	r0,c
	move.w	b,a
	jsr	SPLITAC		; calls setdec
	clr.1	#3,hst
	jsr	DIVF
	jsr	___ckinf
	jsr	___packsb
	sethex			; restore hexmode
	move.1	#7,p		; restore p
	pop
	move.a	c,d1		; restore frame-pointer
	move.w	a,c
	rtn


___cmpdf
	;
	; Floating point compare A.W and C.W according to P
	;
	; P == 1:	CARRY := (A.W <  C.W)
	; P == 2:	CARRY := (A.W == C.W)
	; P == 3:	CARRY := (A.W <= C.W)
	; P == 4:	CARRY := (A.W >  C.W)
	; P == 6:	CARRY := (A.W >= C.W)
	; P == 13:	CARRY := (A.W != C.W)
	;
	exg.a   c,d1
	push			; save frame-pointer
	exg.a   c,d1
	jsr	SPLITAC		; calls setdec
	jsr	TST15
	sethex			; restore hexmode
	move.1 #7,p		; restore p
	pop
	move.a	c,d1		; restore frame-pointer
	rtn

___ltdf2
	move.1	#1,p		; compare is <
	bra.3	___cmpdf

___eqdf2
	move.1	#2,p		; compare is ==
	bra.3	___cmpdf

___ledf2
	move.1	#3,p		; compare is <=
	bra.3	___cmpdf

___gtdf2
	move.1	#4,p		; compare is >
	bra.3	___cmpdf

___gedf2
	move.1	#6,p		; compare is >=
	bra.3	___cmpdf

___nedf2
	move.1	#13,p		; compare is !=
	bra.3	___cmpdf


___fixdfsi
	;
	; Convert real in B.W to signed short in C.A
	;
	move.s	b,c
	move.s	c,r0		; remember sign
	clr.s	b
	bsr.3	___fixunsdfti
	clr.w	a
	move.1	#0,p
	move.4	#7fff,a.p
	move.1	#7,p
	ble.w	c,a,.t0		; c <= SHORT_MAX
	move.a	a,c		; c == SHORT_MAX
.t0
	move.s	r0,a
	rteq.s	a,0		; positive, return
	move.1	#0,p
	move.4	#$8000,a.p
	move.1	#7,p
	or.a	a,c		; c = -c
	rtn

___fixunsdfsi
	;
	; Convert real in B.W to unsigned short in C.A
	;
	bsr.3	___fixunsdfti
	clr.w	a
	move.1	#0,p
	move.4	#$ffff,a.p
	move.1	#7,p
	rtle.w	c,a		; return C.A
	move.a	a,c		; return USHORT_MAX
	rtn

___fixdfdi
	;
	; Convert real in B.W to signed long in C.WP
	;
	move.s	b,c
	move.s	c,r0		; remember sign
	clr.s	b
	bsr.3	___fixunsdfti
	move.w	c,a
	add.p	a,a		; n > LONG_MAX
	bcs	.t0
	clr.wp	a
	beq.w	a,0,.t1		; n <= LONG_MAX
.t0
	clr.wp	c
	not.wp	c
	lsr.wp	#1,c		; C.WP = LONG_MAX
.t1
	move.s	r0,a
	rteq.s	a,0		; positive, return c
	neg.wp	c		; c = -c
	rtn

___fixunsdfdi
	;
	; Convert real in B.W to unsigned long in C.WP
	;
	bsr.3	___fixunsdfti
	move.w	c,a		; test high 32 bits in C.W
	clr.wp	a
	rteq.w	a,0		; n <= ULONG_MAX, return C.WP
	clr.wp	c		; n > ULONG_MAX, return ULONG_MAX
	not.wp	c
	rtn

___fixdfti
	;
	; Convert real in B.W to signed long long in C.W
	;
	move.s	b,c
	move.s	c,r0		; remember sign
	clr.s	b
	bsr.3	___fixunsdfti
	move.w	c,a
	add.s	a,a
	bcc	.t0
	clr.w	c		; C > LONG_LONG_MAX
	not.w	c
	lsr.w	#1,c
.t0
	move.s	r0,a
	rteq.s	a,0		; positive, return
	neg.w	c		; c = -c
	rtn

___fixunsdfti
	;
	; Convert real in B.W to unsigned long long in C.W
	;
	sethex
	beq.s	b,0,.t3
.t1				; n <= 0
	move.1	#7,p
	clr.w	c
	rtn
.t2				; n >= ULONG_LONG_MAX
	move.1	#7,p
	clr.w	c
	dec.w	c
	rtn
.t3
	beq.w	b,0,.t1		; n == 0
	move.1	#0,p
	move.3	#$500,c.p
	bgt.x	c,b,.t4		; exp >= 0, calculate n
	move.3	#$999,c.p
	bne.x	b,c,.t1		; exp < -1, return 0
	move.1	#14,p
	move.1	#4,c.p
	sub.p	b,c
	bcc	.t1		; n < 0.5, return 0
	move.1	#7,p
	clr.w	c		; 0.5 <= n < 1, return 1
	inc.w	c
	rtn
.t4
	move.3	#$19,c.p
	blt.x	c,b,.t2		; exp > 19, return ULONG_LONG_MAX
	move.x	b,c
	move.x	b,a
	clr.x	b
	move.1	#14,p
	setdec
.t5
	dec.1	p		; locate the floating point
	beq.1	#2,p,.t6
	dec.x	c
	bcc	.t5
	add.p	b,b		; round B.W at floating point position
	clr.wp	b		; truncate
	bcc	.t6
	dec.wp	b		; round up
	inc.w	b
	beq.s	b,0,.t6
	lsr.w	#4,b		; overflow in mantissa, B.W /= 10,
	inc.x	a		; increment exponent
.t6
	clr.w	c
	sethex
.t7				; calculate C.W *= 10
	add.w	c,c		; C.W = 2 * c
	bcc	.t9
.t8
	bra.3	.t2		; overflow in C.W, return ULONG_LONG_MAX
.t9
	move.w	c,d		; D.W = 2 * c
	add.w	c,c		; C.W = 4 * c
	bcs	.t8
	add.w	c,c		; C.W = 8 * c
	bcs	.t8
	add.w	d,c		; C.W = 10 * c
	bcs	.t8
	lsl.w	#4,b
	beq.s	b,0,.t11
.t10				; C.W += high mantissa nibble
	inc.w	c
	dec.s	b
	bne.s	b,0,.t10
.t11
	setdec
	dec.x	a
	sethex
	bcc	.t7
	move.1	#7,p
	rtn