1 /* 2 * "@(#)c_sqrt.c 1.1" 3 */ 4 5 #include "complex" 6 7 c_sqrt(r, z) 8 complex *r, *z; 9 { 10 double mag, sqrt(), cabs(); 11 12 if( (mag = cabs(z->real, z->imag)) == 0.) 13 r->real = r->imag = 0.; 14 else if(z->real > 0) 15 { 16 r->real = sqrt(0.5 * (mag + z->real) ); 17 r->imag = z->imag / r->real / 2; 18 } 19 else 20 { 21 r->imag = sqrt(0.5 * (mag - z->real) ); 22 if(z->imag < 0) 23 r->imag = - r->imag; 24 r->real = z->imag / r->imag /2; 25 } 26 } 27