From Prof
. Kahan at UC at Berkeley
* Copyright (c) 1985 Regents of the University of California.
* Use and reproduction of this software are granted in accordance with
* the terms and conditions specified in the Berkeley Software License
* Agreement (in particular, this entails acknowledgement of the programs'
* source, and inclusion of this notice) with the additional understanding
* that 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 "sendbug 4bsd-bugs@BERKELEY", to the authors.
static char sccsid
[] = "@(#)cbrt.c 1.1 (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.
static 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-s is exact */
/* retore the sign bit */