* 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
->real
, z
->imag
)) == 0.)
r
->real
= sqrt(0.5 * (mag
+ z
->real
) );
r
->imag
= z
->imag
/ r
->real
/ 2;
r
->imag
= sqrt(0.5 * (mag
- z
->real
) );
*(unsigned long*)&(r
->imag
) ^= SIGN_BIT
;
r
->real
= z
->imag
/ r
->imag
/2;