xref: /netbsd/lib/libc/arch/mips/gen/ldexp.S (revision 6550d01e) 
1/*	$NetBSD: ldexp.S,v 1.9 2009/12/14 01:07:42 matt 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 <mips/asm.h>
36
37#if defined(LIBC_SCCS) && !defined(lint)
38#if 0
39	RCSID("from: @(#)ldexp.s	8.1 (Berkeley) 6/4/93")
40#else
41	RCSID("$NetBSD: ldexp.S,v 1.9 2009/12/14 01:07:42 matt Exp $")
42#endif
43#endif /* LIBC_SCCS and not lint */
44
45#define DEXP_INF	0x7ff
46#define DEXP_BIAS	1023
47#define DEXP_MIN	-1022
48#define DEXP_MAX	1023
49#define DFRAC_BITS	52
50#define DIMPL_ONE	0x00100000
51#define DLEAD_ZEROS	31 - 20
52#define STICKYBIT	1
53#define GUARDBIT	0x80000000
54#define DSIGNAL_NAN	0x00040000
55#define DQUIET_NAN0	0x0007ffff
56#define DQUIET_NAN1	0xffffffff
57
58/*
59 * double ldexp(x, N)
60 *	double x; int N;
61 *
62 * Return x * (2**N), for integer values N.
63 */
64LEAF(ldexp)
65	mfc1	v1, $f13		# get MSW of x
66	mfc1	t3, $f12		# get LSW of x
67	sll	t1, v1, 1		# get x exponent
68	srl	t1, t1, 32 - 11
69	beq	t1, DEXP_INF, 9f	# is it a NAN or infinity?
70	beq	t1, zero, 1f		# zero or denormalized number?
71	addu	t1, t1, a2		# scale exponent
72	sll	v0, a2, 20		# position N for addition
73	bge	t1, DEXP_INF, 8f	# overflow?
74	addu	v0, v0, v1		# multiply by (2**N)
75	ble	t1, zero, 4f		# underflow?
76	mtc1	v0, $f1			# save MSW of result
77	mtc1	t3, $f0			# save LSW of result
78	j	ra
791:
80	sll	t2, v1, 32 - 20		# get x fraction
81	srl	t2, t2, 32 - 20
82	srl	t0, v1, 31		# get x sign
83	bne	t2, zero, 1f
84	beq	t3, zero, 9f		# result is zero
851:
86/*
87 * Find out how many leading zero bits are in t2,t3 and put in t9.
88 */
89	move	v0, t2
90	move	t9, zero
91	bne	t2, zero, 1f
92	move	v0, t3
93	addu	t9, 32
941:
95	srl	ta0, v0, 16
96	bne	ta0, zero, 1f
97	addu	t9, 16
98	sll	v0, 16
991:
100	srl	ta0, v0, 24
101	bne	ta0, zero, 1f
102	addu	t9, 8
103	sll	v0, 8
1041:
105	srl	ta0, v0, 28
106	bne	ta0, zero, 1f
107	addu	t9, 4
108	sll	v0, 4
1091:
110	srl	ta0, v0, 30
111	bne	ta0, zero, 1f
112	addu	t9, 2
113	sll	v0, 2
1141:
115	srl	ta0, v0, 31
116	bne	ta0, zero, 1f
117	addu	t9, 1
118/*
119 * Now shift t2,t3 the correct number of bits.
120 */
1211:
122	subu	t9, t9, DLEAD_ZEROS	# dont count normal leading zeros
123	li	t1, DEXP_MIN + DEXP_BIAS
124	subu	t1, t1, t9		# adjust exponent
125	addu	t1, t1, a2		# scale exponent
126	li	v0, 32
127	blt	t9, v0, 1f
128	subu	t9, t9, v0		# shift fraction left >= 32 bits
129	sll	t2, t3, t9
130	move	t3, zero
131	b	2f
1321:
133	subu	v0, v0, t9		# shift fraction left < 32 bits
134	sll	t2, t2, t9
135	srl	ta0, t3, v0
136	or	t2, t2, ta0
137	sll	t3, t3, t9
1382:
139	bge	t1, DEXP_INF, 8f	# overflow?
140	ble	t1, zero, 4f		# underflow?
141	sll	t2, t2, 32 - 20		# clear implied one bit
142	srl	t2, t2, 32 - 20
1433:
144	sll	t1, t1, 31 - 11		# reposition exponent
145	sll	t0, t0, 31		# reposition sign
146	or	t0, t0, t1		# put result back together
147	or	t0, t0, t2
148	mtc1	t0, $f1			# save MSW of result
149	mtc1	t3, $f0			# save LSW of result
150	j	ra
1514:
152	li	v0, 0x80000000
153	ble	t1, -52, 7f		# is result too small for denorm?
154	sll	t2, v1, 31 - 20		# clear exponent, extract fraction
155	or	t2, t2, v0		# set implied one bit
156	blt	t1, -30, 2f		# will all bits in t3 be shifted out?
157	srl	t2, t2, 31 - 20		# shift fraction back to normal position
158	subu	t1, t1, 1
159	sll	ta0, t2, t1		# shift right t2,t3 based on exponent
160	srl	t8, t3, t1		# save bits shifted out
161	negu	t1
162	srl	t3, t3, t1
163	or	t3, t3, ta0
164	srl	t2, t2, t1
165	bge	t8, zero, 1f		# does result need to be rounded?
166	addu	t3, t3, 1		# round result
167	sltu	ta0, t3, 1
168	sll	t8, t8, 1
169	addu	t2, t2, ta0
170	bne	t8, zero, 1f		# round result to nearest
171	and	t3, t3, ~1
1721:
173	mtc1	t3, $f0			# save denormalized result (LSW)
174	mtc1	t2, $f1			# save denormalized result (MSW)
175	bge	v1, zero, 1f		# should result be negative?
176	neg.d	$f0, $f0		# negate result
1771:
178	j	ra
1792:
180	mtc1	zero, $f1		# exponent and upper fraction
181	addu	t1, t1, 20		# compute amount to shift right by
182	sll	t8, t2, t1		# save bits shifted out
183	negu	t1
184	srl	t3, t2, t1
185	bge	t8, zero, 1f		# does result need to be rounded?
186	addu	t3, t3, 1		# round result
187	sltu	ta0, t3, 1
188	sll	t8, t8, 1
189	mtc1	ta0, $f1			# exponent and upper fraction
190	bne	t8, zero, 1f		# round result to nearest
191	and	t3, t3, ~1
1921:
193	mtc1	t3, $f0
194	bge	v1, zero, 1f		# is result negative?
195	neg.d	$f0, $f0		# negate result
1961:
197	j	ra
1987:
199	mtc1	zero, $f0		# result is zero
200	mtc1	zero, $f1
201	beq	t0, zero, 1f		# is result positive?
202	neg.d	$f0, $f0		# negate result
2031:
204	j	ra
2058:
206	li	t1, 0x7ff00000		# result is infinity (MSW)
207	mtc1	t1, $f1
208	mtc1	zero, $f0		# result is infinity (LSW)
209	bge	v1, zero, 1f		# should result be negative infinity?
210	neg.d	$f0, $f0		# result is negative infinity
2111:
212	add.d	$f0, $f0		# cause overflow faults if enabled
213	j	ra
2149:
215	mov.d	$f0, $f12		# yes, result is just x
216	j	ra
217END(ldexp)
218