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