xref: /csrg-svn/lib/libm/tahoe/support.s (revision 31842) 
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:
16	.asciz	"@(#)support.s	1.2	(ucb.elefunt)	07/13/87"
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
32	.align	2
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
421:	ret
43/*
44 * double logb(x)
45 * double x;
46 */
47	.globl	_logb
48	.text
49	.align	2
50_logb:
51	.word	0x0000			# save nothing
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
591:	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
642:	ret
65/*
66 * long finite(x)
67 * double x;
68 */
69	.globl	_finite
70	.text
71	.align	2
72_finite:
73	.word	0x0000			# save nothing
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
791:	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
88	.align	2
89_scalb:
90	.word	0x000c			# save r2-r3
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
1081:	ret
1092:	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
1143:	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
125	.align	2
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
1401:	ret
141
1422:	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
1483:	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
1544:	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
160Loop:	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
1745:	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-n*t1)-n*(t-t1)
188	std	r0			# r0:r1 = (x-n*t1)-n*(t-t1)
189	brb	Loop
190
1916:	movl	r12,r6			# r6 = nx
192	beql	7f			# if nx == 0 goto 7
193	addl2	r6,r0			# x <- x*2**57:scale x up by nx
194	clrl	r12			# clear nx
195	brb	Loop
196
1977:	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
2068:	ldd	r0			# acc = x
207	subd	r2			# acc = x-y
208	std	r0			# r0:r1 = x-y
2099:	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
2180:	orl2	r8,r0			# put r8 into x's exponent field
219	ret
220