|
|
1.1 root 1: /*
2: * "@(#)z_sqrt.c 1.1"
3: */
4:
5: #include "complex"
6:
7: /* THIS IS BASED ON THE TAHOE FP REPRESENTATION */
8: #include <FP.h>
9:
10: z_sqrt(r, z)
11: dcomplex *r, *z;
12: {
13: double mag, sqrt(), zabs();
14:
15: if( (mag = zabs(z->dreal, z->dimag)) == 0.)
16: r->dreal = r->dimag = 0.;
17: else if(z->dreal > 0)
18: {
19: r->dreal = sqrt(0.5 * (mag + z->dreal) );
20: r->dimag = z->dimag / r->dreal / 2;
21: }
22: else
23: {
24: r->dimag = sqrt(0.5 * (mag - z->dreal) );
25: if(z->dimag < 0)
26: *((long int *)&r->dimag) ^= SIGN_BIT;
27: r->dreal = z->dimag / r->dimag / 2;
28: }
29: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.