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