* Copyright (c) 1980 Regents of the University of California.
* All rights reserved. The Berkeley software License Agreement
* specifies the terms and conditions for redistribution.
#include <tahoemath/FP.h>
double mag
, sqrt(), cabs();
if( (mag
= cabs(z
->dreal
, z
->dimag
)) == 0.)
r
->dreal
= r
->dimag
= 0.;
r
->dreal
= sqrt(0.5 * (mag
+ z
->dreal
) );
r
->dimag
= z
->dimag
/ r
->dreal
/ 2;
r
->dimag
= sqrt(0.5 * (mag
- z
->dreal
) );
*((long int *)&r
->dimag
) ^= SIGN_BIT
;
r
->dreal
= z
->dimag
/ r
->dimag
/ 2;