| 1 | #include "f2c.h" |
| 2 | |
| 3 | #ifdef KR_headers |
| 4 | extern double sqrt(), f__cabs(); |
| 5 | |
| 6 | VOID c_sqrt(r, z) complex *r, *z; |
| 7 | #else |
| 8 | #undef abs |
| 9 | #include "math.h" |
| 10 | extern double f__cabs(double, double); |
| 11 | |
| 12 | void c_sqrt(complex *r, complex *z) |
| 13 | #endif |
| 14 | { |
| 15 | double mag, t; |
| 16 | |
| 17 | if( (mag = f__cabs(z->r, z->i)) == 0.) |
| 18 | r->r = r->i = 0.; |
| 19 | else if(z->r > 0) |
| 20 | { |
| 21 | r->r = t = sqrt(0.5 * (mag + z->r) ); |
| 22 | t = z->i / t; |
| 23 | r->i = 0.5 * t; |
| 24 | } |
| 25 | else |
| 26 | { |
| 27 | t = sqrt(0.5 * (mag - z->r) ); |
| 28 | if(z->i < 0) |
| 29 | t = -t; |
| 30 | r->i = t; |
| 31 | t = z->i / t; |
| 32 | r->r = 0.5 * t; |
| 33 | } |
| 34 | } |