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

/* 
 * Copyright (c) 1987 Regents of the University of California.
 * 
 * Use and reproduction of this software are granted  in  accordance  with
 * the terms and conditions specified in  the  Berkeley  Software  License
 * Agreement (in particular, this entails acknowledgement of the programs'
 * source, and inclusion of this notice) with the additional understanding
 * that  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 "sendbug 4bsd-bugs@BERKELEY", to the authors.
 */
	.data
	.align	2
_sccsid:
	.asciz	"@(#)support.s	1.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 4/9/85.
 * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87.
 */
/*
 * double copysign(x,y)
 * double x,y;
 */
	.globl	_copysign
	.text
	.align	2
_copysign:
	.word	0x0004			# save r2
	movl	8(fp),r1
	movl	4(fp),r0		# r0:r1 = x
	andl3	$0x7f800000,r0,r2	# r2 = biased exponent of x
	beql	1f			# if 0 or reserved op then return x
	andl3	$0x80000000,12(fp),r2	# r2 = sign bit of y at bit-31
	andl2	$0x7fffffff,r0		# replace x by |x|
	orl2	r2,r0			# copy the sign bit of y to x
1:	ret
/*
 * double logb(x)
 * double x;
 */
	.globl	_logb
	.text
	.align	2
_logb:
	.word	0x0000			# save nothing
	andl3	$0x7f800000,4(fp),r0	# r0[b23:b30] = biased exponent of x
	beql    1f
	shrl	$23,r0,r0		# r0[b0:b7] = biased exponent of x
	subl2	$129,r0			# r0 = unbiased exponent of x
	cvld	r0			# acc = unbiased exponent of x (double)
	std	r0			# r0 =  unbiased exponent of x (double)
	ret
1:	movl	8(fp),r1		# 8(fp) must be moved first
	movl	4(fp),r0		# r0:r1 = x (zero or reserved op)
	blss	2f			# simply return if reserved op
	movl	$0xfe000000,r1
	movl	$0xcfffffff,r0		# -2147483647.0
2:	ret
/*
 * long finite(x)
 * double x;
 */
	.globl	_finite
	.text
	.align	2
_finite:
	.word	0x0000			# save nothing
	andl3	$0xff800000,4(fp),r0	# r0 = sign of x & its biased exponent
	cmpl	r0,$0x80000000		# is x a 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	2
_scalb:
	.word	0x000c			# save r2-r3
	movl	8(fp),r1
	movl	4(fp),r0		# r0:r1 = x (-128 <= Ex <= 126)
	andl3	$0x7f800000,r0,r3	# r3[b23:b30] = biased exponent of x
	beql	1f			# is x a 0 or a reserved operand?
	movl	12(fp),r2		# r2 = N
	cmpl	r2,$0xff		# if N >= 255
	bgeq	2f			# then the result must overflow
	cmpl	r2,$-0xff		# if N <= -255
	bleq	3f			# then the result must underflow
	shrl	$23,r3,r3		# r3[b0:b7] = biased exponent of x
	addl2	r2,r3			# r3 = biased exponent of the result
	bleq	3f			# if <= 0 then the result underflows
	cmpl	r3,$0x100		# if >= 256 then the result overflows
	bgeq	2f
	shll	$23,r3,r3		# r3[b23:b30] = biased exponent of res.
	andl2	$0x807fffff,r0
	orl2	r3,r0			# r0:r1 = x*2^N
1:	ret
2:	pushl	$ERANGE			# if the result would overflow
	callf	$8,_infnan		# and _infnan returns
	andl3	$0x80000000,4(fp),r2	# get the sign of input arg
	orl2	r2,r0			# re-attach the sign to r0:r1
	ret
3:	clrl	r1			# if the result would underflow
	clrl	r0			# then return 0
	ret
/*
 * double drem(x,y)
 * double x,y;
 * Returns x-n*y where n=[x/y] rounded (to even in the half way case).
 */
	.globl	_drem
	.set	EDOM,33
	.text
	.align	2
_drem:
	.word	0x1ffc			# save r2-r12
	movl	16(fp),r3
	movl	12(fp),r2		# r2:r3 = y
	movl	8(fp),r1
	movl	4(fp),r0		# r0:r1 = x
	andl3	$0xff800000,r0,r4
	cmpl	r4,$0x80000000		# is x a reserved operand?
	beql	1f			# if yes then propagate x and return
	andl3	$0xff800000,r2,r4
	cmpl	r4,$0x80000000		# is y a reserved operand?
	bneq	2f
	movl	r3,r1
	movl	r2,r0			# if yes then propagate y and return
1:	ret

2:	tstl	r4			# is y a 0?
	bneq	3f
	pushl	$EDOM			# if so then generate reserved op fault
	callf	$8,_infnan
	ret

3:	andl2	$0x7fffffff,r2		# r2:r3 = y <- |y|
	clrl	r12			# r12 = nx := 0
	cmpl	r2,$0x1c800000		# Ey ? 57 
	bgtr	4f			# if Ey > 57 goto 4
	addl2	$0x1c800000,r2		# scale up y by 2**57
	movl	$0x1c800000,r12		# r12[b23:b30] = nx = 57
4:	pushl	r12			# pushed onto stack: nf := nx
	andl3	$0x80000000,r0,-(sp)	# pushed onto stack: sign of x
	andl2	$0x7fffffff,r0		# r0:r1 = x <- |x|
	movl	r3,r11			# r10:r11 = y1 = y w/ last 27 bits 0
	andl3	$0xf8000000,r10,r11	# clear last 27 bits of y1

Loop:	cmpd2	r0,r2			# x ? y
	bleq	6f			# if x <= y goto 6
 /* 					# begin argument reduction */
	movl	r3,r5
	movl	r2,r4			# r4:r5 = t = y
	movl	r11,r7
	movl	r10,r6			# r6:r7 = t1 = y1
	andl3	$0x7f800000,r0,r8	# r8[b23:b30] = Ex:biased exponent of x
	andl3	$0x7f800000,r2,r9	# r9[b23:b30] = Ey:biased exponent of y
	subl2	r9,r8			# r8[b23:b30] = Ex-Ey
	subl2	$0x0c800000,r8		# r8[b23:b30] = k = Ex-Ey-25
	blss	5f			# if k < 0 goto 5
	addl2	r8,r4			# t += k	
	addl2	r8,r6			# t1 += k, scale up t and t1
5:	ldd	r0			# acc = x
	divd	r4			# acc = x/t
	cvdl	r8			# r8 = n = [x/t] truncated
	cvld	r8			# acc = dble(n)
	std	r8			# r8:r9 = dble(n)
	ldd	r4			# acc = t
	subd	r6			# acc = t-t1
	muld	r8			# acc = n*(t-t1)
	std	r4			# r4:r5 = n*(t-t1)
	ldd	r6			# acc = t1
	muld	r8			# acc = n*t1
	subd	r0			# acc = n*t1-x
	negd				# acc = x-n*t1
	subd	r4			# acc = (x-n*t1)-n*(t-t1)
	std	r0			# r0:r1 = (x-n*t1)-n*(t-t1)
	brb	Loop

6:	movl	r12,r6			# r6 = nx
	beql	7f			# if nx == 0 goto 7
	addl2	r6,r0			# x <- x*2**57:scale x up by nx
	clrl	r12			# clear nx
	brb	Loop

7:	movl	r3,r5
	movl	r2,r4			# r4:r5 = y
	subl2	$0x800000,r4		# r4:r5 = y/2
	cmpd2	r0,r4			# x ? y/2
	blss	9f			# if x < y/2 goto 9
	bgtr	8f			# if x > y/2 goto 8
	ldd	r8			# acc = dble(n)
	cvdl	r8			# r8 = ifix(dble(n))
	bbc	$0,r8,9f		# if the last bit is zero, goto 9
8:	ldd	r0			# acc = x
	subd	r2			# acc = x-y
	std	r0			# r0:r1 = x-y
9:	xorl2	(sp)+,r0		# x^sign (exclusive or)
	movl	(sp)+,r6		# r6 = nf
	andl3	$0x7f800000,r0,r8	# r8 = biased exponent of x
	andl2	$0x807fffff,r0		# r0 = x w/ exponent zapped
	subl2	r6,r8			# r8 = Ex-nf
	bgtr	0f			# if Ex-nf > 0 goto 0
	clrl	r8			# clear r8
	clrl	r0
	clrl	r1			# x underflows to zero
0:	orl2	r8,r0			# put r8 into x's exponent field
	ret
Commit	Line	Data
97158cb3 ZAL	1	/*
	2	* Copyright (c) 1987 Regents of the University of California.
	3	*
	4	* Use and reproduction of this software are granted in accordance with
	5	* the terms and conditions specified in the Berkeley Software License
	6	* Agreement (in particular, this entails acknowledgement of the programs'
	7	* source, and inclusion of this notice) with the additional understanding
	8	* that all recipients should regard themselves as participants in an
	9	* ongoing research project and hence should feel obligated to report
	10	* their experiences (good or bad) with these elementary function codes,
	11	* using "sendbug 4bsd-bugs@BERKELEY", to the authors.
	12	*/
	13	.data
	14	.align 2
	15	_sccsid:
f47ef219	16	.asciz "@(#)support.s 1.2 (ucb.elefunt) %G%"
97158cb3 ZAL	17	/*
	18	* copysign(x,y),
	19	* logb(x),
	20	* scalb(x,N),
	21	* finite(x),
	22	* drem(x,y),
	23	* Coded in vax assembly language by K. C. Ng 4/9/85.
	24	* Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87.
	25	*/
	26	/*
	27	* double copysign(x,y)
	28	* double x,y;
	29	*/
	30	.globl _copysign
	31	.text
f47ef219	32	.align 2
97158cb3 ZAL	33	_copysign:
	34	.word 0x0004 # save r2
	35	movl 8(fp),r1
	36	movl 4(fp),r0 # r0:r1 = x
	37	andl3 $0x7f800000,r0,r2 # r2 = biased exponent of x
	38	beql 1f # if 0 or reserved op then return x
	39	andl3 $0x80000000,12(fp),r2 # r2 = sign bit of y at bit-31
	40	andl2 $0x7fffffff,r0 # replace x by \|x\|
	41	orl2 r2,r0 # copy the sign bit of y to x
	42	1: ret
	43	/*
	44	* double logb(x)
	45	* double x;
	46	*/
	47	.globl _logb
	48	.text
f47ef219	49	.align 2
97158cb3	50	_logb:
f47ef219	51	.word 0x0000 # save nothing
97158cb3 ZAL	52	andl3 $0x7f800000,4(fp),r0 # r0[b23:b30] = biased exponent of x
	53	beql 1f
	54	shrl $23,r0,r0 # r0[b0:b7] = biased exponent of x
	55	subl2 $129,r0 # r0 = unbiased exponent of x
	56	cvld r0 # acc = unbiased exponent of x (double)
	57	std r0 # r0 = unbiased exponent of x (double)
	58	ret
	59	1: movl 8(fp),r1 # 8(fp) must be moved first
	60	movl 4(fp),r0 # r0:r1 = x (zero or reserved op)
	61	blss 2f # simply return if reserved op
	62	movl $0xfe000000,r1
	63	movl $0xcfffffff,r0 # -2147483647.0
	64	2: ret
	65	/*
	66	* long finite(x)
	67	* double x;
	68	*/
	69	.globl _finite
	70	.text
f47ef219	71	.align 2
97158cb3	72	_finite:
f47ef219	73	.word 0x0000 # save nothing
97158cb3 ZAL	74	andl3 $0xff800000,4(fp),r0 # r0 = sign of x & its biased exponent
	75	cmpl r0,$0x80000000 # is x a reserved op?
	76	beql 1f # if so, return FALSE (0)
	77	movl $1,r0 # else return TRUE (1)
	78	ret
	79	1: clrl r0
	80	ret
	81	/*
	82	* double scalb(x,N)
	83	* double x; int N;
	84	*/
	85	.globl _scalb
	86	.set ERANGE,34
	87	.text
f47ef219	88	.align 2
97158cb3	89	_scalb:
f47ef219	90	.word 0x000c # save r2-r3
97158cb3 ZAL	91	movl 8(fp),r1
	92	movl 4(fp),r0 # r0:r1 = x (-128 <= Ex <= 126)
	93	andl3 $0x7f800000,r0,r3 # r3[b23:b30] = biased exponent of x
	94	beql 1f # is x a 0 or a reserved operand?
	95	movl 12(fp),r2 # r2 = N
	96	cmpl r2,$0xff # if N >= 255
	97	bgeq 2f # then the result must overflow
	98	cmpl r2,$-0xff # if N <= -255
	99	bleq 3f # then the result must underflow
	100	shrl $23,r3,r3 # r3[b0:b7] = biased exponent of x
	101	addl2 r2,r3 # r3 = biased exponent of the result
	102	bleq 3f # if <= 0 then the result underflows
	103	cmpl r3,$0x100 # if >= 256 then the result overflows
	104	bgeq 2f
	105	shll $23,r3,r3 # r3[b23:b30] = biased exponent of res.
	106	andl2 $0x807fffff,r0
	107	orl2 r3,r0 # r0:r1 = x*2^N
	108	1: ret
	109	2: pushl $ERANGE # if the result would overflow
	110	callf $8,_infnan # and _infnan returns
	111	andl3 $0x80000000,4(fp),r2 # get the sign of input arg
	112	orl2 r2,r0 # re-attach the sign to r0:r1
	113	ret
	114	3: clrl r1 # if the result would underflow
	115	clrl r0 # then return 0
	116	ret
	117	/*
	118	* double drem(x,y)
	119	* double x,y;
	120	* Returns x-n*y where n=[x/y] rounded (to even in the half way case).
	121	*/
	122	.globl _drem
	123	.set EDOM,33
	124	.text
f47ef219	125	.align 2
97158cb3 ZAL	126	_drem:
	127	.word 0x1ffc # save r2-r12
	128	movl 16(fp),r3
	129	movl 12(fp),r2 # r2:r3 = y
	130	movl 8(fp),r1
	131	movl 4(fp),r0 # r0:r1 = x
	132	andl3 $0xff800000,r0,r4
	133	cmpl r4,$0x80000000 # is x a reserved operand?
	134	beql 1f # if yes then propagate x and return
	135	andl3 $0xff800000,r2,r4
	136	cmpl r4,$0x80000000 # is y a reserved operand?
	137	bneq 2f
	138	movl r3,r1
	139	movl r2,r0 # if yes then propagate y and return
	140	1: ret
	141
	142	2: tstl r4 # is y a 0?
	143	bneq 3f
	144	pushl $EDOM # if so then generate reserved op fault
	145	callf $8,_infnan
	146	ret
	147
	148	3: andl2 $0x7fffffff,r2 # r2:r3 = y <- \|y\|
	149	clrl r12 # r12 = nx := 0
	150	cmpl r2,$0x1c800000 # Ey ? 57
	151	bgtr 4f # if Ey > 57 goto 4
	152	addl2 $0x1c800000,r2 # scale up y by 2**57
	153	movl $0x1c800000,r12 # r12[b23:b30] = nx = 57
	154	4: pushl r12 # pushed onto stack: nf := nx
	155	andl3 $0x80000000,r0,-(sp) # pushed onto stack: sign of x
	156	andl2 $0x7fffffff,r0 # r0:r1 = x <- \|x\|
	157	movl r3,r11 # r10:r11 = y1 = y w/ last 27 bits 0
	158	andl3 $0xf8000000,r10,r11 # clear last 27 bits of y1
	159
	160	Loop: cmpd2 r0,r2 # x ? y
	161	bleq 6f # if x <= y goto 6
	162	/* # begin argument reduction */
	163	movl r3,r5
	164	movl r2,r4 # r4:r5 = t = y
	165	movl r11,r7
	166	movl r10,r6 # r6:r7 = t1 = y1
	167	andl3 $0x7f800000,r0,r8 # r8[b23:b30] = Ex:biased exponent of x
	168	andl3 $0x7f800000,r2,r9 # r9[b23:b30] = Ey:biased exponent of y
	169	subl2 r9,r8 # r8[b23:b30] = Ex-Ey
	170	subl2 $0x0c800000,r8 # r8[b23:b30] = k = Ex-Ey-25
	171	blss 5f # if k < 0 goto 5
	172	addl2 r8,r4 # t += k
	173	addl2 r8,r6 # t1 += k, scale up t and t1
	174	5: ldd r0 # acc = x
	175	divd r4 # acc = x/t
	176	cvdl r8 # r8 = n = [x/t] truncated
	177	cvld r8 # acc = dble(n)
	178	std r8 # r8:r9 = dble(n)
	179	ldd r4 # acc = t
	180	subd r6 # acc = t-t1
	181	muld r8 # acc = n*(t-t1)
	182	std r4 # r4:r5 = n*(t-t1)
	183	ldd r6 # acc = t1
	184	muld r8 # acc = n*t1
	185	subd r0 # acc = n*t1-x
	186	negd # acc = x-n*t1
	187	subd r4 # acc = (x-nt1)-n(t-t1)
	188	std r0 # r0:r1 = (x-nt1)-n(t-t1)
	189	brb Loop
190
191	6: movl r12,r6 # r6 = nx
192	beql 7f # if nx == 0 goto 7
193	addl2 r6,r0 # x <- x2*57:scale x up by nx
194	clrl r12 # clear nx
195	brb Loop
196
197	7: movl r3,r5
198	movl r2,r4 # r4:r5 = y
199	subl2 $0x800000,r4 # r4:r5 = y/2
200	cmpd2 r0,r4 # x ? y/2
201	blss 9f # if x < y/2 goto 9
202	bgtr 8f # if x > y/2 goto 8
203	ldd r8 # acc = dble(n)
204	cvdl r8 # r8 = ifix(dble(n))
205	bbc $0,r8,9f # if the last bit is zero, goto 9
206	8: ldd r0 # acc = x
207	subd r2 # acc = x-y
208	std r0 # r0:r1 = x-y
209	9: xorl2 (sp)+,r0 # x^sign (exclusive or)
210	movl (sp)+,r6 # r6 = nf
211	andl3 $0x7f800000,r0,r8 # r8 = biased exponent of x
212	andl2 $0x807fffff,r0 # r0 = x w/ exponent zapped
213	subl2 r6,r8 # r8 = Ex-nf
214	bgtr 0f # if Ex-nf > 0 goto 0
215	clrl r8 # clear r8
216	clrl r0
217	clrl r1 # x underflows to zero
218	0: orl2 r8,r0 # put r8 into x's exponent field
219	ret