* Copyright (c) 1985 Regents of the University of California.
* Redistribution and use in source and binary forms are permitted
* provided that the above copyright notice and this paragraph are
* duplicated in all such forms and that any documentation,
* advertising materials, and other materials related to such
* distribution and use acknowledge that the software was developed
* by the University of California, Berkeley. The name of the
* University may not be used to endorse or promote products derived
* from this software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
* All recipients should regard themselves as participants in an ongoing
* research project and hence should feel obligated to report their
* experiences (good or bad) with these elementary function codes, using
* the sendbug(8) program, to the authors.
static char sccsid
[] = "@(#)cbrt.c 5.5 (Berkeley) %G%";
/* kahan's cube root (53 bits IEEE double precision)
* coded in C by K.C. Ng, 4/30/85
* 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 */
unsigned long *px
= (unsigned long *) &x
,
*pt
= (unsigned long *) &t
,
#ifdef national /* ordering of words in a floating points number */
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
;
/* new cbrt to 23 bits, may be implemented in single precision */
/* 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
=(r
-t
)/(w
+r
); /* r-t is exact */
/* retore the sign bit */