diff options
Diffstat (limited to 'sysdeps/ieee754/cbrt.c')
| -rw-r--r-- | sysdeps/ieee754/cbrt.c | 120 |
1 files changed, 120 insertions, 0 deletions
diff --git a/sysdeps/ieee754/cbrt.c b/sysdeps/ieee754/cbrt.c new file mode 100644 index 0000000000..fe5fb95511 --- /dev/null +++ b/sysdeps/ieee754/cbrt.c @@ -0,0 +1,120 @@ +/* + * Copyright (c) 1985, 1993 + * The Regents of the University of California. All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. All advertising materials mentioning features or use of this software + * must display the following acknowledgement: + * This product includes software developed by the University of + * California, Berkeley and its contributors. + * 4. Neither the name of the University nor the names of its contributors + * may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#ifndef lint +static char sccsid[] = "@(#)cbrt.c 8.1 (Berkeley) 6/4/93"; +#endif /* not lint */ + +#include <sys/cdefs.h> + +/* kahan's cube root (53 bits IEEE double precision) + * for IEEE machines only + * coded in C by K.C. Ng, 4/30/85 + * + * Accuracy: + * better than 0.667 ulps according to an error analysis. Maximum + * error observed was 0.666 ulps in an 1,000,000 random arguments test. + * + * Warning: this code is semi machine dependent; the ordering of words in + * a floating point number must be known in advance. I assume that the + * long interger at the address of a floating point number will be the + * leading 32 bits of that floating point number (i.e., sign, exponent, + * and the 20 most significant bits). + * On a National machine, it has different ordering; therefore, this code + * must be compiled with flag -DNATIONAL. + */ +#if !defined(vax)&&!defined(tahoe) + +static const unsigned long + B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ + B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ +static const double + C= 19./35., + D= -864./1225., + E= 99./70., + F= 45./28., + G= 5./14.; + +double cbrt(x) +double x; +{ + double r,s,t=0.0,w; + unsigned long *px = (unsigned long *) &x, + *pt = (unsigned long *) &t, + mexp,sign; + +#ifdef national /* ordering of words in a floating points number */ + const int n0=1,n1=0; +#else /* national */ + const int n0=0,n1=1; +#endif /* national */ + + mexp=px[n0]&0x7ff00000; + if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */ + if(x==0.0) return(x); /* cbrt(0) is itself */ + + sign=px[n0]&0x80000000; /* sign= sign(x) */ + px[n0] ^= sign; /* x=|x| */ + + + /* rough cbrt to 5 bits */ + if(mexp==0) /* subnormal number */ + {pt[n0]=0x43500000; /* set t= 2**54 */ + t*=x; pt[n0]=pt[n0]/3+B2; + } + else + pt[n0]=px[n0]/3+B1; + + + /* 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) */ + pt[n1]=0; pt[n0]+=0x00000001; + + + /* 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-t is exact */ + t=t+t*r; + + + /* retore the sign bit */ + pt[n0] |= sign; + return(t); +} +#endif |