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 #define n0 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, sx; 43 44 if (!finitel(x)) 45 return (x + x); 46 if (iszerol(x)) 47 return (x); 48 sx = signbitl(x); 49 x = fabsl(x); 50 n = ilogbl(x); 51 m = n / 3; 52 m3 = m + m + m; 53 y = scalbnl(x, -m3); 54 dx = (double) y; 55 dy = cbrt(dx); 56 py[1 - n0] += 2; 57 if (py[1 - n0] == 0) 58 py[n0] += 1; 59 60 /* one step newton iteration to 113 bits with error < 0.667ulps */ 61 t = (long double) dy; 62 t = scalbnl(t, m); 63 s = t * t; 64 r = x / s; 65 w = t + t; 66 r = (r - t) / (w + r); 67 t += t * r; 68 69 return (sx == 0 ? t : -t); 70 } 71