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