diff options
Diffstat (limited to 'third_party/js-1.7/fdlibm/s_cbrt.c')
| -rw-r--r-- | third_party/js-1.7/fdlibm/s_cbrt.c | 133 |
1 files changed, 133 insertions, 0 deletions
diff --git a/third_party/js-1.7/fdlibm/s_cbrt.c b/third_party/js-1.7/fdlibm/s_cbrt.c new file mode 100644 index 0000000..4aed19b --- /dev/null +++ b/third_party/js-1.7/fdlibm/s_cbrt.c @@ -0,0 +1,133 @@ +/* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*- + * + * ***** BEGIN LICENSE BLOCK ***** + * Version: MPL 1.1/GPL 2.0/LGPL 2.1 + * + * The contents of this file are subject to the Mozilla Public License Version + * 1.1 (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * Software distributed under the License is distributed on an "AS IS" basis, + * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License + * for the specific language governing rights and limitations under the + * License. + * + * The Original Code is Mozilla Communicator client code, released + * March 31, 1998. + * + * The Initial Developer of the Original Code is + * Sun Microsystems, Inc. + * Portions created by the Initial Developer are Copyright (C) 1998 + * the Initial Developer. All Rights Reserved. + * + * Contributor(s): + * + * Alternatively, the contents of this file may be used under the terms of + * either of the GNU General Public License Version 2 or later (the "GPL"), + * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), + * in which case the provisions of the GPL or the LGPL are applicable instead + * of those above. If you wish to allow use of your version of this file only + * under the terms of either the GPL or the LGPL, and not to allow others to + * use your version of this file under the terms of the MPL, indicate your + * decision by deleting the provisions above and replace them with the notice + * and other provisions required by the GPL or the LGPL. If you do not delete + * the provisions above, a recipient may use your version of this file under + * the terms of any one of the MPL, the GPL or the LGPL. + * + * ***** END LICENSE BLOCK ***** */ + +/* @(#)s_cbrt.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +#include "fdlibm.h" + +/* cbrt(x) + * Return cube root of x + */ +#ifdef __STDC__ +static const unsigned +#else +static unsigned +#endif + B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ + B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ + +#ifdef __STDC__ +static const double +#else +static double +#endif +C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */ +D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */ +E = 1.41428571428571436819e+00, /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */ +F = 1.60714285714285720630e+00, /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */ +G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */ + +#ifdef __STDC__ + double fd_cbrt(double x) +#else + double fd_cbrt(x) + double x; +#endif +{ + fd_twoints u; + int hx; + double r,s,t=0.0,w; + unsigned sign; + + u.d = x; + hx = __HI(u); /* high word of x */ + sign=hx&0x80000000; /* sign= sign(x) */ + hx ^=sign; + if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */ + if((hx|__LO(u))==0) { + x = u.d; + return(x); /* cbrt(0) is itself */ + } + u.d = x; + __HI(u) = hx; /* x <- |x| */ + x = u.d; + /* rough cbrt to 5 bits */ + if(hx<0x00100000) /* subnormal number */ + {u.d = t; __HI(u)=0x43500000; t=u.d; /* set t= 2**54 */ + t*=x; __HI(u)=__HI(u)/3+B2; + } + else { + u.d = t; __HI(u)=hx/3+B1; t = u.d; + } + + + /* new cbrt to 23 bits, may be implemented in single precision */ + r=t*t/x; + s=C+r*t; + t*=G+F/(s+E+D/s); + + /* chopped to 20 bits and make it larger than cbrt(x) */ + u.d = t; + __LO(u)=0; __HI(u)+=0x00000001; + t = u.d; + + /* one step newton iteration to 53 bits with error less than 0.667 ulps */ + s=t*t; /* t*t is exact */ + r=x/s; + w=t+t; + r=(r-t)/(w+r); /* r-s is exact */ + t=t+t*r; + + /* retore the sign bit */ + u.d = t; + __HI(u) |= sign; + t = u.d; + return(t); +} |