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