1#include "f2c.h"
2
3#ifdef KR_headers
4extern double sqrt(), f__cabs();
5
6VOID c_sqrt(r, z) complex *r, *z;
7#else
8#undef abs
9#include "math.h"
10extern double f__cabs(double, double);
11
12void c_sqrt(complex *r, complex *z)
13#endif
14{
15	double mag, t;
16	double zi = z->i, zr = z->r;
17
18	if( (mag = f__cabs(zr, zi)) == 0.)
19		r->r = r->i = 0.;
20	else if(zr > 0)
21		{
22		r->r = t = sqrt(0.5 * (mag + zr) );
23		t = zi / t;
24		r->i = 0.5 * t;
25		}
26	else
27		{
28		t = sqrt(0.5 * (mag - zr) );
29		if(zi < 0)
30			t = -t;
31		r->i = t;
32		t = zi / t;
33		r->r = 0.5 * t;
34		}
35	}
36