1/*	$NetBSD: ldexp.S,v 1.8 2003/08/07 16:42:15 agc Exp $	*/
2
3/*-
4 * Copyright (c) 1991, 1993
5 *	The Regents of the University of California.  All rights reserved.
6 *
7 * This code is derived from software contributed to Berkeley by
8 * Ralph Campbell.
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 * 1. Redistributions of source code must retain the above copyright
14 *    notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 *    notice, this list of conditions and the following disclaimer in the
17 *    documentation and/or other materials provided with the distribution.
18 * 3. Neither the name of the University nor the names of its contributors
19 *    may be used to endorse or promote products derived from this software
20 *    without specific prior written permission.
21 *
22 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
23 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
26 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
28 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
29 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
30 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
31 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
32 * SUCH DAMAGE.
33 */
34
35#include <machine/asm.h>
36__FBSDID("$FreeBSD$");
37
38#if defined(LIBC_SCCS) && !defined(lint)
39	ASMSTR("from: @(#)ldexp.s	8.1 (Berkeley) 6/4/93")
40	ASMSTR("$NetBSD: ldexp.S,v 1.8 2003/08/07 16:42:15 agc Exp $")
41#endif /* LIBC_SCCS and not lint */
42
43#ifdef __ABICALLS__
44	.abicalls
45#endif
46
47#define DEXP_INF	0x7ff
48#define DEXP_BIAS	1023
49#define DEXP_MIN	-1022
50#define DEXP_MAX	1023
51#define DFRAC_BITS	52
52#define DIMPL_ONE	0x00100000
53#define DLEAD_ZEROS	31 - 20
54#define STICKYBIT	1
55#define GUARDBIT	0x80000000
56#define DSIGNAL_NAN	0x00040000
57#define DQUIET_NAN0	0x0007ffff
58#define DQUIET_NAN1	0xffffffff
59
60/*
61 * double ldexp(x, N)
62 *	double x; int N;
63 *
64 * Return x * (2**N), for integer values N.
65 */
66LEAF(ldexp)
67	mfc1	v1, $f13		# get MSW of x
68	mfc1	t3, $f12		# get LSW of x
69	sll	t1, v1, 1		# get x exponent
70	srl	t1, t1, 32 - 11
71	beq	t1, DEXP_INF, 9f	# is it a NAN or infinity?
72	beq	t1, zero, 1f		# zero or denormalized number?
73	addu	t1, t1, a2		# scale exponent
74	sll	v0, a2, 20		# position N for addition
75	bge	t1, DEXP_INF, 8f	# overflow?
76	addu	v0, v0, v1		# multiply by (2**N)
77	ble	t1, zero, 4f		# underflow?
78	mtc1	v0, $f1			# save MSW of result
79	mtc1	t3, $f0			# save LSW of result
80	j	ra
811:
82	sll	t2, v1, 32 - 20		# get x fraction
83	srl	t2, t2, 32 - 20
84	srl	t0, v1, 31		# get x sign
85	bne	t2, zero, 1f
86	beq	t3, zero, 9f		# result is zero
871:
88/*
89 * Find out how many leading zero bits are in t2,t3 and put in t9.
90 */
91	move	v0, t2
92	move	t9, zero
93	bne	t2, zero, 1f
94	move	v0, t3
95	addu	t9, 32
961:
97	srl	ta0, v0, 16
98	bne	ta0, zero, 1f
99	addu	t9, 16
100	sll	v0, 16
1011:
102	srl	ta0, v0, 24
103	bne	ta0, zero, 1f
104	addu	t9, 8
105	sll	v0, 8
1061:
107	srl	ta0, v0, 28
108	bne	ta0, zero, 1f
109	addu	t9, 4
110	sll	v0, 4
1111:
112	srl	ta0, v0, 30
113	bne	ta0, zero, 1f
114	addu	t9, 2
115	sll	v0, 2
1161:
117	srl	ta0, v0, 31
118	bne	ta0, zero, 1f
119	addu	t9, 1
120/*
121 * Now shift t2,t3 the correct number of bits.
122 */
1231:
124	subu	t9, t9, DLEAD_ZEROS	# dont count normal leading zeros
125	li	t1, DEXP_MIN + DEXP_BIAS
126	subu	t1, t1, t9		# adjust exponent
127	addu	t1, t1, a2		# scale exponent
128	li	v0, 32
129	blt	t9, v0, 1f
130	subu	t9, t9, v0		# shift fraction left >= 32 bits
131	sll	t2, t3, t9
132	move	t3, zero
133	b	2f
1341:
135	subu	v0, v0, t9		# shift fraction left < 32 bits
136	sll	t2, t2, t9
137	srl	ta0, t3, v0
138	or	t2, t2, ta0
139	sll	t3, t3, t9
1402:
141	bge	t1, DEXP_INF, 8f	# overflow?
142	ble	t1, zero, 4f		# underflow?
143	sll	t2, t2, 32 - 20		# clear implied one bit
144	srl	t2, t2, 32 - 20
1453:
146	sll	t1, t1, 31 - 11		# reposition exponent
147	sll	t0, t0, 31		# reposition sign
148	or	t0, t0, t1		# put result back together
149	or	t0, t0, t2
150	mtc1	t0, $f1			# save MSW of result
151	mtc1	t3, $f0			# save LSW of result
152	j	ra
1534:
154	li	v0, 0x80000000
155	ble	t1, -52, 7f		# is result too small for denorm?
156	sll	t2, v1, 31 - 20		# clear exponent, extract fraction
157	or	t2, t2, v0		# set implied one bit
158	blt	t1, -30, 2f		# will all bits in t3 be shifted out?
159	srl	t2, t2, 31 - 20		# shift fraction back to normal position
160	subu	t1, t1, 1
161	sll	ta0, t2, t1		# shift right t2,t3 based on exponent
162	srl	t8, t3, t1		# save bits shifted out
163	negu	t1
164	srl	t3, t3, t1
165	or	t3, t3, ta0
166	srl	t2, t2, t1
167	bge	t8, zero, 1f		# does result need to be rounded?
168	addu	t3, t3, 1		# round result
169	sltu	ta0, t3, 1
170	sll	t8, t8, 1
171	addu	t2, t2, ta0
172	bne	t8, zero, 1f		# round result to nearest
173	and	t3, t3, ~1
1741:
175	mtc1	t3, $f0			# save denormalized result (LSW)
176	mtc1	t2, $f1			# save denormalized result (MSW)
177	bge	v1, zero, 1f		# should result be negative?
178	neg.d	$f0, $f0		# negate result
1791:
180	j	ra
1812:
182	mtc1	zero, $f1		# exponent and upper fraction
183	addu	t1, t1, 20		# compute amount to shift right by
184	sll	t8, t2, t1		# save bits shifted out
185	negu	t1
186	srl	t3, t2, t1
187	bge	t8, zero, 1f		# does result need to be rounded?
188	addu	t3, t3, 1		# round result
189	sltu	ta0, t3, 1
190	sll	t8, t8, 1
191	mtc1	ta0, $f1			# exponent and upper fraction
192	bne	t8, zero, 1f		# round result to nearest
193	and	t3, t3, ~1
1941:
195	mtc1	t3, $f0
196	bge	v1, zero, 1f		# is result negative?
197	neg.d	$f0, $f0		# negate result
1981:
199	j	ra
2007:
201	mtc1	zero, $f0		# result is zero
202	mtc1	zero, $f1
203	beq	t0, zero, 1f		# is result positive?
204	neg.d	$f0, $f0		# negate result
2051:
206	j	ra
2078:
208	li	t1, 0x7ff00000		# result is infinity (MSW)
209	mtc1	t1, $f1
210	mtc1	zero, $f0		# result is infinity (LSW)
211	bge	v1, zero, 1f		# should result be negative infinity?
212	neg.d	$f0, $f0		# result is negative infinity
2131:
214	add.d	$f0, $f0		# cause overflow faults if enabled
215	j	ra
2169:
217	mov.d	$f0, $f12		# yes, result is just x
218	j	ra
219END(ldexp)
220