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 __csqrt = csqrt 31 32 /* INDENT OFF */ 33 /* 34 * dcomplex csqrt(dcomplex z); 35 * 36 * 2 2 2 37 * Let w=r+i*s = sqrt(x+iy). Then (r + i s) = r - s + i 2sr = x + i y. 38 * 39 * Hence x = r*r-s*s, y = 2sr. 40 * 41 * Note that x*x+y*y = (s*s+r*r)**2. Thus, we have 42 * ________ 43 * 2 2 / 2 2 44 * (1) r + s = \/ x + y , 45 * 46 * 2 2 47 * (2) r - s = x 48 * 49 * (3) 2sr = y. 50 * 51 * Perform (1)-(2) and (1)+(2), we obtain 52 * 53 * 2 54 * (4) 2 r = hypot(x,y)+x, 55 * 56 * 2 57 * (5) 2*s = hypot(x,y)-x 58 * ________ 59 * / 2 2 60 * where hypot(x,y) = \/ x + y . 61 * 62 * In order to avoid numerical cancellation, we use formula (4) for 63 * positive x, and (5) for negative x. The other component is then 64 * computed by formula (3). 65 * 66 * 67 * ALGORITHM 68 * ------------------ 69 * 70 * (assume x and y are of medium size, i.e., no over/underflow in squaring) 71 * 72 * If x >=0 then 73 * ________ 74 * / 2 2 75 * 2 \/ x + y + x y 76 * r = ---------------------, s = -------; (6) 77 * 2 2 r 78 * 79 * (note that we choose sign(s) = sign(y) to force r >=0). 80 * Otherwise, 81 * ________ 82 * / 2 2 83 * 2 \/ x + y - x y 84 * s = ---------------------, r = -------; (7) 85 * 2 2 s 86 * 87 * EXCEPTION: 88 * 89 * One may use the polar coordinate of a complex number to justify the 90 * following exception cases: 91 * 92 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)): 93 * csqrt(+-0+ i 0 ) = 0 + i 0 94 * csqrt( x + i inf ) = inf + i inf for all x (including NaN) 95 * csqrt( x + i NaN ) = NaN + i NaN with invalid for finite x 96 * csqrt(-inf+ iy ) = 0 + i inf for finite positive-signed y 97 * csqrt(+inf+ iy ) = inf + i 0 for finite positive-signed y 98 * csqrt(-inf+ i NaN) = NaN +-i inf 99 * csqrt(+inf+ i NaN) = inf + i NaN 100 * csqrt(NaN + i y ) = NaN + i NaN for finite y 101 * csqrt(NaN + i NaN) = NaN + i NaN 102 */ 103 /* INDENT ON */ 104 105 #include "libm.h" /* fabs/sqrt */ 106 #include "complex_wrapper.h" 107 108 /* INDENT OFF */ 109 static const double 110 two300 = 2.03703597633448608627e+90, 111 twom300 = 4.90909346529772655310e-91, 112 two599 = 2.07475778444049647926e+180, 113 twom601 = 1.20495993255144205887e-181, 114 two = 2.0, 115 zero = 0.0, 116 half = 0.5; 117 /* INDENT ON */ 118 119 dcomplex 120 csqrt(dcomplex z) { 121 dcomplex ans; 122 double x, y, t, ax, ay; 123 int n, ix, iy, hx, hy, lx, ly; 124 125 x = D_RE(z); 126 y = D_IM(z); 127 hx = HI_WORD(x); 128 lx = LO_WORD(x); 129 hy = HI_WORD(y); 130 ly = LO_WORD(y); 131 ix = hx & 0x7fffffff; 132 iy = hy & 0x7fffffff; 133 ay = fabs(y); 134 ax = fabs(x); 135 if (ix >= 0x7ff00000 || iy >= 0x7ff00000) { 136 /* x or y is Inf or NaN */ 137 if (ISINF(iy, ly)) 138 D_IM(ans) = D_RE(ans) = ay; 139 else if (ISINF(ix, lx)) { 140 if (hx > 0) { 141 D_RE(ans) = ax; 142 D_IM(ans) = ay * zero; 143 } else { 144 D_RE(ans) = ay * zero; 145 D_IM(ans) = ax; 146 } 147 } else 148 D_IM(ans) = D_RE(ans) = ax + ay; 149 } else if ((iy | ly) == 0) { /* y = 0 */ 150 if (hx >= 0) { 151 D_RE(ans) = sqrt(ax); 152 D_IM(ans) = zero; 153 } else { 154 D_IM(ans) = sqrt(ax); 155 D_RE(ans) = zero; 156 } 157 } else if (ix >= iy) { 158 n = (ix - iy) >> 20; 159 if (n >= 30) { /* x >> y or y=0 */ 160 t = sqrt(ax); 161 } else if (ix >= 0x5f300000) { /* x > 2**500 */ 162 ax *= twom601; 163 y *= twom601; 164 t = two300 * sqrt(ax + sqrt(ax * ax + y * y)); 165 } else if (iy < 0x20b00000) { /* y < 2**-500 */ 166 ax *= two599; 167 y *= two599; 168 t = twom300 * sqrt(ax + sqrt(ax * ax + y * y)); 169 } else 170 t = sqrt(half * (ax + sqrt(ax * ax + ay * ay))); 171 if (hx >= 0) { 172 D_RE(ans) = t; 173 D_IM(ans) = ay / (t + t); 174 } else { 175 D_IM(ans) = t; 176 D_RE(ans) = ay / (t + t); 177 } 178 } else { 179 n = (iy - ix) >> 20; 180 if (n >= 30) { /* y >> x */ 181 if (n >= 60) 182 t = sqrt(half * ay); 183 else if (iy >= 0x7fe00000) 184 t = sqrt(half * ay + half * ax); 185 else if (ix <= 0x00100000) 186 t = half * sqrt(two * (ay + ax)); 187 else 188 t = sqrt(half * (ay + ax)); 189 } else if (iy >= 0x5f300000) { /* y > 2**500 */ 190 ax *= twom601; 191 y *= twom601; 192 t = two300 * sqrt(ax + sqrt(ax * ax + y * y)); 193 } else if (ix < 0x20b00000) { /* x < 2**-500 */ 194 ax *= two599; 195 y *= two599; 196 t = twom300 * sqrt(ax + sqrt(ax * ax + y * y)); 197 } else 198 t = sqrt(half * (ax + sqrt(ax * ax + ay * ay))); 199 if (hx >= 0) { 200 D_RE(ans) = t; 201 D_IM(ans) = ay / (t + t); 202 } else { 203 D_IM(ans) = t; 204 D_RE(ans) = ay / (t + t); 205 } 206 } 207 if (hy < 0) 208 D_IM(ans) = -D_IM(ans); 209 return (ans); 210 } 211