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