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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak sqrtf = __sqrtf 30 31 #include "libm.h" 32 33 #ifdef __INLINE 34 35 extern float __inline_sqrtf(float); 36 37 float 38 sqrtf(float x) { 39 return (__inline_sqrtf(x)); 40 } 41 42 #else /* defined(__INLINE) */ 43 44 static const float huge = 1.0e35F, tiny = 1.0e-35F, zero = 0.0f; 45 46 float 47 sqrtf(float x) { 48 float dz, w; 49 int *pw = (int *)&w; 50 int ix, j, r, q, m, n, s, t; 51 52 w = x; 53 ix = pw[0]; 54 if (ix <= 0) { 55 /* x is <= 0 or nan */ 56 j = ix & 0x7fffffff; 57 if (j == 0) 58 return (w); 59 return ((w * zero) / zero); 60 } 61 62 if ((ix & 0x7f800000) == 0x7f800000) { 63 /* x is +inf or nan */ 64 return (w * w); 65 } 66 67 m = ir_ilogb_(&w); 68 n = -m; 69 w = r_scalbn_(&w, (int *)&n); 70 ix = (pw[0] & 0x007fffff) | 0x00800000; 71 n = m / 2; 72 if ((n + n) != m) { 73 ix = ix + ix; 74 m -= 1; 75 n = m / 2; 76 } 77 78 /* generate sqrt(x) bit by bit */ 79 ix <<= 1; 80 q = s = 0; 81 r = 0x01000000; 82 for (j = 1; j <= 25; j++) { 83 t = s + r; 84 if (t <= ix) { 85 s = t + r; 86 ix -= t; 87 q += r; 88 } 89 ix <<= 1; 90 r >>= 1; 91 } 92 if (ix == 0) 93 goto done; 94 95 /* raise inexact and determine the ambient rounding mode */ 96 dz = huge - tiny; 97 if (dz < huge) 98 goto done; 99 dz = huge + tiny; 100 if (dz > huge) 101 q += 1; 102 q += (q & 1); 103 104 done: 105 pw[0] = (q >> 1) + 0x3f000000; 106 return (r_scalbn_(&w, (int *)&n)); 107 } 108 109 #endif /* defined(__INLINE) */ 110