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 cbrtl = __cbrtl 32 #endif 33 34 #include "libm.h" 35 #include "longdouble.h" 36 37 static const double d_one = 1.0; 38 39 long double 40 cbrtl(long double x) { 41 long double s, t, r, w, y; 42 double dx, dy; 43 int *py = (int *) &dy; 44 int n, m, m3, n0, sx; 45 46 if (!finitel(x)) 47 return (x + x); 48 if (iszerol(x)) 49 return (x); 50 n0 = 0; 51 if (*((int *) &d_one) == 0) 52 n0 = 1; 53 sx = signbitl(x); 54 x = fabsl(x); 55 n = ilogbl(x); 56 m = n / 3; 57 m3 = m + m + m; 58 y = scalbnl(x, -m3); 59 dx = (double) y; 60 dy = cbrt(dx); 61 py[1 - n0] += 2; 62 if (py[1 - n0] == 0) 63 py[n0] += 1; 64 65 /* one step newton iteration to 113 bits with error < 0.667ulps */ 66 t = (long double) dy; 67 t = scalbnl(t, m); 68 s = t * t; 69 r = x / s; 70 w = t + t; 71 r = (r - t) / (w + r); 72 t += t * r; 73 74 return (sx == 0 ? t : -t); 75 } 76