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