|
|
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.