1*05a0b428SJohn Marino /* $OpenBSD: s_csqrt.c,v 1.6 2013/07/03 04:46:36 espie Exp $ */
2*05a0b428SJohn Marino /*
3*05a0b428SJohn Marino * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
4*05a0b428SJohn Marino *
5*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this software for any
6*05a0b428SJohn Marino * purpose with or without fee is hereby granted, provided that the above
7*05a0b428SJohn Marino * copyright notice and this permission notice appear in all copies.
8*05a0b428SJohn Marino *
9*05a0b428SJohn Marino * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
10*05a0b428SJohn Marino * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
11*05a0b428SJohn Marino * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
12*05a0b428SJohn Marino * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
13*05a0b428SJohn Marino * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
14*05a0b428SJohn Marino * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
15*05a0b428SJohn Marino * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
16*05a0b428SJohn Marino */
17*05a0b428SJohn Marino
18*05a0b428SJohn Marino /* csqrt()
19*05a0b428SJohn Marino *
20*05a0b428SJohn Marino * Complex square root
21*05a0b428SJohn Marino *
22*05a0b428SJohn Marino *
23*05a0b428SJohn Marino *
24*05a0b428SJohn Marino * SYNOPSIS:
25*05a0b428SJohn Marino *
26*05a0b428SJohn Marino * double complex csqrt();
27*05a0b428SJohn Marino * double complex z, w;
28*05a0b428SJohn Marino *
29*05a0b428SJohn Marino * w = csqrt (z);
30*05a0b428SJohn Marino *
31*05a0b428SJohn Marino *
32*05a0b428SJohn Marino *
33*05a0b428SJohn Marino * DESCRIPTION:
34*05a0b428SJohn Marino *
35*05a0b428SJohn Marino *
36*05a0b428SJohn Marino * If z = x + iy, r = |z|, then
37*05a0b428SJohn Marino *
38*05a0b428SJohn Marino * 1/2
39*05a0b428SJohn Marino * Re w = [ (r + x)/2 ] ,
40*05a0b428SJohn Marino *
41*05a0b428SJohn Marino * 1/2
42*05a0b428SJohn Marino * Im w = [ (r - x)/2 ] .
43*05a0b428SJohn Marino *
44*05a0b428SJohn Marino * Cancellation error in r-x or r+x is avoided by using the
45*05a0b428SJohn Marino * identity 2 Re w Im w = y.
46*05a0b428SJohn Marino *
47*05a0b428SJohn Marino * Note that -w is also a square root of z. The root chosen
48*05a0b428SJohn Marino * is always in the right half plane and Im w has the same sign as y.
49*05a0b428SJohn Marino *
50*05a0b428SJohn Marino *
51*05a0b428SJohn Marino *
52*05a0b428SJohn Marino * ACCURACY:
53*05a0b428SJohn Marino *
54*05a0b428SJohn Marino * Relative error:
55*05a0b428SJohn Marino * arithmetic domain # trials peak rms
56*05a0b428SJohn Marino * DEC -10,+10 25000 3.2e-17 9.6e-18
57*05a0b428SJohn Marino * IEEE -10,+10 1,000,000 2.9e-16 6.1e-17
58*05a0b428SJohn Marino *
59*05a0b428SJohn Marino */
60*05a0b428SJohn Marino
61*05a0b428SJohn Marino #include <complex.h>
62*05a0b428SJohn Marino #include <float.h>
63*05a0b428SJohn Marino #include <math.h>
64*05a0b428SJohn Marino
65*05a0b428SJohn Marino double complex
csqrt(double complex z)66*05a0b428SJohn Marino csqrt(double complex z)
67*05a0b428SJohn Marino {
68*05a0b428SJohn Marino double complex w;
69*05a0b428SJohn Marino double x, y, r, t, scale;
70*05a0b428SJohn Marino
71*05a0b428SJohn Marino x = creal (z);
72*05a0b428SJohn Marino y = cimag (z);
73*05a0b428SJohn Marino
74*05a0b428SJohn Marino if (y == 0.0) {
75*05a0b428SJohn Marino if (x == 0.0) {
76*05a0b428SJohn Marino w = 0.0 + y * I;
77*05a0b428SJohn Marino }
78*05a0b428SJohn Marino else {
79*05a0b428SJohn Marino r = fabs (x);
80*05a0b428SJohn Marino r = sqrt (r);
81*05a0b428SJohn Marino if (x < 0.0) {
82*05a0b428SJohn Marino w = 0.0 + r * I;
83*05a0b428SJohn Marino }
84*05a0b428SJohn Marino else {
85*05a0b428SJohn Marino w = r + y * I;
86*05a0b428SJohn Marino }
87*05a0b428SJohn Marino }
88*05a0b428SJohn Marino return (w);
89*05a0b428SJohn Marino }
90*05a0b428SJohn Marino if (x == 0.0) {
91*05a0b428SJohn Marino r = fabs (y);
92*05a0b428SJohn Marino r = sqrt (0.5*r);
93*05a0b428SJohn Marino if (y > 0)
94*05a0b428SJohn Marino w = r + r * I;
95*05a0b428SJohn Marino else
96*05a0b428SJohn Marino w = r - r * I;
97*05a0b428SJohn Marino return (w);
98*05a0b428SJohn Marino }
99*05a0b428SJohn Marino /* Rescale to avoid internal overflow or underflow. */
100*05a0b428SJohn Marino if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) {
101*05a0b428SJohn Marino x *= 0.25;
102*05a0b428SJohn Marino y *= 0.25;
103*05a0b428SJohn Marino scale = 2.0;
104*05a0b428SJohn Marino }
105*05a0b428SJohn Marino else {
106*05a0b428SJohn Marino x *= 1.8014398509481984e16; /* 2^54 */
107*05a0b428SJohn Marino y *= 1.8014398509481984e16;
108*05a0b428SJohn Marino scale = 7.450580596923828125e-9; /* 2^-27 */
109*05a0b428SJohn Marino #if 0
110*05a0b428SJohn Marino x *= 4.0;
111*05a0b428SJohn Marino y *= 4.0;
112*05a0b428SJohn Marino scale = 0.5;
113*05a0b428SJohn Marino #endif
114*05a0b428SJohn Marino }
115*05a0b428SJohn Marino w = x + y * I;
116*05a0b428SJohn Marino r = cabs(w);
117*05a0b428SJohn Marino if (x > 0) {
118*05a0b428SJohn Marino t = sqrt(0.5 * r + 0.5 * x);
119*05a0b428SJohn Marino r = scale * fabs((0.5 * y) / t);
120*05a0b428SJohn Marino t *= scale;
121*05a0b428SJohn Marino }
122*05a0b428SJohn Marino else {
123*05a0b428SJohn Marino r = sqrt( 0.5 * r - 0.5 * x );
124*05a0b428SJohn Marino t = scale * fabs( (0.5 * y) / r );
125*05a0b428SJohn Marino r *= scale;
126*05a0b428SJohn Marino }
127*05a0b428SJohn Marino if (y < 0)
128*05a0b428SJohn Marino w = t - r * I;
129*05a0b428SJohn Marino else
130*05a0b428SJohn Marino w = t + r * I;
131*05a0b428SJohn Marino return (w);
132*05a0b428SJohn Marino }
133*05a0b428SJohn Marino
134*05a0b428SJohn Marino #if LDBL_MANT_DIG == DBL_MANT_DIG
135*05a0b428SJohn Marino __strong_alias(csqrtl, csqrt);
136*05a0b428SJohn Marino #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
137