/*							cmplx.c * *	Complex number arithmetic *      This version is for C9X. * * * * SYNOPSIS: * * typedef struct { *      double r;     real part *      double i;     imaginary part *     }cmplx; * * cmplx a, b, c; * * c = cadd( a, b );     c = b + a * c = csub( a, b );     c = b - a * c = cmul( a, b );     c = b * a * c = cdiv( a, b );     c = b / a * c = cneg( a );           c = -a * cmov( b, c );        c = b * * * * DESCRIPTION: * * Addition: *    c.r  =  b.r + a.r *    c.i  =  b.i + a.i * * Subtraction: *    c.r  =  b.r - a.r *    c.i  =  b.i - a.i * * Multiplication: *    c.r  =  b.r * a.r  -  b.i * a.i *    c.i  =  b.r * a.i  +  b.i * a.r * * Division: *    d    =  a.r * a.r  +  a.i * a.i *    c.r  = (b.r * a.r  + b.i * a.i)/d *    c.i  = (b.i * a.r  -  b.r * a.i)/d * ACCURACY: * * In DEC arithmetic, the test (1/z) * z = 1 had peak relative * error 3.1e-17, rms 1.2e-17.  The test (y/z) * (z/y) = 1 had * peak relative error 8.3e-17, rms 2.1e-17. * * Tests in the rectangle {-10,+10}: *                      Relative error: * arithmetic   function  # trials      peak         rms *    DEC        cadd       10000       1.4e-17     3.4e-18 *    IEEE       cadd      100000       1.1e-16     2.7e-17 *    DEC        csub       10000       1.4e-17     4.5e-18 *    IEEE       csub      100000       1.1e-16     3.4e-17 *    DEC        cmul        3000       2.3e-17     8.7e-18 *    IEEE       cmul      100000       2.1e-16     6.9e-17 *    DEC        cdiv       18000       4.9e-17     1.3e-17 *    IEEE       cdiv      100000       3.7e-16     1.1e-16 *//*				cmplx.c * complex number arithmetic *//*Cephes Math Library Release 2.3:  March, 1995Copyright 1984, 1995 by Stephen L. Moshier*/#include "complex.h"#include "mconf.h"#ifndef ANSIPROTdouble fabs(), cabs(), sqrt(), atan2(), cos(), sin();double sqrt(), frexp(), ldexp();#endifint isnan();extern double MAXNUM, MACHEP, PI, PIO2, INFINITY;double complex czero = 0.0;double complex cone = 1.0;/*	c = b + a	*/double complexcadd( a, b )     double complex a, b;{  return (creal (b) + creal (a) + (cimag (b) + cimag (a)) * I);}/*	c = b - a	*/double complexcsub( a, b )     double complex a, b;{  return (creal (b) - creal (a) + (cimag (b) - cimag (a)) * I);}/*	c = b * a */double complexcmul( a, b )     double complex a, b;{  return ((creal (b) * creal (a) - cimag (b) * cimag (a))	  + (creal (b) * cimag (a) + cimag (b) * creal (a)) * I);}/*	c = b / a */double complexcdiv( a, b )     double complex a, b;{  double y, p, q, w;  y = creal (a) * creal (a) + cimag (a) * cimag (a);  p = creal (b) * creal (a) + cimag (b) * cimag (a);  q = cimag (b) * creal (a) - creal (b) * cimag (a);  if( y < 1.0 )    {      w = MAXNUM * y;      if ((fabs(p) > w) || (fabs(q) > w) || (y == 0.0))	{	  mtherr( "cdiv", OVERFLOW );	  return (MAXNUM + MAXNUM * I);	}    }  return (p/y + (q/y) * I);}/*							cabs() * *	Complex absolute value * * * * SYNOPSIS: * * double cabs(); * double complex z; * double a; * * a = cabs( z ); * * * * DESCRIPTION: * * * If z = x + iy * * then * *       a = sqrt( x^2 + y^2 ). *  * Overflow and underflow are avoided by testing the magnitudes * of x and y before squaring.  If either is outside half of * the floating point full scale range, both are rescaled. * * * ACCURACY: * *                      Relative error: * arithmetic   domain     # trials      peak         rms *    DEC       -30,+30     30000       3.2e-17     9.2e-18 *    IEEE      -10,+10    100000       2.7e-16     6.9e-17 *//*Cephes Math Library Release 2.1:  January, 1989Copyright 1984, 1987, 1989 by Stephen L. MoshierDirect inquiries to 30 Frost Street, Cambridge, MA 02140*/#ifdef UNK#define PREC 27#define MAXEXPD 1024#define MINEXPD -1077#endif#ifdef DEC#define PREC 29#define MAXEXPD 128#define MINEXPD -128#endif#ifdef IBMPC#define PREC 27#define MAXEXPD 1024#define MINEXPD -1077#endif#ifdef MIEEE#define PREC 27#define MAXEXPD 1024#define MINEXPD -1077#endif#if 1doublecabs( z )     double complex z;{  double x, y, b, re, im;  int ex, ey, e;#ifdef INFINITIES/* Note, cabs(INFINITY,NAN) = INFINITY. */  if(creal (z) == INFINITY || cimag (z) == INFINITY     || creal (z) == -INFINITY || cimag (z) == -INFINITY )    return( INFINITY );#endif#ifdef NANS  if (isnan(creal(z)))    return (creal(z));  if(isnan(cimag(z)))    return(cimag(z));#endif  re = fabs (creal(z));  im = fabs (cimag(z));  if (re == 0.0)    return (im);  if (im == 0.0)    return (re);  /* Get the exponents of the numbers */  x = frexp( re, &ex );  y = frexp( im, &ey );  /* Check if one number is tiny compared to the other */  e = ex - ey;  if (e > PREC)    return (re);  if (e < -PREC)    return (im);  /* Find approximate exponent e of the geometric mean. */  e = (ex + ey) >> 1;  /* Rescale so mean is about 1 */  x = ldexp( re, -e );  y = ldexp( im, -e );		  /* Hypotenuse of the right triangle */  b = sqrt( x * x  +  y * y );  /* Compute the exponent of the answer. */  y = frexp( b, &ey );  ey = e + ey;  /* Check it for overflow and underflow. */  if (ey > MAXEXPD)    {      mtherr ("cabs", OVERFLOW);      return (INFINITY);    }  if (ey < MINEXPD)    return (0.0);  /* Undo the scaling */  b = ldexp (b, e);  return (b);}#endif /* 1 *//*							csqrt() * *	Complex square root * * * * SYNOPSIS: * * double complex csqrt(); * double complex z, w; * * w = csqrt (z); * * * * DESCRIPTION: * * * If z = x + iy,  r = |z|, then * *                       1/2 * Re w  =  [ (r + x)/2 ]   , * *                       1/2 * Im w  =  [ (r - x)/2 ]   . * * Cancellation error in r-x or r+x is avoided by using the * identity  2 Re w Im w  =  y. * * Note that -w is also a square root of z.  The root chosen * is always in the right half plane and Im w has the same sign as y. * * * * ACCURACY: * *                      Relative error: * arithmetic   domain     # trials      peak         rms *    DEC       -10,+10     25000       3.2e-17     9.6e-18 *    IEEE      -10,+10   1,000,000     2.9e-16     6.1e-17 * */double complexcsqrt (z)     double complex z;{  double complex w;  double x, y, r, t, scale;  x = creal (z);  y = cimag (z);  if (y == 0.0)    {      if (x == 0.0)	{	  w = 0.0 + y * I;	}      else	{	  r = fabs (x);	  r = sqrt (r);	  if (x < 0.0)	    {	      w = 0.0 + r * I;	    }	  else	    {	      w = r + y * I;	    }	}      return (w);    }  if (x == 0.0)    {      r = fabs (y);      r = sqrt (0.5*r);      if (y > 0)	w = r + r * I;      else	w = r - r * I;      return (w);    } /* Rescale to avoid internal overflow or underflow.  */ if ((fabs(x) > 4.0) || (fabs(y) > 4.0))   {     x *= 0.25;     y *= 0.25;     scale = 2.0;   } else   {#if 1     x *= 1.8014398509481984e16;  /* 2^54 */     y *= 1.8014398509481984e16;     scale = 7.450580596923828125e-9; /* 2^-27 */#else     x *= 4.0;     y *= 4.0;     scale = 0.5;#endif   }  w = x + y * I;  r = cabs(w);  if( x > 0 )    {      t = sqrt( 0.5 * r + 0.5 * x );      r = scale * fabs( (0.5 * y) / t );      t *= scale;    }  else    {      r = sqrt( 0.5 * r - 0.5 * x );      t = scale * fabs( (0.5 * y) / r );      r *= scale;    }  if (y < 0)    w = t - r * I;  else    w = t + r * I;  return (w);}doublehypot( x, y )     double x, y;{  double complex z;  z = x + y * I;return (cabs(z));}