cmplx.c

gridgen是一款强大的网格生成程序

/*							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, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/


#include "complex.h"
#include "mconf.h"

#ifndef ANSIPROT
double fabs(), cabs(), sqrt(), atan2(), cos(), sin();
double sqrt(), frexp(), ldexp();
#endif

int isnan();

extern double MAXNUM, MACHEP, PI, PIO2, INFINITY;
double complex czero = 0.0;
double complex cone = 1.0;

/*	c = b + a	*/

double complex
cadd( a, b )
     double complex a, b;
{
  return (creal (b) + creal (a) + (cimag (b) + cimag (a)) * I);
}


/*	c = b - a	*/

double complex
csub( a, b )
     double complex a, b;
{
  return (creal (b) - creal (a) + (cimag (b) - cimag (a)) * I);
}

/*	c = b * a */

double complex
cmul( 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 complex
cdiv( 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, 1989
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct 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 1
double
cabs( 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 complex
csqrt (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);
}


double
hypot( x, y )
     double x, y;
{
  double complex z;

  z = x + y * I;
return (cabs(z));
}