1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak __csqrtl = csqrtl
31
32 #include "libm.h" /* fabsl/isinfl/sqrtl */
33 #include "complex_wrapper.h"
34 #include "longdouble.h"
35
36 /* INDENT OFF */
37 static const long double
38 twom9001 = 2.6854002716003034957421765100615693043656e-2710L,
39 twom4500 = 2.3174987687592429423263242862381544149252e-1355L,
40 two8999 = 9.3095991180122343502582347372163290310934e+2708L,
41 two4500 = 4.3149968987270974283777803545571722250806e+1354L,
42 zero = 0.0L,
43 half = 0.5L,
44 two = 2.0L;
45 /* INDENT ON */
46
47 ldcomplex
csqrtl(ldcomplex z)48 csqrtl(ldcomplex z) {
49 ldcomplex ans;
50 long double x, y, t, ax, ay;
51 int n, ix, iy, hx, hy;
52
53 x = LD_RE(z);
54 y = LD_IM(z);
55 hx = HI_XWORD(x);
56 hy = HI_XWORD(y);
57 ix = hx & 0x7fffffff;
58 iy = hy & 0x7fffffff;
59 ay = fabsl(y);
60 ax = fabsl(x);
61 if (ix >= 0x7fff0000 || iy >= 0x7fff0000) {
62 /* x or y is Inf or NaN */
63 if (isinfl(y))
64 LD_IM(ans) = LD_RE(ans) = ay;
65 else if (isinfl(x)) {
66 if (hx > 0) {
67 LD_RE(ans) = ax;
68 LD_IM(ans) = ay * zero;
69 } else {
70 LD_RE(ans) = ay * zero;
71 LD_IM(ans) = ax;
72 }
73 } else
74 LD_IM(ans) = LD_RE(ans) = ax + ay;
75 } else if (y == zero) {
76 if (hx >= 0) {
77 LD_RE(ans) = sqrtl(ax);
78 LD_IM(ans) = zero;
79 } else {
80 LD_IM(ans) = sqrtl(ax);
81 LD_RE(ans) = zero;
82 }
83 } else if (ix >= iy) {
84 n = (ix - iy) >> 16;
85 #if defined(__x86) /* 64 significant bits */
86 if (n >= 35)
87 #else /* 113 significant bits */
88 if (n >= 60)
89 #endif
90 t = sqrtl(ax);
91 else if (ix >= 0x5f3f0000) { /* x > 2**8000 */
92 ax *= twom9001;
93 y *= twom9001;
94 t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
95 } else if (iy <= 0x20bf0000) { /* y < 2**-8000 */
96 ax *= two8999;
97 y *= two8999;
98 t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
99 } else
100 t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
101
102 if (hx >= 0) {
103 LD_RE(ans) = t;
104 LD_IM(ans) = ay / (t + t);
105 } else {
106 LD_IM(ans) = t;
107 LD_RE(ans) = ay / (t + t);
108 }
109 } else {
110 n = (iy - ix) >> 16;
111 #if defined(__x86) /* 64 significant bits */
112 if (n >= 35) { /* } */
113 #else /* 113 significant bits */
114 if (n >= 60) {
115 #endif
116 if (n >= 120)
117 t = sqrtl(half * ay);
118 else if (iy >= 0x7ffe0000)
119 t = sqrtl(half * ay + half * ax);
120 else if (ix <= 0x00010000)
121 t = half * (sqrtl(two * (ax + ay)));
122 else
123 t = sqrtl(half * (ax + ay));
124 } else if (iy >= 0x5f3f0000) { /* y > 2**8000 */
125 ax *= twom9001;
126 y *= twom9001;
127 t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
128 } else if (ix <= 0x20bf0000) {
129 ax *= two8999;
130 y *= two8999;
131 t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
132 } else
133 t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
134
135 if (hx >= 0) {
136 LD_RE(ans) = t;
137 LD_IM(ans) = ay / (t + t);
138 } else {
139 LD_IM(ans) = t;
140 LD_RE(ans) = ay / (t + t);
141 }
142 }
143 if (hy < 0)
144 LD_IM(ans) = -LD_IM(ans);
145 return (ans);
146 }
147