; Copyright (c) 1985 Regents of the University of California.
; All rights reserved.
;
; Redistribution and use in source and binary forms, with or without
; modification, are permitted provided that the following conditions
; are met:
; 1. Redistributions of source code must retain the above copyright
;    notice, this list of conditions and the following disclaimer.
; 2. Redistributions in binary form must reproduce the above copyright
;    notice, this list of conditions and the following disclaimer in the
;    documentation and/or other materials provided with the distribution.
; 3. All advertising materials mentioning features or use of this software
;    must display the following acknowledgement:
;	This product includes software developed by the University of
;	California, Berkeley and its contributors.
; 4. Neither the name of the University nor the names of its contributors
;    may be used to endorse or promote products derived from this software
;    without specific prior written permission.
;
; THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
; ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
; IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
; ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
; FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
; OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
; HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
; LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
; OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
; SUCH DAMAGE.
;
;	from: @(#)support.s	5.4 (Berkeley) 10/9/90
;	$Id: support.S,v 1.2 1993/08/01 18:34:22 mycroft Exp $
;

; IEEE recommended functions
; 
; double copysign(x,y)
; double x,y;
; IEEE 754 recommended function, return x*sign(y)
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
;
	.vers	2
	.text
	.align	2
	.globl	_copysign
_copysign:
	movl	4(sp),f0
	movd	8(sp),r0
	movd	16(sp),r1
	xord	r0,r1
	andd	0x80000000,r1
	cmpqd	0,r1
	beq	end
	negl	f0,f0
end:	ret	0

; 
; double logb(x)
; double x;
; IEEE p854 recommended function, return the exponent of x (return float(N) 
; such that 1 <= x*2**-N < 2, even for subnormal number.
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of. 
;
	.vers	2
	.text
	.align	2
	.globl	_logb
_logb:
;
; extract the exponent of x
; glossaries:	r0 = high part of x
;		r1 = unbias exponent of x
;		r2 = 20 (first exponent bit position)
;
	movd	8(sp),r0
	movd	20,r2
	extd	r2,r0,r1,11	; extract the exponent of x
	cmpqd	0,r1		; if exponent bits = 0, goto L3
	beq	L3
	cmpd	0x7ff,r1
	beq	L2		; if exponent bits = 0x7ff, goto L2
L1:	subd	1023,r1		; unbias the exponent
	movdl	r1,f0		; convert the exponent to floating value
	ret	0
;
; x is INF or NaN, simply return x
;
L2:
	movl	4(sp),f0	; logb(+inf)=+inf, logb(NaN)=NaN
	ret	0
;
; x is 0 or subnormal
;
L3:
	movl	4(sp),f0
	cmpl	0f0,f0
	beq	L5		; x is 0 , goto L5 (return -inf)
;
; Now x is subnormal
;
	mull	L64,f0		; scale up f0 with 2**64
	movl	f0,tos
	movd	tos,r0
	movd	tos,r0		; now r0 = new high part of x
	extd	r2,r0,r1,11	; extract the exponent of x to r1
	subd	1087,r1		; unbias the exponent with correction 
	movdl	r1,f0		; convert the exponent to floating value
	ret	0
;
; x is 0, return logb(0)= -INF
;
L5:
	movl	0f1.0e300,f0
	mull	0f-1.0e300,f0	; multiply two big numbers to get -INF
	ret	0
; 
; double rint(x)
; double x;
; ... delivers integer nearest x in direction of prevailing rounding
; ... mode
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of. 
;
	.vers	2
	.text
	.align	2
	.globl	_rint
_rint:
;
	movd	8(sp),r0
	movd	20,r2
	extd	r2,r0,r1,11	; extract the exponent of x
	cmpd	0x433,r1
	ble	itself
	movl	L52,f2		; f2 = L = 2**52
	cmpqd	0,r0
	ble	L1
	negl	f2,f2		; f2 = s = copysign(L,x)
L1:	addl	f2,f0		; f0 = x + s
	subl	f2,f0		; f0 = f0 - s
	ret	0
itself:	movl	4(sp),f0
	ret	0
L52:	.double	0x0,0x43300000	; L52=2**52
; 
; int finite(x)
; double x;
; IEEE 754 recommended function, return 0 if x is NaN or INF, else 0
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
;
	.vers	2
	.text
	.align	2
	.globl	_finite
_finite:
	movd	4(sp),r1
	andd	0x800fffff,r1
	cmpd	0x7ff00000,r1
	sned	r0		; r0=0 if exponent(x) = 0x7ff
	ret	0
; 
; double scalb(x,N)
; double x; int N;
; IEEE 754 recommended function, return x*2**N by adjusting 
; exponent of x.
; Coded by K.C. Ng in National 32k assembler, 11/9/85. 
; Note: subnormal number (if implemented) will be taken care of 
;
	.vers	2
	.text
	.align	2
	.globl	_scalb
_scalb:
;
; if x=0 return 0
;
	movl	4(sp),f0
	cmpl	0f0,f0
	beq	end		; scalb(0,N) is x itself
;
; extract the exponent of x
; glossaries:	r0 = high part of x, 
;		r1 = unbias exponent of x,
;		r2 = 20 (first exponent bit position).
;
	movd	8(sp),r0	; r0 = high part of x
	movd	20,r2		; r2 = 20
	extd	r2,r0,r1,11	; extract the exponent of x in r1
	cmpd	0x7ff,r1	
;
; if exponent of x is 0x7ff, then x is NaN or INF; simply return x  
;
	beq	end		
	cmpqd	0,r1
;
; if exponent of x is zero, then x is subnormal; goto L19
;
	beq	L19		
	addd	12(sp),r1	; r1 = (exponent of x) + N
	bfs	inof		; if integer overflows, goto inof
	cmpqd	0,r1		; if new exponent <= 0, goto underflow
	bge	underflow
	cmpd	2047,r1		; if new exponent >= 2047 goto overflow
	ble	overflow
	insd	r2,r1,r0,11	; insert the new exponent 
	movd	r0,tos
	movd	8(sp),tos
	movl	tos,f0		; return x*2**N
end:	ret	0
inof:	bcs	underflow	; negative int overflow if Carry bit is set
overflow:
	andd	0x80000000,r0	; keep the sign of x
	ord	0x7fe00000,r0	; set x to a huge number
	movd	r0,tos
	movqd	0,tos
	movl	tos,f0
	mull	0f1.0e300,f0	; multiply two huge number to get overflow
	ret	0
underflow:
	addd	64,r1		; add 64 to exonent to see if it is subnormal
	cmpqd	0,r1
	bge	zero		; underflow to zero
	insd	r2,r1,r0,11	; insert the new exponent 
	movd	r0,tos
	movd	8(sp),tos
	movl	tos,f0
	mull	L30,f0		; readjust x by multiply it with 2**-64
	ret	0
zero:	andd	0x80000000,r0	; keep the sign of x
	ord	0x00100000,r0	; set x to a tiny number
	movd	r0,tos
	movqd	0,tos
	movl	tos,f0
	mull	0f1.0e-300,f0	; underflow to 0  by multipling two tiny nos.
	ret	0
L19:		; subnormal number
	mull	L32,f0		; scale up x by 2**64
	movl	f0,tos
	movd	tos,r0
	movd	tos,r0		; get the high part of new x
	extd	r2,r0,r1,11	; extract the exponent of x in r1
	addd	12(sp),r1	; exponent of x + N
	subd	64,r1		; adjust it by subtracting 64
	cmpqd	0,r1
	bge	underflow
	cmpd	2047,r1
	ble	overflow
	insd	r2,r1,r0,11	; insert back the incremented exponent 
	movd	r0,tos
	movd	8(sp),tos
	movl	tos,f0
end:	ret	0
L30:	.double	0x0,0x3bf00000	; floating point 2**-64
L32:	.double	0x0,0x43f00000	; floating point 2**64
; 
; double drem(x,y)
; double x,y;
; IEEE double remainder function, return x-n*y, where n=x/y rounded to 
; nearest integer (half way case goes to even). Result exact.
; Coded by K.C. Ng in National 32k assembly, 11/19/85.
;
	.vers	2
	.text
	.align	2
	.globl	_drem
_drem:
;
; glossaries:	
;		r2 = high part of x
;		r3 = exponent of x
;		r4 = high part of y
;		r5 = exponent of y
;		r6 = sign of x
;		r7 = constant 0x7ff00000
;
;  16(fp) : y
;   8(fp) : x
; -12(fp) : adjustment on y when y is subnormal
; -16(fp) : fsr
; -20(fp) : nx
; -28(fp) : t
; -36(fp) : t1
; -40(fp) : nf
; 
;
	enter	[r3,r4,r5,r6,r7],40
	movl	f6,tos
	movl	f4,tos
	movl	0f0,-12(fp)
	movd	0,-20(fp)
	movd	0,-40(fp)
	movd	0x7ff00000,r7	; initialize r7=0x7ff00000
	movd	12(fp),r2	; r2 = high(x)
	movd	r2,r3
	andd	r7,r3		; r3 = xexp
	cmpd	r7,r3
; if x is NaN or INF goto L1
	beq	L1		
	movd	20(fp),r4
	bicd	[31],r4		; r4 = high part of |y|
	movd	r4,20(fp)	; y = |y|
	movd	r4,r5
	andd	r7,r5		; r5 = yexp
	cmpd	r7,r5
	beq	L2		; if y is NaN or INF goto L2
	cmpd	0x04000000,r5	; 
	bgt	L3		; if y is tiny goto L3
;
; now y != 0 , x is finite
;
L10:
	movd	r2,r6
	andd	0x80000000,r6	; r6 = sign(x)
	bicd	[31],r2		; x <- |x|
	sfsr	r1
	movd	r1,-16(fp)	; save fsr in -16(fp)
	bicd	[5],r1
	lfsr	r1		; disable inexact interupt
	movd	16(fp),r0	; r0 = low part of y
	movd	r0,r1		; r1 = r0 = low part of y
	andd	0xf8000000,r1	; mask off the lsb 27 bits of y

	movd	r2,12(fp)	; update x to |x|
	movd	r0,-28(fp)	; 
	movd	r4,-24(fp)	; t  = y
	movd	r4,-32(fp)	; 
	movd	r1,-36(fp)	; t1 = y with trialing 27 zeros
	movd	0x01900000,r1	; r1 = 25 in exponent field
LOOP:
	movl	8(fp),f0	; f0 = x
	movl	16(fp),f2	; f2 = y
	cmpl	f0,f2
	ble	fnad		; goto fnad (final adjustment) if x <= y
	movd	r4,-32(fp)
	movd	r3,r0
	subd	r5,r0		; xexp - yexp
	subd	r1,r0		; r0 = xexp - yexp - m25
	cmpqd	0,r0		; r0 > 0 ?
	bge	1f
	addd	r4,r0		; scale up (high) y
	movd	r0,-24(fp)	; scale up t
	movl	-28(fp),f2	; t
	movd	r0,-32(fp)	; scale up t1
1:
	movl	-36(fp),f4	; t1
	movl	f0,f6
	divl	f2,f6		; f6 = x/t
	floorld	f6,r0		; r0 = [x/t]
	movdl	r0,f6		; f6 = n
	subl	f4,f2		; t = t - t1 (tail of t1)
	mull	f6,f4		; f4 = n*t1	...exact
	subl	f4,f0		; x = x - n*t1
	mull	f6,f2		; n*(t-t1)	...exact
	subl	f2,f0		; x = x - n*(t-t1)
; update xexp
	movl	f0,8(fp)
	movd	12(fp),r3	
	andd	r7,r3
	jump	LOOP
fnad:
	cmpqd	0,-20(fp)	; 0 = nx?
	beq	final
	mull	-12(fp),8(fp)	; scale up x the same amount as y
	movd	0,-20(fp)
	movd	12(fp),r2
	movd	r2,r3
	andd	r7,r3		; update exponent of x
	jump	LOOP

final:
	movl	16(fp),f2	; f2 = y (f0=x, r0=n)
	subd	0x100000,r4	; high y /2
	movd	r4,-24(fp)
	movl	-28(fp),f4	; f4 = y/2
	cmpl	f0,f4		; x > y/2 ?
	bgt	1f
	bne	2f
	andd	1,r0		; n is odd or even
	cmpqd	0,r0
	beq	2f
1:
	subl	f2,f0		; x = x - y
2:
	cmpqd	0,-40(fp)
	beq	3f
	divl	-12(fp),f0	; scale down the answer
3:
	movl	f0,tos
	xord	r6,tos
	movl	tos,f0
	movd	-16(fp),r0
	lfsr	r0		; restore the fsr

end:	movl	tos,f4
	movl	tos,f6
	exit	[r3,r4,r5,r6,r7]
	ret	0
;
; y is NaN or INF
;
L2:	
	movd	16(fp),r0	; r0 = low part of y
	andd	0xfffff,r4	; r4 = high part of y & 0x000fffff
	ord	r4,r0
	cmpqd	0,r0
	beq	L4
	movl	16(fp),f0	; y is NaN, return y
	jump	end
L4:	movl	8(fp),f0	; y is inf, return x
	jump	end
;
; exponent of y is less than 64, y may be zero or subnormal
;
L3:
	movl	16(fp),f0
	cmpl	0f0,f0
	bne	L5
	divl	f0,f0		; y is 0, return NaN by doing 0/0
	jump	end
;
; subnormal y or tiny y
;
L5:	
	movd	0x04000000,-20(fp)	; nx = 64 in exponent field
	movl	L64,f2
	movl	f2,-12(fp)
	mull	f2,f0
	cmpl	f0,LTINY
	bgt	L6
	mull	f2,f0
	addd	0x04000000,-20(fp)	; nx = nx + 64 in exponent field
	mull	f2,-12(fp)
L6:
	movd	-20(fp),-40(fp)
	movl	f0,16(fp)
	movd	20(fp),r4
	movd	r4,r5
	andd	r7,r5		; exponent of new y
	jump	L10
;
; x is NaN or INF, return x-x
;
L1:
	movl	8(fp),f0
	subl	f0,f0		; if x is INF, then INF-INF is NaN
	ret	0
L64:	.double 0x0,0x43f00000	; L64 = 2**64
LTINY:	.double 0x0,0x04000000	; LTINY = 2**-959