home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / unix / volume10 / complex-lib / cxsqrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1987-07-14  |  1.3 KB  |  60 lines
  1. /*
  2.     CxSqrt -- compute square root of complex number
  3.  
  4.     CxSqrt( &c )    replaces  c  by  sqrt(c)  and returns  &c
  5.  
  6.     Note:    This is a double-valued function; the result of
  7.         CxSqrt() always has nonnegative imaginary part.
  8.  
  9.     inspired by Jeff Hanes' version
  10.  
  11.     last edit:    86/01/04    D A Gwyn
  12.  
  13.     SCCS ID:    @(#)cxsqrt.c    1.1 (modified for public version)
  14. */
  15.  
  16. #include    <math.h>
  17.  
  18. #include    <complex.h>
  19.  
  20. #define    Sgn( x )    ((x) == 0 ? 0 : (x) > 0 ? 1 : -1)
  21.  
  22. complex *
  23. CxSqrt( cp )
  24.     register complex    *cp;
  25.     {
  26.     /* record signs of original real & imaginary parts */
  27.     int            re_sign = Sgn( cp->re );
  28.     int            im_sign = Sgn( cp->im );
  29.  
  30.     /* special cases are not necessary; they are here for speed */
  31.  
  32.     if ( re_sign == 0 )
  33.         if ( im_sign == 0 )
  34.             ;        /* (0,0) already there */
  35.         else if ( im_sign > 0 )
  36.             cp->re = cp->im = sqrt( cp->im / 2.0 );
  37.         else            /* im_sign < 0 */
  38.             cp->re = -(cp->im = sqrt( -cp->im / 2.0 ));
  39.     else if ( im_sign == 0 )
  40.         if ( re_sign > 0 )
  41.             cp->re = sqrt( cp->re );
  42. /*            cp->im = 0.0;    /* 0 already there */
  43.         else    {        /* re_sign < 0 */
  44.             cp->im = sqrt( -cp->re );
  45.             cp->re = 0.0;
  46.             }
  47.     else    {            /* no shortcuts */
  48.         double    ampl = CxAmpl( cp );
  49.  
  50.         cp->im = sqrt( (ampl - cp->re) /2.0 );
  51.  
  52.         if ( im_sign > 0 )
  53.             cp->re = sqrt( (ampl + cp->re) / 2.0 );
  54.         else            /* im_sign < 0 */
  55.             cp->re = -sqrt( (ampl + cp->re) / 2.0 );
  56.         }
  57.  
  58.     return cp;
  59.     }
  60.