1*05a0b428SJohn Marino /* @(#)s_cbrt.c 5.1 93/09/24 */
2*05a0b428SJohn Marino /*
3*05a0b428SJohn Marino * ====================================================
4*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5*05a0b428SJohn Marino *
6*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
7*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
8*05a0b428SJohn Marino * software is freely granted, provided that this notice
9*05a0b428SJohn Marino * is preserved.
10*05a0b428SJohn Marino * ====================================================
11*05a0b428SJohn Marino */
12*05a0b428SJohn Marino
13*05a0b428SJohn Marino #include <float.h>
14*05a0b428SJohn Marino #include <math.h>
15*05a0b428SJohn Marino
16*05a0b428SJohn Marino #include "math_private.h"
17*05a0b428SJohn Marino
18*05a0b428SJohn Marino /* cbrt(x)
19*05a0b428SJohn Marino * Return cube root of x
20*05a0b428SJohn Marino */
21*05a0b428SJohn Marino static const u_int32_t
22*05a0b428SJohn Marino B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */
23*05a0b428SJohn Marino B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */
24*05a0b428SJohn Marino
25*05a0b428SJohn Marino static const double
26*05a0b428SJohn Marino C = 5.42857142857142815906e-01, /* 19/35 = 0x3FE15F15, 0xF15F15F1 */
27*05a0b428SJohn Marino D = -7.05306122448979611050e-01, /* -864/1225 = 0xBFE691DE, 0x2532C834 */
28*05a0b428SJohn Marino E = 1.41428571428571436819e+00, /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */
29*05a0b428SJohn Marino F = 1.60714285714285720630e+00, /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */
30*05a0b428SJohn Marino G = 3.57142857142857150787e-01; /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
31*05a0b428SJohn Marino
32*05a0b428SJohn Marino double
cbrt(double x)33*05a0b428SJohn Marino cbrt(double x)
34*05a0b428SJohn Marino {
35*05a0b428SJohn Marino int32_t hx;
36*05a0b428SJohn Marino double r,s,t=0.0,w;
37*05a0b428SJohn Marino u_int32_t sign;
38*05a0b428SJohn Marino u_int32_t high,low;
39*05a0b428SJohn Marino
40*05a0b428SJohn Marino GET_HIGH_WORD(hx,x);
41*05a0b428SJohn Marino sign=hx&0x80000000; /* sign= sign(x) */
42*05a0b428SJohn Marino hx ^=sign;
43*05a0b428SJohn Marino if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
44*05a0b428SJohn Marino GET_LOW_WORD(low,x);
45*05a0b428SJohn Marino if((hx|low)==0)
46*05a0b428SJohn Marino return(x); /* cbrt(0) is itself */
47*05a0b428SJohn Marino
48*05a0b428SJohn Marino SET_HIGH_WORD(x,hx); /* x <- |x| */
49*05a0b428SJohn Marino /* rough cbrt to 5 bits */
50*05a0b428SJohn Marino if(hx<0x00100000) /* subnormal number */
51*05a0b428SJohn Marino {SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */
52*05a0b428SJohn Marino t*=x; GET_HIGH_WORD(high,t); SET_HIGH_WORD(t,high/3+B2);
53*05a0b428SJohn Marino }
54*05a0b428SJohn Marino else
55*05a0b428SJohn Marino SET_HIGH_WORD(t,hx/3+B1);
56*05a0b428SJohn Marino
57*05a0b428SJohn Marino
58*05a0b428SJohn Marino /* new cbrt to 23 bits, may be implemented in single precision */
59*05a0b428SJohn Marino r=t*t/x;
60*05a0b428SJohn Marino s=C+r*t;
61*05a0b428SJohn Marino t*=G+F/(s+E+D/s);
62*05a0b428SJohn Marino
63*05a0b428SJohn Marino /* chopped to 20 bits and make it larger than cbrt(x) */
64*05a0b428SJohn Marino GET_HIGH_WORD(high,t);
65*05a0b428SJohn Marino INSERT_WORDS(t,high+0x00000001,0);
66*05a0b428SJohn Marino
67*05a0b428SJohn Marino
68*05a0b428SJohn Marino /* one step newton iteration to 53 bits with error less than 0.667 ulps */
69*05a0b428SJohn Marino s=t*t; /* t*t is exact */
70*05a0b428SJohn Marino r=x/s;
71*05a0b428SJohn Marino w=t+t;
72*05a0b428SJohn Marino r=(r-t)/(w+r); /* r-s is exact */
73*05a0b428SJohn Marino t=t+t*r;
74*05a0b428SJohn Marino
75*05a0b428SJohn Marino /* retore the sign bit */
76*05a0b428SJohn Marino GET_HIGH_WORD(high,t);
77*05a0b428SJohn Marino SET_HIGH_WORD(t,high|sign);
78*05a0b428SJohn Marino return(t);
79*05a0b428SJohn Marino }
80*05a0b428SJohn Marino
81*05a0b428SJohn Marino #if LDBL_MANT_DIG == DBL_MANT_DIG
82*05a0b428SJohn Marino __strong_alias(cbrtl, cbrt);
83*05a0b428SJohn Marino #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
84