s_cbrt.c revision 153520
12116Sjkh/* @(#)s_cbrt.c 5.1 93/09/24 */
22116Sjkh/*
32116Sjkh * ====================================================
42116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
52116Sjkh *
62116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business.
72116Sjkh * Permission to use, copy, modify, and distribute this
88870Srgrimes * software is freely granted, provided that this notice
92116Sjkh * is preserved.
102116Sjkh * ====================================================
11153386Sbde *
12153386Sbde * Optimized by Bruce D. Evans.
132116Sjkh */
142116Sjkh
152116Sjkh#ifndef lint
1650476Speterstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/s_cbrt.c 153520 2005-12-19 00:22:03Z bde $";
172116Sjkh#endif
182116Sjkh
192116Sjkh#include "math.h"
202116Sjkh#include "math_private.h"
212116Sjkh
222116Sjkh/* cbrt(x)
232116Sjkh * Return cube root of x
242116Sjkh */
252116Sjkhstatic const u_int32_t
26153306Sbde	B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
27153306Sbde	B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
282116Sjkh
29153520Sbde/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */
302116Sjkhstatic const double
31153520SbdeP0 =  1.87595182427177009643,		/* 0x3ffe03e6, 0x0f61e692 */
32153520SbdeP1 = -1.88497979543377169875,		/* 0xbffe28e0, 0x92f02420 */
33153520SbdeP2 =  1.621429720105354466140,		/* 0x3ff9f160, 0x4a49d6c2 */
34153520SbdeP3 = -0.758397934778766047437,		/* 0xbfe844cb, 0xbee751d9 */
35153520SbdeP4 =  0.145996192886612446982;		/* 0x3fc2b000, 0xd4e4edd7 */
362116Sjkh
3797413Salfreddouble
3897413Salfredcbrt(double x)
392116Sjkh{
402116Sjkh	int32_t	hx;
41153517Sbde	union {
42153517Sbde	    double value;
43153517Sbde	    uint64_t bits;
44153517Sbde	} u;
452116Sjkh	double r,s,t=0.0,w;
46153517Sbde	uint64_t bits;
472116Sjkh	u_int32_t sign;
482116Sjkh	u_int32_t high,low;
492116Sjkh
502116Sjkh	GET_HIGH_WORD(hx,x);
512116Sjkh	sign=hx&0x80000000; 		/* sign= sign(x) */
522116Sjkh	hx  ^=sign;
532116Sjkh	if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
542116Sjkh	GET_LOW_WORD(low,x);
558870Srgrimes	if((hx|low)==0)
56153517Sbde	    return(x);			/* cbrt(0) is itself */
572116Sjkh
58153306Sbde    /*
59153306Sbde     * Rough cbrt to 5 bits:
60153306Sbde     *    cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
61153306Sbde     * where e is integral and >= 0, m is real and in [0, 1), and "/" and
62153306Sbde     * "%" are integer division and modulus with rounding towards minus
63153306Sbde     * infinity.  The RHS is always >= the LHS and has a maximum relative
64153306Sbde     * error of about 1 in 16.  Adding a bias of -0.03306235651 to the
65153306Sbde     * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
66153306Sbde     * floating point representation, for finite positive normal values,
67153306Sbde     * ordinary integer divison of the value in bits magically gives
68153306Sbde     * almost exactly the RHS of the above provided we first subtract the
69153306Sbde     * exponent bias (1023 for doubles) and later add it back.  We do the
70153306Sbde     * subtraction virtually to keep e >= 0 so that ordinary integer
71153306Sbde     * division rounds towards minus infinity; this is also efficient.
72153306Sbde     */
73153382Sbde	if(hx<0x00100000) { 		/* subnormal number */
74153382Sbde	    SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */
75153382Sbde	    t*=x;
76153382Sbde	    GET_HIGH_WORD(high,t);
77153386Sbde	    SET_HIGH_WORD(t,sign|((high&0x7fffffff)/3+B2));
78153382Sbde	} else
79153386Sbde	    SET_HIGH_WORD(t,sign|(hx/3+B1));
802116Sjkh
81153447Sbde    /*
82153520Sbde     * New cbrt to 23 bits:
83153520Sbde     *    cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x)
84153520Sbde     * where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r)
85153520Sbde     * to within 2**-23.5 when |r - 1| < 1/10.  The rough approximation
86153520Sbde     * has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this
87153520Sbde     * gives us bounds for r = t**3/x.
88153520Sbde     *
89153520Sbde     * Try to optimize for parallel evaluation as in k_tanf.c.
90153447Sbde     */
91153520Sbde	r=(t*t)*(t/x);
92153520Sbde	t=t*((P0+r*(P1+r*P2))+((r*r)*r)*(P3+r*P4));
932116Sjkh
94153517Sbde    /*
95153520Sbde     * Round t away from zero to 23 bits (sloppily except for ensuring that
96153517Sbde     * the result is larger in magnitude than cbrt(x) but not much more than
97153520Sbde     * 2 23-bit ulps larger).  With rounding towards zero, the error bound
98153520Sbde     * would be ~5/6 instead of ~4/6.  With a maximum error of 2 23-bit ulps
99153517Sbde     * in the rounded t, the infinite-precision error in the Newton
100153517Sbde     * approximation barely affects third digit in the the final error
101153520Sbde     * 0.667; the error in the rounded t can be up to about 3 23-bit ulps
102153517Sbde     * before the final error is larger than 0.667 ulps.
103153517Sbde     */
104153517Sbde	u.value=t;
105153520Sbde	u.bits=(u.bits+0x80000000)&0xffffffffc0000000ULL;
106153517Sbde	t=u.value;
1072116Sjkh
108153517Sbde    /* one step Newton iteration to 53 bits with error < 0.667 ulps */
109153517Sbde	s=t*t;				/* t*t is exact */
110153517Sbde	r=x/s;				/* error <= 0.5 ulps; |r| < |t| */
111153517Sbde	w=t+t;				/* t+t is exact */
112153517Sbde	r=(r-t)/(w+r);			/* r-t is exact; w+r ~= 3*t */
113153517Sbde	t=t+t*r;			/* error <= 0.5 + 0.5/3 + epsilon */
1142116Sjkh
1152116Sjkh	return(t);
1162116Sjkh}
117