system/freebsd-src
system/freebsd-src
1
0
Fork 0
mirror of https://github.com/freebsd/freebsd-src synced 2024-10-05 08:00:30 +00:00
Code Issues Packages Projects Releases Wiki Activity

Added comments about the magic behind

Browse source 
<cbrt(x) in bits> ~= <x in bits>/3 + BIAS.
Keep the large comments only in the double version as usual.

Fixed some style bugs (mainly grammar and spelling errors in comments).
This commit is contained in:
Bruce Evans 2005-12-11 19:51:30 +00:00
parent 288a8c86cb
commit af7f99131d
Notes: svn2git 2020-12-20 02:59:44 +00:00 
svn path=/head/; revision=153306
2 changed files with 25 additions and 14 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

32
lib/msun/src/s_cbrt.c
View file

@ -21,8 +21,8 @@ static char rcsid[] = "$FreeBSD$";
* Return cube root of x
*/
static const u_int32_t
B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
static const double
C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */
@ -48,7 +48,21 @@ cbrt(double x)
return(x); /* cbrt(0) is itself */
SET_HIGH_WORD(x,hx); /* x <- |x| */
/* rough cbrt to 5 bits */
/*
* Rough cbrt to 5 bits:
* cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
* where e is integral and >= 0, m is real and in [0, 1), and "/" and
* "%" are integer division and modulus with rounding towards minus
* infinity. The RHS is always >= the LHS and has a maximum relative
* error of about 1 in 16. Adding a bias of -0.03306235651 to the
* (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
* floating point representation, for finite positive normal values,
* ordinary integer divison of the value in bits magically gives
* almost exactly the RHS of the above provided we first subtract the
* exponent bias (1023 for doubles) and later add it back. We do the
* subtraction virtually to keep e >= 0 so that ordinary integer
* division rounds towards minus infinity; this is also efficient.
*/
if(hx<0x00100000) /* subnormal number */
{SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */
t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2);
@ -56,25 +70,23 @@ cbrt(double x)
else
SET_HIGH_WORD(t,hx/3+B1);
/* new cbrt to 23 bits, may be implemented in single precision */
/* 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) */
/* chop t to 20 bits and make it larger than cbrt(x) */
GET_HIGH_WORD(high,t);
INSERT_WORDS(t,high+0x00000001,0);
/* one step newton iteration to 53 bits with error less than 0.667 ulps */
/* 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 */
r=(r-t)/(w+r); /* r-t is exact */
t=t+t*r;
/* retore the sign bit */
/* restore the sign bit */
GET_HIGH_WORD(high,t);
SET_HIGH_WORD(t,high|sign);
return(t);

7
lib/msun/src/s_cbrtf.c
View file

@ -25,8 +25,8 @@ static char rcsid[] = "$FreeBSD$";
* Return cube root of x
*/
static const unsigned
B1 = 709958130, /* B1 = (84+2/3-0.03306235651)*2**23 */
B2 = 642849266; /* B2 = (76+2/3-0.03306235651)*2**23 */
B1 = 709958130, /* B1 = (127-127.0/3-0.03306235651)*2**23 */
B2 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
static const float
C = 5.4285717010e-01, /* 19/35 = 0x3f0af8b0 */
@ -59,7 +59,6 @@ cbrtf(float x)
else
SET_FLOAT_WORD(t,hx/3+B1);
/* new cbrt to 23 bits */
r=t*t/x;
s=C+r*t;
@ -76,7 +75,7 @@ cbrtf(float x)
r=(r-t)/(w+r); /* r-t is exact */
t=t+t*r;
/* retore the sign bit */
/* restore the sign bit */
GET_FLOAT_WORD(high,t);
SET_FLOAT_WORD(t,high|sign);
return(t);