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 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __csqrtl = csqrtl
32
33 #include "libm.h" /* fabsl/isinfl/sqrtl */
34 #include "complex_wrapper.h"
35 #include "longdouble.h"
36
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
46
47 ldcomplex
48 csqrtl(ldcomplex z)
49 {
50 ldcomplex ans;
51 long double x, y, t, ax, ay;
52 int n, ix, iy, hx, hy;
53
54 x = LD_RE(z);
55 y = LD_IM(z);
56 hx = HI_XWORD(x);
57 hy = HI_XWORD(y);
58 ix = hx & 0x7fffffff;
59 iy = hy & 0x7fffffff;
60 ay = fabsl(y);
61 ax = fabsl(x);
62
63 if (ix >= 0x7fff0000 || iy >= 0x7fff0000) {
64 /* x or y is Inf or NaN */
65 if (isinfl(y)) {
66 LD_IM(ans) = LD_RE(ans) = ay;
67 } else if (isinfl(x)) {
68 if (hx > 0) {
69 LD_RE(ans) = ax;
70 LD_IM(ans) = ay * zero;
71 } else {
72 LD_RE(ans) = ay * zero;
73 LD_IM(ans) = ax;
74 }
75 } else {
76 LD_IM(ans) = LD_RE(ans) = ax + ay;
77 }
78 } else if (y == zero) {
79 if (hx >= 0) {
80 LD_RE(ans) = sqrtl(ax);
81 LD_IM(ans) = zero;
82 } else {
83 LD_IM(ans) = sqrtl(ax);
84 LD_RE(ans) = zero;
85 }
86 } else if (ix >= iy) {
87 n = (ix - iy) >> 16;
88 #if defined(__x86) /* 64 significant bits */
89 if (n >= 35)
90 #else /* 113 significant bits */
91 if (n >= 60)
92 #endif
93 {
94 t = sqrtl(ax);
95 } else if (ix >= 0x5f3f0000) { /* x > 2**8000 */
96 ax *= twom9001;
97 y *= twom9001;
98 t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
99 } else if (iy <= 0x20bf0000) { /* y < 2**-8000 */
100 ax *= two8999;
101 y *= two8999;
102 t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
103 } else {
104 t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
105 }
106
107 if (hx >= 0) {
108 LD_RE(ans) = t;
109 LD_IM(ans) = ay / (t + t);
110 } else {
111 LD_IM(ans) = t;
112 LD_RE(ans) = ay / (t + t);
113 }
114 } else {
115 n = (iy - ix) >> 16;
116 #if defined(__x86) /* 64 significant bits */
117 if (n >= 35) { /* } */
118 #else /* 113 significant bits */
119 if (n >= 60) {
120 #endif
121
122 if (n >= 120)
123 t = sqrtl(half * ay);
124 else if (iy >= 0x7ffe0000)
125 t = sqrtl(half * ay + half * ax);
126 else if (ix <= 0x00010000)
127 t = half * (sqrtl(two * (ax + ay)));
128 else
129 t = sqrtl(half * (ax + ay));
130 } else if (iy >= 0x5f3f0000) { /* y > 2**8000 */
131 ax *= twom9001;
132 y *= twom9001;
133 t = two4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
134 } else if (ix <= 0x20bf0000) {
135 ax *= two8999;
136 y *= two8999;
137 t = twom4500 * sqrtl(ax + sqrtl(ax * ax + y * y));
138 } else {
139 t = sqrtl(half * (ax + sqrtl(ax * ax + y * y)));
140 }
141
142 if (hx >= 0) {
143 LD_RE(ans) = t;
144 LD_IM(ans) = ay / (t + t);
145 } else {
146 LD_IM(ans) = t;
147 LD_RE(ans) = ay / (t + t);
148 }
149 }
150
151 if (hy < 0)
152 LD_IM(ans) = -LD_IM(ans);
153
154 return (ans);
155 }