xref: /illumos-gate/usr/src/lib/libm/common/m9x/frexp.c (revision 25c28e83beb90e7c80452a7c818c5e6f73a07dc8)
1 /*
2  * CDDL HEADER START
3  *
4  * The contents of this file are subject to the terms of the
5  * Common Development and Distribution License (the "License").
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
17  * information: Portions Copyright [yyyy] [name of copyright owner]
18  *
19  * CDDL HEADER END
20  */
21 
22 /*
23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
24  */
25 /*
26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27  * Use is subject to license terms.
28  */
29 
30 #if defined(ELFOBJ)
31 #pragma weak frexp = __frexp
32 #endif
33 
34 /*
35  * frexp(x, exp) returns the normalized significand of x and sets
36  * *exp so that x = r*2^(*exp) where r is the return value.  If x
37  * is finite and nonzero, 1/2 <= |r| < 1.
38  *
39  * If x is zero, infinite or NaN, frexp returns x and sets *exp = 0.
40  * (The relevant standards do not specify *exp when x is infinite or
41  * NaN, but this code sets it anyway.)
42  *
43  * If x is a signaling NaN, this code returns x without attempting
44  * to raise the invalid operation exception.  If x is subnormal,
45  * this code treats it as nonzero regardless of nonstandard mode.
46  */
47 
48 #include "libm.h"
49 
50 double
51 __frexp(double x, int *exp) {
52 	union {
53 		unsigned i[2];
54 		double d;
55 	} xx, yy;
56 	double t;
57 	unsigned hx;
58 	int e;
59 
60 	xx.d = x;
61 	hx = xx.i[HIWORD] & ~0x80000000;
62 
63 	if (hx >= 0x7ff00000) { /* x is infinite or NaN */
64 		*exp = 0;
65 		return (x);
66 	}
67 
68 	e = 0;
69 	if (hx < 0x00100000) { /* x is subnormal or zero */
70 		if ((hx | xx.i[LOWORD]) == 0) {
71 			*exp = 0;
72 			return (x);
73 		}
74 
75 		/*
76 		 * normalize x by regarding it as an integer
77 		 *
78 		 * Here we use 32-bit integer arithmetic to avoid trapping
79 		 * or emulating 64-bit arithmetic.  If 64-bit arithmetic is
80 		 * available (e.g., in SPARC V9), do this instead:
81 		 *
82 		 *  long lx = ((long) hx << 32) | xx.i[LOWORD];
83 		 *  xx.d = (xx.i[HIWORD] < 0)? -lx : lx;
84 		 *
85 		 * If subnormal arithmetic doesn't trap, just multiply x by
86 		 * a power of two.
87 		 */
88 		yy.i[HIWORD] = 0x43300000 | hx;
89 		yy.i[LOWORD] = xx.i[LOWORD];
90 		t = yy.d;
91 		yy.i[HIWORD] = 0x43300000;
92 		yy.i[LOWORD] = 0;
93 		t -= yy.d; /* t = |x| scaled */
94 		xx.d = ((int)xx.i[HIWORD] < 0)? -t : t;
95 		hx = xx.i[HIWORD] & ~0x80000000;
96 		e = -1074;
97 	}
98 
99 	/* now xx.d is normal */
100 	xx.i[HIWORD] = (xx.i[HIWORD] & ~0x7ff00000) | 0x3fe00000;
101 	*exp = e + (hx >> 20) - 0x3fe;
102 	return (xx.d);
103 }
104