#include "complex" c_sqrt(r, z) complex *r, *z; { double mag, sqrt(), cabs(); if( (mag = cabs(z->real, z->imag)) == 0.) r->real = r->imag = 0.; else if(z->real > 0) { r->real = sqrt(0.5 * (mag + z->real) ); r->imag = z->imag / r->real / 2; } else { r->imag = sqrt(0.5 * (mag - z->real) ); if(z->imag < 0) r->imag = - r->imag; r->real = z->imag / r->imag /2; } }