| 667 | |
| 668 | extern double f__cabs(double, double); |
| 669 | void z_sqrt(doublecomplex *r, doublecomplex *z) |
| 670 | #endif |
| 671 | { |
| 672 | double mag; |
| 673 | |
| 674 | if( (mag = f__cabs(z->r, z->i)) == 0.) |
| 675 | r->r = r->i = 0.; |
| 676 | else if(z->r > 0) |
| 677 | { |
| 678 | r->r = sqrt(0.5 * (mag + z->r) ); |
| 679 | r->i = z->i / r->r / 2; |
| 680 | } |
| 681 | else |
| 682 | { |
| 683 | r->i = sqrt(0.5 * (mag - z->r) ); |
| 684 | if(z->i < 0) |
| 685 | r->i = - r->i; |
| 686 | r->r = z->i / r->i / 2; |
| 687 | } |
| 688 | } |
| 689 | #ifdef __cplusplus |
| 690 | extern "C" { |
| 691 | #endif |