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 #pragma weak frexp = __frexp 31 32 /* 33 * frexp(x, exp) returns the normalized significand of x and sets 34 * *exp so that x = r*2^(*exp) where r is the return value. If x 35 * is finite and nonzero, 1/2 <= |r| < 1. 36 * 37 * If x is zero, infinite or NaN, frexp returns x and sets *exp = 0. 38 * (The relevant standards do not specify *exp when x is infinite or 39 * NaN, but this code sets it anyway.) 40 * 41 * If x is a signaling NaN, this code returns x without attempting 42 * to raise the invalid operation exception. If x is subnormal, 43 * this code treats it as nonzero regardless of nonstandard mode. 44 */ 45 46 #include "libm.h" 47 48 double 49 __frexp(double x, int *exp) { 50 union { 51 unsigned i[2]; 52 double d; 53 } xx, yy; 54 double t; 55 unsigned hx; 56 int e; 57 58 xx.d = x; 59 hx = xx.i[HIWORD] & ~0x80000000; 60 61 if (hx >= 0x7ff00000) { /* x is infinite or NaN */ 62 *exp = 0; 63 return (x); 64 } 65 66 e = 0; 67 if (hx < 0x00100000) { /* x is subnormal or zero */ 68 if ((hx | xx.i[LOWORD]) == 0) { 69 *exp = 0; 70 return (x); 71 } 72 73 /* 74 * normalize x by regarding it as an integer 75 * 76 * Here we use 32-bit integer arithmetic to avoid trapping 77 * or emulating 64-bit arithmetic. If 64-bit arithmetic is 78 * available (e.g., in SPARC V9), do this instead: 79 * 80 * long lx = ((long) hx << 32) | xx.i[LOWORD]; 81 * xx.d = (xx.i[HIWORD] < 0)? -lx : lx; 82 * 83 * If subnormal arithmetic doesn't trap, just multiply x by 84 * a power of two. 85 */ 86 yy.i[HIWORD] = 0x43300000 | hx; 87 yy.i[LOWORD] = xx.i[LOWORD]; 88 t = yy.d; 89 yy.i[HIWORD] = 0x43300000; 90 yy.i[LOWORD] = 0; 91 t -= yy.d; /* t = |x| scaled */ 92 xx.d = ((int)xx.i[HIWORD] < 0)? -t : t; 93 hx = xx.i[HIWORD] & ~0x80000000; 94 e = -1074; 95 } 96 97 /* now xx.d is normal */ 98 xx.i[HIWORD] = (xx.i[HIWORD] & ~0x7ff00000) | 0x3fe00000; 99 *exp = e + (hx >> 20) - 0x3fe; 100 return (xx.d); 101 } 102