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 scalbln = __scalbln 32 #endif 33 34 #include "libm.h" 35 #include <float.h> /* DBL_MAX, DBL_MIN */ 36 37 static const double twom54 = 5.5511151231257827021181583404541015625e-17; 38 #if defined(USE_FPSCALE) || defined(__x86) 39 static const double two52 = 4503599627370496.0; 40 #else 41 /* 42 * Normalize non-zero subnormal x and return biased exponent of x in [-51,0] 43 */ 44 static int 45 ilogb_biased(unsigned *px) { 46 int s = 52; 47 unsigned v = px[HIWORD] & ~0x80000000, w = px[LOWORD], t = v; 48 49 if (t) 50 s -= 32; 51 else 52 t = w; 53 if (t & 0xffff0000) 54 s -= 16, t >>= 16; 55 if (t & 0xff00) 56 s -= 8, t >>= 8; 57 if (t & 0xf0) 58 s -= 4, t >>= 4; 59 t <<= 1; 60 s -= (0xffffaa50 >> t) & 0x3; 61 if (s < 32) { 62 v = (v << s) | w >> (32 - s); 63 w <<= s; 64 } else { 65 v = w << (s - 32); 66 w = 0; 67 } 68 px[HIWORD] = (px[HIWORD] & 0x80000000) | v; 69 px[LOWORD] = w; 70 return (1 - s); 71 } 72 #endif /* defined(USE_FPSCALE) */ 73 74 double 75 scalbln(double x, long n) { 76 int *px = (int *) &x, ix, k; 77 78 ix = px[HIWORD] & ~0x80000000; 79 k = ix >> 20; 80 if (k == 0x7ff) 81 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN) 82 return ((px[HIWORD] & 0x80000) != 0 ? x : x + x); 83 /* assumes sparc-like QNaN */ 84 #else 85 return (x + x); 86 #endif 87 if ((px[LOWORD] | ix) == 0 || n == 0) 88 return (x); 89 if (k == 0) { 90 #if defined(USE_FPSCALE) || defined(__x86) 91 x *= two52; 92 k = ((px[HIWORD] & ~0x80000000) >> 20) - 52; 93 #else 94 k = ilogb_biased((unsigned *) px); 95 #endif 96 } 97 k += (int) n; 98 if (n > 5000 || k > 0x7fe) 99 return (DBL_MAX * copysign(DBL_MAX, x)); 100 if (n < -5000 || k <= -54) 101 return (DBL_MIN * copysign(DBL_MIN, x)); 102 if (k > 0) { 103 px[HIWORD] = (px[HIWORD] & ~0x7ff00000) | (k << 20); 104 return (x); 105 } 106 k += 54; 107 px[HIWORD] = (px[HIWORD] & ~0x7ff00000) | (k << 20); 108 return (x * twom54); 109 } 110