1 #include <complex.h>
2 #include <math.h>
3 
4 // http://www.mathpropress.com/stan/bibliography/complexSquareRoot.pdf
csqrt(double complex z)5 double complex csqrt(double complex z)
6 {
7 	// 1 div sqrt(2)
8 	static double sqrt2 = 0.7071067811865475;
9 	double a = creal(z);
10 	double b = cimag(z);
11 	double sqrt_a_b = sqrt(a*a + b*b);
12 
13 	double realpart = sqrt(a + sqrt_a_b) * sqrt2;
14 	double imagpart = sqrt(sqrt_a_b - a) / sqrt2;
15 	if (b < 0)
16 		imagpart = -imagpart;
17 
18 	return realpart + imagpart * (double complex) I;
19 }
20