[unix-history] / usr / src / lib / libm / vax / support.s

/*
 * Copyright (c) 1985 Regents of the University of California.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms are permitted
 * provided that this notice is preserved and that due credit is given
 * to the University of California at Berkeley. The name of the University
 * may not be used to endorse or promote products derived from this
 * software without specific prior written permission. This software
 * is provided ``as is'' without express or implied warranty.
 *
 * All recipients should regard themselves as participants in an ongoing
 * research project and hence should feel obligated to report their
 * experiences (good or bad) with these elementary function codes, using
 * the sendbug(8) program, to the authors.
 *
 *	@(#)support.s	5.2 (Berkeley) %G%
 */
	.data
	.align	2
_sccsid:
.asciz	"@(#)support.s	1.3 (Berkeley) 8/21/85; 5.2 (ucb.elefunt) %G%"

/*
 * copysign(x,y),
 * logb(x),
 * scalb(x,N),
 * finite(x),
 * drem(x,y),
 * Coded in vax assembly language by K.C. Ng,  3/14/85.
 * Revised by K.C. Ng on 4/9/85.
 */

/*
 * double copysign(x,y)
 * double x,y;
 */
	.globl	_copysign
	.text
	.align	1
_copysign:
	.word	0x4
	movq	4(ap),r0		# load x into r0
	bicw3	$0x807f,r0,r2		# mask off the exponent of x
	beql	Lz			# if zero or reserved op then return x
	bicw3	$0x7fff,12(ap),r2	# copy the sign bit of y into r2
	bicw2	$0x8000,r0		# replace x by |x|
	bisw2	r2,r0			# copy the sign bit of y to x
Lz:	ret

/*
 * double logb(x)
 * double x;
 */
	.globl	_logb
	.text
	.align	1
_logb:
	.word	0x0
	bicl3	$0xffff807f,4(ap),r0	# mask off the exponent of x
	beql    Ln
	ashl	$-7,r0,r0		# get the bias exponent
	subl2	$129,r0			# get the unbias exponent
	cvtld	r0,r0			# return the answer in double
	ret
Ln:	movq	4(ap),r0		# r0:1 = x (zero or reserved op)
	bneq	1f			# simply return if reserved op
	movq 	$0x0000fe00ffffcfff,r0  # -2147483647.0
1:	ret

/*
 * long finite(x)
 * double x;
 */
	.globl	_finite
	.text
	.align	1
_finite:
	.word	0x0000
	bicw3	$0x7f,4(ap),r0		# mask off the mantissa
	cmpw	r0,$0x8000		# to see if x is the reserved op
	beql	1f			# if so, return FALSE (0)
	movl	$1,r0			# else return TRUE (1)
	ret
1:	clrl	r0
	ret

/*
 * double scalb(x,N)
 * double x; int N;
 */
	.globl	_scalb
	.set	ERANGE,34
	.text
	.align	1
_scalb:
	.word	0xc
	movq	4(ap),r0
	bicl3	$0xffff807f,r0,r3
	beql	ret1			# 0 or reserved operand
	movl	12(ap),r2
	cmpl	r2,$0x12c
	bgeq	ovfl
	cmpl	r2,$-0x12c
	bleq	unfl
	ashl	$7,r2,r2
	addl2	r2,r3
	bleq	unfl
	cmpl	r3,$0x8000
	bgeq	ovfl
	addl2	r2,r0
	ret
ovfl:	pushl	$ERANGE
	calls	$1,_infnan		# if it returns
	bicw3	$0x7fff,4(ap),r2	# get the sign of input arg
	bisw2	r2,r0			# re-attach the sign to r0/1
	ret
unfl:	movq	$0,r0
ret1:	ret

/*
 * DREM(X,Y)
 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
 * DOUBLE PRECISION (VAX D format 56 bits)
 * CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/8/85.
 */
	.globl	_drem
	.set	EDOM,33
	.text
	.align	1
_drem:
	.word	0xffc
	subl2	$12,sp	
	movq	4(ap),r0		#r0=x
	movq	12(ap),r2		#r2=y
	jeql	Rop			#if y=0 then generate reserved op fault
	bicw3	$0x007f,r0,r4		#check if x is Rop
	cmpw	r4,$0x8000
	jeql	Ret			#if x is Rop then return Rop
	bicl3	$0x007f,r2,r4		#check if y is Rop
	cmpw	r4,$0x8000
	jeql	Ret			#if y is Rop then return Rop
	bicw2	$0x8000,r2		#y  := |y|
	movw	$0,-4(fp)		#-4(fp) = nx := 0
	cmpw	r2,$0x1c80		#yexp ? 57 
	bgtr	C1			#if yexp > 57 goto C1
	addw2	$0x1c80,r2		#scale up y by 2**57
	movw	$0x1c80,-4(fp)		#nx := 57 (exponent field)
C1:
	movw	-4(fp),-8(fp)		#-8(fp) = nf := nx
	bicw3	$0x7fff,r0,-12(fp)	#-12(fp) = sign of x
	bicw2	$0x8000,r0		#x  := |x|
	movq	r2,r10			#y1 := y
	bicl2	$0xffff07ff,r11		#clear the last 27 bits of y1
loop:
	cmpd	r0,r2			#x ? y
	bleq	E1			#if x <= y goto E1
 /* begin argument reduction */
	movq	r2,r4			#t =y
	movq	r10,r6			#t1=y1
	bicw3	$0x807f,r0,r8		#xexp= exponent of x
	bicw3	$0x807f,r2,r9		#yexp= exponent fo y
	subw2	r9,r8			#xexp-yexp
	subw2	$0x0c80,r8		#k=xexp-yexp-25(exponent bit field)
	blss	C2			#if k<0 goto C2
	addw2	r8,r4			#t +=k	
	addw2	r8,r6			#t1+=k, scale up t and t1
C2:
	divd3	r4,r0,r8		#x/t
	cvtdl	r8,r8			#n=[x/t] truncated
	cvtld	r8,r8			#float(n)
	subd2	r6,r4			#t:=t-t1
	muld2	r8,r4			#n*(t-t1)
	muld2	r8,r6			#n*t1
	subd2	r6,r0			#x-n*t1
	subd2	r4,r0			#(x-n*t1)-n*(t-t1)
	brb	loop
E1:
	movw	-4(fp),r6		#r6=nx
	beql	C3			#if nx=0 goto C3
	addw2	r6,r0			#x:=x*2**57 scale up x by nx
	movw	$0,-4(fp)		#clear nx
	brb	loop
C3:
	movq	r2,r4			#r4 = y
	subw2	$0x80,r4		#r4 = y/2
	cmpd	r0,r4			#x:y/2
	blss	E2			#if x < y/2 goto E2
	bgtr	C4			#if x > y/2 goto C4
	cvtdl	r8,r8			#ifix(float(n))
	blbc	r8,E2			#if the last bit is zero, goto E2
C4:
	subd2	r2,r0			#x-y
E2:
	xorw2	-12(fp),r0		#x^sign (exclusive or)
	movw	-8(fp),r6		#r6=nf
	bicw3	$0x807f,r0,r8		#r8=exponent of x
	bicw2	$0x7f80,r0		#clear the exponent of x
	subw2	r6,r8			#r8=xexp-nf
	bgtr	C5			#if xexp-nf is positive goto C5
	movw	$0,r8			#clear r8
	movq	$0,r0			#x underflow to zero
C5:
	bisw2	r8,r0			#put r8 into x's exponent field
	ret
Rop:					#Reserved operand
	pushl	$EDOM
	calls	$1,_infnan		#generate reserved op fault
	ret
Ret:
	movq	$0x8000,r0		#propagate reserved op
	ret
Commit	Line	Data
fe5be67b	1	/*
04bb685f	2	* Copyright (c) 1985 Regents of the University of California.
fe5be67b KB	3	* All rights reserved.
	4	*
	5	* Redistribution and use in source and binary forms are permitted
	6	* provided that this notice is preserved and that due credit is given
	7	* to the University of California at Berkeley. The name of the University
	8	* may not be used to endorse or promote products derived from this
	9	* software without specific prior written permission. This software
	10	* is provided ``as is'' without express or implied warranty.
	11	*
	12	* All recipients should regard themselves as participants in an ongoing
	13	* research project and hence should feel obligated to report their
	14	* experiences (good or bad) with these elementary function codes, using
	15	* the sendbug(8) program, to the authors.
	16	*
	17	* @(#)support.s 5.2 (Berkeley) %G%
7c0a3811 GK	18	*/
	19	.data
	20	.align 2
	21	_sccsid:
fe5be67b	22	.asciz "@(#)support.s 1.3 (Berkeley) 8/21/85; 5.2 (ucb.elefunt) %G%"
7c0a3811 GK	23
7c0a3811 GK	24	/*
04bb685f ZAL	25	* copysign(x,y),
	26	* logb(x),
	27	* scalb(x,N),
	28	* finite(x),
	29	* drem(x,y),
	30	* Coded in vax assembly language by K.C. Ng, 3/14/85.
	31	* Revised by K.C. Ng on 4/9/85.
	32	*/
	33
	34	/*
	35	* double copysign(x,y)
	36	* double x,y;
	37	*/
	38	.globl _copysign
	39	.text
	40	.align 1
	41	_copysign:
	42	.word 0x4
	43	movq 4(ap),r0 # load x into r0
	44	bicw3 $0x807f,r0,r2 # mask off the exponent of x
	45	beql Lz # if zero or reserved op then return x
	46	bicw3 $0x7fff,12(ap),r2 # copy the sign bit of y into r2
	47	bicw2 $0x8000,r0 # replace x by \|x\|
	48	bisw2 r2,r0 # copy the sign bit of y to x
	49	Lz: ret
	50
	51	/*
	52	* double logb(x)
	53	* double x;
	54	*/
	55	.globl _logb
	56	.text
	57	.align 1
	58	_logb:
	59	.word 0x0
	60	bicl3 $0xffff807f,4(ap),r0 # mask off the exponent of x
	61	beql Ln
	62	ashl $-7,r0,r0 # get the bias exponent
	63	subl2 $129,r0 # get the unbias exponent
	64	cvtld r0,r0 # return the answer in double
	65	ret
	66	Ln: movq 4(ap),r0 # r0:1 = x (zero or reserved op)
	67	bneq 1f # simply return if reserved op
	68	movq $0x0000fe00ffffcfff,r0 # -2147483647.0
	69	1: ret
	70
	71	/*
	72	* long finite(x)
	73	* double x;
	74	*/
	75	.globl _finite
	76	.text
	77	.align 1
	78	_finite:
	79	.word 0x0000
	80	bicw3 $0x7f,4(ap),r0 # mask off the mantissa
	81	cmpw r0,$0x8000 # to see if x is the reserved op
	82	beql 1f # if so, return FALSE (0)
	83	movl $1,r0 # else return TRUE (1)
	84	ret
	85	1: clrl r0
	86	ret
	87
	88	/*
89	* double scalb(x,N)
90	* double x; int N;
91	*/
92	.globl _scalb
93	.set ERANGE,34
94	.text
95	.align 1
96	_scalb:
97	.word 0xc
98	movq 4(ap),r0
99	bicl3 $0xffff807f,r0,r3
100	beql ret1 # 0 or reserved operand
101	movl 12(ap),r2
102	cmpl r2,$0x12c
103	bgeq ovfl
104	cmpl r2,$-0x12c
105	bleq unfl
106	ashl $7,r2,r2
107	addl2 r2,r3
108	bleq unfl
109	cmpl r3,$0x8000
110	bgeq ovfl
111	addl2 r2,r0
112	ret
113	ovfl: pushl $ERANGE
114	calls $1,_infnan # if it returns
115	bicw3 $0x7fff,4(ap),r2 # get the sign of input arg
116	bisw2 r2,r0 # re-attach the sign to r0/1
117	ret
118	unfl: movq $0,r0
119	ret1: ret
120
121	/*
122	* DREM(X,Y)
123	* RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
124	* DOUBLE PRECISION (VAX D format 56 bits)
125	* CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/8/85.
126	*/
127	.globl _drem
128	.set EDOM,33
129	.text
130	.align 1
131	_drem:
132	.word 0xffc
133	subl2 $12,sp
134	movq 4(ap),r0 #r0=x
135	movq 12(ap),r2 #r2=y
136	jeql Rop #if y=0 then generate reserved op fault
137	bicw3 $0x007f,r0,r4 #check if x is Rop
138	cmpw r4,$0x8000
139	jeql Ret #if x is Rop then return Rop
140	bicl3 $0x007f,r2,r4 #check if y is Rop
141	cmpw r4,$0x8000
142	jeql Ret #if y is Rop then return Rop
143	bicw2 $0x8000,r2 #y := \|y\|
144	movw $0,-4(fp) #-4(fp) = nx := 0
145	cmpw r2,$0x1c80 #yexp ? 57
146	bgtr C1 #if yexp > 57 goto C1
147	addw2 $0x1c80,r2 #scale up y by 2**57
148	movw $0x1c80,-4(fp) #nx := 57 (exponent field)
149	C1:
150	movw -4(fp),-8(fp) #-8(fp) = nf := nx
151	bicw3 $0x7fff,r0,-12(fp) #-12(fp) = sign of x
152	bicw2 $0x8000,r0 #x := \|x\|
153	movq r2,r10 #y1 := y
154	bicl2 $0xffff07ff,r11 #clear the last 27 bits of y1
155	loop:
156	cmpd r0,r2 #x ? y
157	bleq E1 #if x <= y goto E1
158	/* begin argument reduction */
159	movq r2,r4 #t =y
160	movq r10,r6 #t1=y1
161	bicw3 $0x807f,r0,r8 #xexp= exponent of x
162	bicw3 $0x807f,r2,r9 #yexp= exponent fo y
163	subw2 r9,r8 #xexp-yexp
164	subw2 $0x0c80,r8 #k=xexp-yexp-25(exponent bit field)
165	blss C2 #if k<0 goto C2
166	addw2 r8,r4 #t +=k
167	addw2 r8,r6 #t1+=k, scale up t and t1
168	C2:
169	divd3 r4,r0,r8 #x/t
170	cvtdl r8,r8 #n=[x/t] truncated
171	cvtld r8,r8 #float(n)
172	subd2 r6,r4 #t:=t-t1
173	muld2 r8,r4 #n*(t-t1)
174	muld2 r8,r6 #n*t1
175	subd2 r6,r0 #x-n*t1
176	subd2 r4,r0 #(x-nt1)-n(t-t1)
177	brb loop
178	E1:
179	movw -4(fp),r6 #r6=nx
180	beql C3 #if nx=0 goto C3
181	addw2 r6,r0 #x:=x2*57 scale up x by nx
182	movw $0,-4(fp) #clear nx
183	brb loop
184	C3:
185	movq r2,r4 #r4 = y
186	subw2 $0x80,r4 #r4 = y/2
187	cmpd r0,r4 #x:y/2
188	blss E2 #if x < y/2 goto E2
189	bgtr C4 #if x > y/2 goto C4
190	cvtdl r8,r8 #ifix(float(n))
191	blbc r8,E2 #if the last bit is zero, goto E2
192	C4:
193	subd2 r2,r0 #x-y
194	E2:
195	xorw2 -12(fp),r0 #x^sign (exclusive or)
196	movw -8(fp),r6 #r6=nf
197	bicw3 $0x807f,r0,r8 #r8=exponent of x
198	bicw2 $0x7f80,r0 #clear the exponent of x
199	subw2 r6,r8 #r8=xexp-nf
200	bgtr C5 #if xexp-nf is positive goto C5
201	movw $0,r8 #clear r8
202	movq $0,r0 #x underflow to zero
203	C5:
204	bisw2 r8,r0 #put r8 into x's exponent field
205	ret
206	Rop: #Reserved operand
207	pushl $EDOM
208	calls $1,_infnan #generate reserved op fault
209	ret
210	Ret:
211	movq $0x8000,r0 #propagate reserved op
212	ret