1 /* Copyright JS Foundation and other contributors, http://js.foundation
2 *
3 * Licensed under the Apache License, Version 2.0 (the "License");
4 * you may not use this file except in compliance with the License.
5 * You may obtain a copy of the License at
6 *
7 * http://www.apache.org/licenses/LICENSE-2.0
8 *
9 * Unless required by applicable law or agreed to in writing, software
10 * distributed under the License is distributed on an "AS IS" BASIS
11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 * See the License for the specific language governing permissions and
13 * limitations under the License.
14 *
15 * This file is based on work under the following copyright and permission
16 * notice:
17 *
18 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
19 *
20 * Developed at SunSoft, a Sun Microsystems, Inc. business.
21 * Permission to use, copy, modify, and distribute this
22 * software is freely granted, provided that this notice
23 * is preserved.
24 *
25 * @(#)s_cbrt.c 1.3 95/01/18
26 */
27
28 #include "jerry-libm-internal.h"
29
30 /* cbrt(x)
31 * Return cube root of x
32 */
33
34 #define B1 715094163 /* B1 = (682 - 0.03306235651) * 2**20 */
35 #define B2 696219795 /* B2 = (664 - 0.03306235651) * 2**20 */
36 #define C 5.42857142857142815906e-01 /* 19/35 = 0x3FE15F15, 0xF15F15F1 */
37 #define D -7.05306122448979611050e-01 /* -864/1225 = 0xBFE691DE, 0x2532C834 */
38 #define E 1.41428571428571436819e+00 /* 99/70 = 0x3FF6A0EA, 0x0EA0EA0F */
39 #define F 1.60714285714285720630e+00 /* 45/28 = 0x3FF9B6DB, 0x6DB6DB6E */
40 #define G 3.57142857142857150787e-01 /* 5/14 = 0x3FD6DB6D, 0xB6DB6DB7 */
41
42 double
cbrt(double x)43 cbrt (double x)
44 {
45 double r, s, w;
46 double_accessor t, temp;
47 unsigned int sign;
48 t.dbl = 0.0;
49 temp.dbl = x;
50
51 sign = temp.as_int.hi & 0x80000000; /* sign = sign(x) */
52 temp.as_int.hi ^= sign;
53
54 if (temp.as_int.hi >= 0x7ff00000)
55 {
56 /* cbrt(NaN, INF) is itself */
57 return (x + x);
58 }
59 if ((temp.as_int.hi | temp.as_int.lo) == 0)
60 {
61 /* cbrt(0) is itself */
62 return (x);
63 }
64 /* rough cbrt to 5 bits */
65 if (temp.as_int.hi < 0x00100000) /* subnormal number */
66 {
67 t.as_int.hi = 0x43500000; /* set t= 2**54 */
68 t.dbl *= temp.dbl;
69 t.as_int.hi = t.as_int.hi / 3 + B2;
70 }
71 else
72 {
73 t.as_int.hi = temp.as_int.hi / 3 + B1;
74 }
75
76 /* new cbrt to 23 bits, may be implemented in single precision */
77 r = t.dbl * t.dbl / temp.dbl;
78 s = C + r * t.dbl;
79 t.dbl *= G + F / (s + E + D / s);
80
81 /* chopped to 20 bits and make it larger than cbrt(x) */
82 t.as_int.lo = 0;
83 t.as_int.hi += 0x00000001;
84
85 /* one step newton iteration to 53 bits with error less than 0.667 ulps */
86 s = t.dbl * t.dbl; /* t*t is exact */
87 r = temp.dbl / s;
88 w = t.dbl + t.dbl;
89 r = (r - t.dbl) / (w + r); /* r-s is exact */
90 t.dbl = t.dbl + (t.dbl * r);
91
92 /* retore the sign bit */
93 t.as_int.hi |= sign;
94 return (t.dbl);
95 } /* cbrt */
96
97 #undef B1
98 #undef B2
99 #undef C
100 #undef D
101 #undef E
102 #undef F
103 #undef G
104