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