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