complex.c   [plain text]


/*
 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
 *
 * @APPLE_LICENSE_HEADER_START@
 * 
 * The contents of this file constitute Original Code as defined in and
 * are subject to the Apple Public Source License Version 1.1 (the
 * "License").  You may not use this file except in compliance with the
 * License.  Please obtain a copy of the License at
 * http://www.apple.com/publicsource and read it before using this file.
 * 
 * This Original Code and all software distributed under the License are
 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT.  Please see the
 * License for the specific language governing rights and limitations
 * under the License.
 * 
 * @APPLE_LICENSE_HEADER_END@
 */
 /****************************************************************************
**  File: complex.c
**
**  Contains:  C source code for implementations of floating-point
**             (double) complex functions defined in header file
**             "complex.h" for PowerPC Macintoshes in native mode.
**             Transcendental function algorithms are based on the
**             paper "Branch Cuts for Complex Elementary Functions"
**             by W. Kahan, May 17, 1987, and on Pascal and C source
**             code for the SANE 80-/96-bit extended type by Kenton
**             Hanson and Paul Finlayson, respectively.
**
**            
**  Written by: Jon Okada, SANEitation Engineer, ext. 4-4838
**   
**  Copyright:  c 1987-1993 by Apple Computer, Inc., all rights reserved.
**   
**  Change History (most recent first):
**
**     25 Aug 93   ali    Changed clog to cLog to avoid clashing with the
**                        stream i/o definition clog.
**     14 Jul 93   ali    Added #pragma fenv_access on
**     22 Feb 93   ali    Added a nomaf #pragma.
**     05 Feb 93   JPO    Modified calls to feclearexcept and feraiseexcept
**                        to reflect changes in "fenv.h".
**     18 Dec 92   JPO    First created.
**            
****************************************************************************/

#pragma option nomaf
#pragma fenv_access on

#include "math.h"
#include "complex.h"
#include "fenv.h"

#define Real(z) (__real__ z)
#define Imag(z) (__imag__ z)

/****************************************************************************
   CONSTANTS used by complex functions
****************************************************************************/

#if defined(__BIG_ENDIAN__)
#define HEXDOUBLE(hi, lo) { { hi, lo } }
#elif defined(__LITTLE_ENDIAN__)
#define HEXDOUBLE(hi, lo) { { lo, hi } }
#else
#error Unknown endianness
#endif

static const union {              /* +infinity */
   long int ival[2];
   double dval;
   } FPKINF = HEXDOUBLE(0x7ff00000,0x00000000);

static const union {              /* sqrt(2.0) */
   long int ival[2];
   double dval;
   } FPKSQT2 = HEXDOUBLE(0x3ff6a09e,0x667f3bcd);

static const union {              /* sqrt(2.0) + 1.0 to double */
   long int ival[2];
   double dval;
   } FPKR2P1 = HEXDOUBLE(0x4003504f,0x333f9de6);

static const union {              /* sqrt(2.0) + 1.0 - FPKR2P1 to double */
   long int ival[2];
   double dval;
   } FPKT2P1 = HEXDOUBLE(0x3ca21165,0xf626cdd6);

static const union {              /* pi/2.0 */
   long int ival[2];
   double dval;
   } FPKPI2 = HEXDOUBLE(0x3ff921fb,0x54442d18);

static const union {              /* pi */
   long int ival[2];
   double dval;
   } FPKPI = HEXDOUBLE(0x400921fb,0x54442d18);

static const union {              /* underflow threshold / round threshold */
   long int ival[2];
   double dval;
   } FPKLOVERE = HEXDOUBLE(0x03600000,0x00000000);

static const union {              /* sqrt(0.5) */
   long int ival[2];
   double dval;
   } FPKSQRTHALF = HEXDOUBLE(0x3fe6a09e,0x667f3bcd);

static const union {              /* log(2.0) */
   long int ival[2];
   double dval;
   } FPKLOG2 = HEXDOUBLE(0x3fe62e42,0xfefa39ef);

static const union {              /* exp(709.0) */
   long int ival[2];
   double dval;
   } FPKEXP709 = HEXDOUBLE(0x7fdd422d,0x2be5dc9b);

static const union {              /* asinh(nextafterd(+infinity,0.0))/4.0 */
   long int ival[2];
   double dval;
   } FPKASINHOM4 = HEXDOUBLE(0x406633ce,0x8fb9f87e);

static const union {              /* 2.0^(-33) */
   long int ival[2];
   double dval;
   } FPK2M33 = HEXDOUBLE(0x3de00000,0x00000000);

static const union {              /* 2.0^54 */
   long int ival[2];
   double dval;
   } FPK2P54 = HEXDOUBLE(0x43500000,0x00000000);

#if 0
static const union {              /* 2.0^-(1022) */
   long int ival[2];
   double dval;
   } FPKTINY = HEXDOUBLE(0x00100000,0x00000000);
#endif

static const union {              /* sqrt(nextafterd(+infinity,0.0))/4.0 */
   long int ival[2];
   double dval;
   } FPKTHETA = HEXDOUBLE(0x5fcfffff,0xffffffff);

static const union {              /* 4.0/sqrt(nextafterd(+infinity,0.0)) */
   long int ival[2];
   double dval;
   } FPKRHO = HEXDOUBLE(0x20100000,0x00000000);



/****************************************************************************
   double complex xdivc(double x, double complex y) returns (real x) / (complex y).
   
   CONSTANTS:   FPKINF = +infinity
****************************************************************************/

static
double complex xdivc( double x, double complex y )   /* returns (real x) / (complex y) */
{
   double complex      z;
   double   r, denom;
   
   if ( fabs(Real(y)) >= fabs(Imag(y)) ) {     /* |Real(y)| >= |Imag(y)| */
      if (fabs(Real(y)) == FPKINF.dval) {   /* Imag(y) and Real(y) are infinite */
         Real(z) = copysign(0.0,Real(y));
         Imag(z) = copysign(0.0,-Imag(y));
      }
      else {                             /* |Real(y)| >= finite |Imag(y)| */
         r = Imag(y)/Real(y);
         denom = Real(y) + Imag(y)*r;
         Real(z) = x/denom;
         Imag(z) = (-x*r)/denom;
      }
   }
   
   else {                                /* |Real(y)| !>= |Imag(y)| */
      r = Real(y)/Imag(y);
      denom = r*Real(y) + Imag(y);
      Real(z) = (r*x)/denom;
      Imag(z) = -x/denom;
   }
   
   return z;
}

#ifdef __ppc__
static double __cabs ( double complex z )
{
   double a,b,s,t;
   fenv_t env;
   int   clre,clim,ifoo;
   
   clre = fpclassify(Real(z));
   clim = fpclassify(Imag(z));
   
   if ((clre < FP_NORMAL) || (clim < FP_NORMAL)) {
      return (fabs(Real(z)) + fabs(Imag(z))); /* Real(z) or Imag(z) is NaN, INF, or zero */
   }
   
   else {                        /* both components of z are finite, nonzero */
      ifoo = feholdexcept(&env);         /* save environment, clear flags */
      ifoo = fesetround(FE_TONEAREST);   /* set default rounding */
      a = fabs(Real(z));                    /* work with absolute values */
      b = fabs(Imag(z));
      s = 0.0;
      if (a < b) {                       /* order a >= b */
         t = a;
         a = b;
         b = t;
      }
      t = a - b;                         /* magnitude difference */
      
      if (t != a) {                      /* b not negligible relative to a */
         if (t > b) {                    /* a - b > b */
            s = a/b;
            s += sqrt(1.0 + s*s);
         }
         else {                          /* a - b <= b */
            s = t/b;
            t = (2.0 + s)*s;
            s = ((FPKT2P1.dval+t/(FPKSQT2.dval+sqrt(2.0+t)))+s)+FPKR2P1.dval;
         }
         
         s = b/s;                        /* may spuriously underflow */
         feclearexcept(FE_UNDERFLOW);
      }
      
      feupdateenv(&env);                 /* restore environment */
      return (a + s);                    /* deserved overflow occurs here */
   }                                     /* finite, nonzero case */
}   
// XXX Revert __cabs call to cabs and discard __cabs above post-10.3 (chicken/egg situation)
float cabsf( float complex z ){ return (float) __cabs((double complex) z); }
#endif

/****************************************************************************
   double carg(double complex z) returns the argument (in radians) of its
   complex argument z.  The algorithm is from Kahan's paper.
   
   The argument of a complex number z = x + i*y is defined as atan2(y,x)
   for finite x and y.
   
   CONSTANTS:  FPKPI2 = pi/2.0 to double precision
               FPKPI = pi to double precision
   
   Calls:  fpclassify, copysign, fabs, atan
****************************************************************************/

double carg ( double complex z )
{
   double a,b,argr;
   int   clre,clim;
   
   a = Real(z);
   b = Imag(z);
   clre = fpclassify(a);
   clim = fpclassify(b);
   
   if ((clre == FP_ZERO) && (clim == FP_ZERO)) {  /* zero real and imag parts */
      a = copysign(1.0,a);
   }
   
   if ((clre == FP_INFINITE) && (clim == FP_INFINITE)) {   /* both parts INF */
      a = copysign(1.0,a);
      b = copysign(1.0,b);
   }
         
   if (fabs(b) > fabs(a))               /* |imag| > |real| */
      argr = copysign(FPKPI2.dval,b) - atan(a/b);

   else {
      if (a < 0.0)                      /* |real| >= |imag| */
         argr = copysign(FPKPI.dval,b) + atan(b/a);
      else
         argr = atan(b/a);
   }
   
   return argr;
}
   
   

/****************************************************************************
   static double cssqs(double complex z, long int *k) returns rho = |z/(2^*k)|^2
   with scale factor *k set to avoid overflow/underflow.  Algorithm is
   from the Kahan paper.
   
   CONSTANTS:  FPKINF = +infinity
               FPKLOVERE = (double underflow threshold)/(double round threshold)
                         = (4.0*nextafterd(1.0,0.0)/nextafterd(INF,0.0))
						   /(1.0 - nextafterd(1.0,0)
   
   Calls:  fabs, logb, scalb, fmax, feholdexcept, fetestexcept, feclearexcept,
           and feupdateenv.
   
   Called by:  csqrt and clog.
****************************************************************************/

static double cssqs ( double complex z, long int *k)
{
   double a,b,rho;
   fenv_t env;
   long int iscale;
   int ifoo;
   
   iscale = 0;
   a = fabs(Real(z));
   b = fabs(Imag(z));
   ifoo = feholdexcept(&env);          /* save environment, clr flags */
   rho = a*a + b*b;                    /* preliminary result */
   
   if ((a == FPKINF.dval) || (b == FPKINF.dval)) {
      rho = FPKINF.dval;       /* +INF result if Real(z) or Imag(z) is infinite */
   }
   
   else if (fetestexcept(FE_OVERFLOW) || (fetestexcept(FE_UNDERFLOW) && (rho < FPKLOVERE.dval))) {
      iscale = (long int) logb(fmax(a,b));       /* scaling necessary */
      a = scalb(a,-iscale);
      b = scalb(b,-iscale);
      rho = a*a + b*b;                 /* re-calculate scaled square magnitude */
   }
   
   feclearexcept(FE_OVERFLOW + FE_UNDERFLOW);
   feupdateenv(&env);                  /* restore environment */
   *k = iscale;                        /* store scale value */
   return (rho);
}
   
   
   
/****************************************************************************
   double complex csqrt(double complex z) returns the complex square root of its argument.
   The algorithm, which is from the Kahan paper, uses the following 
   identities:
   
      sqrt(x + i*y) = sqrt((|z| + Real(z))/2) + i*sqrt((|z| - Real(z))/2) and
      sqrt(x - i*y) = sqrt((|z| + Real(z))/2) - i*sqrt((|z| - Real(z))/2),
   
   where y is positive and x may be positive or negative.
   
   CONSTANTS:  FPKINF = +infinity
   
   Calls:  cssqs, scalb, fabs, sqrt, copysign.
****************************************************************************/

double complex csqrt ( double complex z )   
{
   double rho,x,y;
   double complex w;
   long int k;

   rho = cssqs(z,&k);                   /* scaled square magnitude */
   
   if (Real(z) == Real(z))
      rho = scalb(fabs(Real(z)),-k) + sqrt(rho);   /* scaled |Real(z)| + |z| */
   
   if (k%2)                             /* k is odd */
      k = (k-1)/2;
   else {                               /* k is even */
      k = k/2 - 1;
      rho = rho + rho;
   }
   
   rho = scalb(sqrt(rho),k);            /* sqrt((|Real(z)| + |z|)/2) */
   x = rho;
   y = Imag(z);
   
   if (rho != 0) {
      if (fabs(y) != FPKINF.dval)
         y = (y/rho)*0.5;          /* signbit(Imag(z))*sqrt((|z|-|Real(z)|)/2 */
      if (Real(z) < 0) {
         x = fabs(y);
         y = copysign(rho,Imag(z));
      }
   }
   
   Real(w) = x;
   Imag(w) = y;
   return w;
}
      


/****************************************************************************
   double complex clog(double complex z) returns the complex natural logarithm of its
   argument.  The algorithm, which is from the Kahan paper, avoids spurious
   underflow or overflow.
   
   CONSTANTS:  FPKSQRTHALF = sqrt(0.5) to double precision
               FPKLOG2 = log(2.0) to double precision

   Calls:  fabs, cssqs, log1p, log, carg. 
****************************************************************************/

double complex clog ( double complex z )   
{
   double rho,dmax,dmin,temp;
   double complex w;
   long int k;
   
   dmax = fabs(Real(z));  /* order real and imaginary parts of z by magnitude */
   dmin = fabs(Imag(z));
   if (dmax < dmin) {
      temp = dmax;
      dmax = dmin;
      dmin = temp;
   }
   
   rho = cssqs(z,&k);                     /* scaled |z*z| */
   
   if ((k == 0) && (dmax > FPKSQRTHALF.dval) && ((dmax <= 1.25) || (rho < 3.0)))
      Real(w) = log1p((dmax - 1.0)*(dmax + 1.0) + dmin*dmin)*0.5; /* |z| near 1.0 */
   else
      Real(w) = log(rho)*0.5 + k*FPKLOG2.dval;   /* more naive approximation */
   
   Imag(w) = carg(z);                        /* imaginary part of logarithm */

   return w;
}
   
   
   
/****************************************************************************
   static double coshmul(double x, double y) returns y*cosh(x) while
   avoiding spurious overflow and invalid exceptions.
   
   CONSTANTS:  FPKEXP709 = exp(709.0) in double precision
   
   Calls:  cosh, exp, fpclassify, fabs, and scalb.
   
   Called by:  csin, ccos, csinh, and ccosh.
****************************************************************************/

static double coshmul ( double x, double y )
{
   double absx, result;
   
   absx = fabs(x);
   if (absx <= 709.0) {                  /* cosh(x) is finite */
      return y*cosh(x);
   }
   
   else if (fpclassify(x) < FP_ZERO) {   /* x is NaN or infinite */
      return (y*absx);
   }
   
   else if (absx > 1460.0) {             /* probable overflow case */
      return (scalb(y,2100));
   }
   
   else {                        /* cosh(x) overflows but y*cosh(x) may not */
      result = (0.5 * FPKEXP709.dval) * y;   /* initialize result to cosh(709) */
      absx -= 709.0;                     /* reduce exponential argument */
      while (absx > 709.0) {             /* exponential reduction loop */
         result *= FPKEXP709.dval;
         absx -= 709.0;
      }
      return (result*exp(absx));         /* final multiplication */
   }
}

   

/****************************************************************************
   static double sinhmul(double x, double y) returns y*sinh(x) while
   avoiding spurious overflow and invalid exceptions.
   
   CONSTANTS:  FPKEXP709 = exp(709.0) in double precision
   
   Calls:  sinh, exp, fpclassify, fabs, and scalb.
   
   Called by:  csin, ccos, csinh, and ccosh.
****************************************************************************/

static double sinhmul ( double x, double y )
{
   double absx, result;
   
   absx = fabs(x);
   if (absx <= 709.0) {                  /* sinh(x) is finite */
      return y*sinh(x);
   }
   
   else if (fpclassify(x) < FP_ZERO) {   /* x is NaN or infinite */
      return (y*x);
   }
   
   else if (absx > 1460.0) {             /* probable overflow case */
      return (scalb(y,2100));
   }
   
   else {                       /* sinh(x) overflows but y*sinh(x) may not */
      result = (0.5*FPKEXP709.dval)*y;  /* initialize result to |sinh(709)| */
      absx -= 709.0;                    /* reduce exponential argument */
      if (signbit(x) != 0) result = -result;   /* take care of sign of result */
      while (absx > 709.0) {            /* exponential reduction loop */
         result *= FPKEXP709.dval;
         absx -= 709.0;
      }
      return (result*exp(absx));        /* final multiplication */
   }
}

   

/****************************************************************************
   double complex csin(double complex z) returns the complex sine of its argument.  The 
   algorithm is based upon the identity:
   
      sin(x + i*y) = sin(x)*cosh(y) + i*cos(x)*sinh(y).
      
   Signaling of spurious overflows, underflows, and invalids is avoided where
   possible.

   Calls:  sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and
           feupdateenv.
****************************************************************************/

double complex csin ( double complex z )   
{
   fenv_t env;
   int ifoo;
   double sinval, cosval;
   double complex w;
   
   ifoo = feholdexcept(&env);         /* save environment, clear flags */
   sinval = sin(Real(z));                /* sine of real part of argument */
   cosval = cos(Real(z));                /* cosine of real part of argument */
   Real(w) = coshmul(Imag(z),sinval);       /* real result = sinval*cosh(Imag(z)) */
   Imag(w) = sinhmul(Imag(z),cosval);       /* imag result = cosval*sinh(Imag(z)) */
   feclearexcept(FE_UNDERFLOW);       /* underflows don't matter */
   feupdateenv(&env);                 /* update environment */
   return w;
}



/****************************************************************************
   double complex ccos(double complex z) returns the complex cosine of its argument.  The 
   algorithm is based upon the identity:
   
      cos(x + i*y) = cos(x)*cosh(y) - i*sin(x)*sinh(y).
      
   Signaling of spurious overflows, underflows, and invalids is avoided where
   possible.

   Calls:  sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and
           feupdateenv.
****************************************************************************/

double complex ccos ( double complex z )   
{
   fenv_t env;
   int ifoo;
   double sinval, cosval;
   double complex w;
   
   ifoo = feholdexcept(&env);         /* save environment, clear flags */
   sinval = sin(Real(z));                /* sine of real part of argument */
   cosval = cos(Real(z));                /* cosine of real part of argument */
   Real(w) = coshmul(Imag(z),cosval);       /* real result = cosval*cosh(Imag(z)) */
   Imag(w) = sinhmul(Imag(z),-sinval);      /* imag result = -sinval*sinh(Imag(z)) */
   feclearexcept(FE_UNDERFLOW);       /* underflows don't matter */
   feupdateenv(&env);                 /* update environment */
   return w;
}



/****************************************************************************
   double complex csinh(double complex z) returns the complex hyperbolic sine of its
   argument.  The algorithm is based upon the identity:
   
      sinh(x + i*y) = cos(y)*sinh(x) + i*sin(y)*cosh(x).
      
   Signaling of spurious overflows, underflows, and invalids is avoided where
   possible.

   Calls:  sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and
           feupdateenv.
****************************************************************************/

double complex csinh ( double complex z )   
{
   fenv_t env;
   int ifoo;
   double sinval, cosval;
   double complex w;
   
   ifoo = feholdexcept(&env);         /* save environment, clear flags */
   sinval = sin(Imag(z));                /* sine of imaginary part of argument */
   cosval = cos(Imag(z));                /* cosine of imaginary part of argument */
   Real(w) = sinhmul(Real(z),cosval);       /* real result = cosval*sinh(Real(z)) */
   Imag(w) = coshmul(Real(z),sinval);       /* imag result = sinval*cosh(Real(z)) */
   feclearexcept(FE_UNDERFLOW);       /* underflows don't matter */
   feupdateenv(&env);                 /* update environment */
   return w;
}



/****************************************************************************
   double complex ccosh(double complex z) returns the complex hyperbolic cosine of its
   argument.  The algorithm is based upon the identity:
   
      cosh(x + i*y) = cos(y)*cosh(x) + i*sin(y)*sinh(x).
      
   Signaling of spurious overflows, underflows, and invalids is avoided where
   possible.

   Calls:  sin, cos, coshmul, sinhmul, feholdexcept, feclearexcept, and
           feupdateenv.
****************************************************************************/

double complex ccosh ( double complex z )   
{
   fenv_t env;
   int ifoo;
   double sinval, cosval;
   double complex w;
   
   ifoo = feholdexcept(&env);         /* save environment, clear flags */
   sinval = sin(Imag(z));                /* sine of imaginary part of argument */
   cosval = cos(Imag(z));                /* cosine of imaginary part of argument */
   Real(w) = coshmul(Real(z),cosval);       /* real result = cosval*cosh(Real(z)) */
   Imag(w) = sinhmul(Real(z),sinval);       /* imag result = sinval*sinh(Real(z)) */
   feclearexcept(FE_UNDERFLOW);       /* underflows don't matter */
   feupdateenv(&env);                 /* update environment */
   return w;
}



/****************************************************************************
   double complex cexp(double complex z) returns the complex exponential of its
   argument.  The algorithm is based upon the identity:
   
      exp(x + i*y) = cos(y)*exp(x) + i*sin(y)*exp(x).
      
   Signaling of spurious overflows and invalids is avoided where
   possible.
   
   CONSTANTS:  FPKEXP709 = exp(709.0) to double precision

   Calls:  sin, cos, exp, fpclassify, and scalb.
   
   Called by:  cepwry, cxpwri, cxpwre, and cxpwry
****************************************************************************/

double complex cexp ( double complex z )   
{
   double sinval, cosval, expval, exparg;
   double complex w;
   
   sinval = sin(Imag(z));                /* sine of imaginary part of argument */
   cosval = cos(Imag(z));                /* cosine of imaginary part of argument */
   
   if (Real(z) <= 709.0) {               /* exp(Real(z)) is finite */
      expval = exp(Real(z));             /* evaluate exponential */
      Real(w) = cosval*expval;           /* real result = cos(Imag(z))*exp(Real(z)) */
      Imag(w) = sinval*expval;           /* imag result = sin(Imag(z))*exp(Real(z)) */
   }
   
   else if (fpclassify(Real(z)) < FP_ZERO) {      /* Real(z) = +INF or a NaN */
      Real(w) = cosval*Real(z);             /* deserved invalid may occur */
      Imag(w) = sinval*Real(z);             /* deserved invalid may occur */
   }
   
   else if (Real(z) > 1460.0) {          /* probable overflow case */
      Real(w) = scalb(cosval,2100);
      Imag(w) = scalb(sinval,2100);
   }
   
   else {           /* exp(Real(z)) overflows but product with sin or cos may not */
      Real(w) = cosval*FPKEXP709.dval;   /* initialize real result */
      Imag(w) = sinval*FPKEXP709.dval;   /* initialize imag result */
      exparg = Real(z) - 709.0;          /* initialize reduced exponent argument */
      while (exparg > 709.0) {        /* exponential reduction loop */
         Real(w) *= FPKEXP709.dval;
         Imag(w) *= FPKEXP709.dval;
         exparg -= 709.0;
      }
      expval = exp(exparg);           /* final exponential value */
      Real(w) *= expval;                 /* final multiplication steps */
      Imag(w) *= expval;
   }
   
   return w;                          /* done */
}
      
/****************************************************************************
   double complex cpow(double complex x, double complex y) returns the complex result of x^y.
   The algorithm is based upon the identity:
   
      x^y = cexp(y*clog(x)).
      
   Calls:  clog, cexp.
****************************************************************************/

double complex cpow ( double complex x, double complex y )     /* (complex x)^(complex y) */
{
   double complex logval,z;
   
   logval = clog(x);                        /* complex logarithm of x */
   Real(z) = Real(y)*Real(logval) - Imag(y)*Imag(logval);  /* multiply by y */
   Imag(z) = Real(y)*Imag(logval) + Imag(y)*Real(logval);
   return (cexp(z));                        /* return complex exponential */
}




/****************************************************************************
   double complex ctanh(double complex x) returns the complex hyperbolic tangent of its
   argument.  The algorithm is from Kahan's paper and is based on the
   identity:
   
      tanh(x+i*y) = (sinh(2*x) + i*sin(2*y))/(cosh(2*x) + cos(2*y))
                  = (cosh(x)*sinh(x)*cscsq + i*tan(y))/(1+cscsq*sinh(x)*sinh(x)),
   
   where cscsq = 1/(cos(y)*cos(y).  For large values of ze.re, spurious
   overflow and invalid signaling is avoided.
   
   CONSTANTS:  FPKASINHOM4 = asinh(nextafterd(+infinity,0.0))/4.0 to double
                  precision
               FPKINF = +infinity
      
   Calls:  tan, sinh, sqrt, fabs, feholdexcept, feraiseexcept, feclearexcept,
           and feupdateenv.
****************************************************************************/

double complex ctanh( double complex z )
{
   fenv_t env;
   int ifoo;
   double tanval, beta, sinhval, coshval, denom;
   double complex w;
   
   ifoo = feholdexcept(&env);            /* save environment and clear flags */
   
   if (fabs(Real(z)) > FPKASINHOM4.dval) {  /* avoid overflow for large |Real(z)| */
      Real(w) = copysign(1.0,Real(z));            /* real result has unit magnitude */
      Imag(w) = copysign(0.0,Imag(z));            /* imag result is signed zero */
      if (fabs(Real(z)) != FPKINF.dval)        /* set inexact for finite Real(z) */
         feraiseexcept(FE_INEXACT);
      feupdateenv(&env);                    /* update environment */
   }                                     /* end large |Real(z)| case */
   
   
   else {                                /* usual case */
      tanval = tan(Imag(z));                   /* evaluate tangent */
      feclearexcept(FE_DIVBYZERO);          /* in case tangent is infinite */
      feupdateenv(&env);                    /* update environment */
      beta = 1.0 + tanval*tanval;           /* 1/(cos(Imag(z)))^2 */
      sinhval = sinh(Real(z));                 /* evaluate sinh */
      coshval = sqrt(1.0+sinhval*sinhval);  /* evaluate cosh */
      
      if (fabs(tanval) == FPKINF.dval) {    /* infinite tangent */
         Real(w) = coshval/sinhval;
         Imag(w) = 1.0/tanval;
      }
      
      else {                                /* finite tangent */
         denom = 1.0 + beta*sinhval*sinhval;
         Real(w) = beta*coshval*sinhval/denom;
         Imag(w) = tanval/denom;
      }
   }                                     /* end usual case */
   
   return w;                             /* done */
}




/****************************************************************************
   double complex ctan(double complex x) returns the complex hyperbolic tangent of its
   argument.  The algorithm is from Kahan's paper and is based on the
   identity:
   
      tan(x + i*y) = (sin(2*x) + i*sinh(2*y))/(cos(2*x) + cosh(2*y))
                   = (tan(x)+i*cosh(y)*sinh(y)*cscsq)/(1+cscsq*sinh(y)*sinh(y)),
   
   where cscsq = 1/(cos(x)*cos(x).  For large values of ze.im, spurious
   overflow and invalid signaling is avoided.
   
   CONSTANTS:  FPKASINHOM4 = asinh(nextafterd(+infinity,0.0))/4.0 to double
                  precision
               FPKINF = +infinity
      
   Calls:  tan, sinh, sqrt, fabs, feholdexcept, feraiseexcept, feclearexcept,
           and feupdateenv.
****************************************************************************/

double complex ctan( double complex z )
{
   fenv_t env;
   int ifoo;
   double tanval, beta, sinhval, coshval, denom;
   double complex w;
   
   ifoo = feholdexcept(&env);            /* save environment and clear flags */
   
   if (fabs(Imag(z)) > FPKASINHOM4.dval) {  /* avoid overflow for large |Imag(z)| */
      Real(w) = copysign(0.0,Real(z));            /* real result has zero magnitude */
      Imag(w) = copysign(1.0,Imag(z));            /* imag result has unit magnitude */
      if (fabs(Imag(z)) != FPKINF.dval)        /* set inexact for finite Real(z) */
         feraiseexcept(FE_INEXACT);
      feupdateenv(&env);                    /* update environment */
   }                                     /* end large |Real(z)| case */
   
   
   else {                                /* usual case */
      tanval = tan(Real(z));                   /* evaluate tangent */
      feclearexcept(FE_DIVBYZERO);          /* in case tangent is infinite */
      feupdateenv(&env);                    /* update environment */
      beta = 1.0 + tanval*tanval;           /* 1/(cos(Real(z)))^2 */
      sinhval = sinh(Imag(z));                 /* evaluate sinh */
      coshval = sqrt(1.0+sinhval*sinhval);  /* evaluate cosh */
      
      if (fabs(tanval) == FPKINF.dval)   {  /* infinite tangent */
         Real(w) = 1.0/tanval;
         Imag(w) = coshval/sinhval;
      }
      
      else {                                /* finite tangent */
         denom = 1.0 + beta*sinhval*sinhval;
         Real(w) = tanval/denom;
         Imag(w) = beta*coshval*sinhval/denom;
      }
   }                                     /* end usual case */
   
   return w;                             /* done */
}




/****************************************************************************
   static double arcsinh(double x) returns the real inverse hyperbolic
   sine of its real argument.  The implementation is based on the identity:
   
      arcsinh(x) = log(x + sqrt(x*x + 1))
                 = log1p(x + x/(1/x + sqrt(1 + 1/(x*x))))
               
   CONSTANTS:  FPK2M33 = 2.0^(-33)
               FPK2P54 = 2.0^54
               FPKLOG2 = log(2.0) to double precision
               FPKTINY = 2.0^(-1022)
   
   Calls:  fpclassify, fabs, log, log1p, sqrt, copysign, feraiseexcept.
   
   Called by:  casinh, casin, cacosh, and cacos.
****************************************************************************/

static double arcsinh ( double x )
{
   double y;
   int iclass;
   
   iclass = fpclassify(x);                  /* classify argument */
   
   if (iclass == FP_NORMAL) {   
      y = fabs(x);                          /* normal argument */
      
      if (y < FPK2M33.dval)                    /* tiny normal argument */
         feraiseexcept(FE_INEXACT);            /* raise inexact and return arg */
   
      else if (y > FPK2P54.dval)               /* huge magnitude argument */
         y = log(y) + FPKLOG2.dval;            /* prevents spurious overflows */
   
      else                                     /* usual case */
         y = log1p(y + y/(1.0/y + sqrt(1.0 + 1.0/(y*y))));
   
      return copysign(y,x);                    /* fix sign */
   }                                        /* end normal case */
   
   else if (iclass < FP_NORMAL)             /* zero, INF, or NaN argument */
      return (x + x);                          /* propagate argument */
      
   else {                                   /* subnormal argument */
      feraiseexcept(FE_UNDERFLOW + FE_INEXACT); /* signal underflow/inex */
      return (x);                           /* return argument */
   }
}
      
   

/****************************************************************************
   double complex casin(double complex z) returns the complex inverse sine of its
   argument.  The algorithm is from Kahan's paper and is based on the
   formulae:
   
       real(casin(z)) = atan (real(z)/real(csqrt(1.0-z)*csqrt(1.0+z)))
       imag(casin(z)) = arcsinh(imag(csqrt(1.0-cconj(z))*csqrt(1.0+z)))
      
   Calls:  arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv.
****************************************************************************/

double complex casin ( double complex z )
{
   double complex zp1, zm, zm1, w;
   fenv_t env;
   int ifoo;
   
   Real(zp1) = 1.0 + Real(z);               /* evaluate zp1 = csqrt(1.0+z) */
   Imag(zp1) = Imag(z);
   zp1 = csqrt(zp1);
   
   Real(zm) = 1.0 - Real(z);                /* evaluate zm1 = csqrt(1.0-z) */
   Imag(zm) = -Imag(z);
   zm1 = csqrt(zm);
   
   ifoo = feholdexcept(&env);         /* save environ., clr flags/enables */
   Real(w) = atan(Real(z)/(Real(zp1)*Real(zm1) -Imag(zp1)*Imag(zm1))); /* real result */
   feclearexcept(FE_DIVBYZERO);       /* clear spurious DIVBYZERO exception */
   feupdateenv(&env);                 /* restore environment with new flags */
   
   Imag(zm) = Imag(z);                      /* evaluate zm1 = csqrt(1.0-cconj(z)) */
   zm1 = csqrt(zm);
   
   Imag(w) = arcsinh(Real(zp1)*Imag(zm1) + Imag(zp1)*Real(zm1));  /* imag result */
   
   return w;
}



/****************************************************************************
   double complex casinh(double complex z) returns the complex inverse hyperbolic sine of
   its argument.  The algorithm is from Kahan's paper and is based on the
   formula:
   
       casinh(z) = -i*casin(i*z).
      
   Calls:  arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv.
****************************************************************************/

double complex casinh ( double complex z )   
{
   double complex zp1, zm, zm1, w;
   fenv_t env;
   int ifoo;
   
   Real(zp1) = 1.0 - Imag(z);                /* evaluate zp1 = csqrt(1.0+i*z) */
   Imag(zp1) = Real(z);
   zp1 = csqrt(zp1);
   
   Real(zm) = 1.0 + Imag(z);                 /* evaluate zm1 = csqrt(1.0-cconj(i*z)) */
   Imag(zm) = Real(z);
   zm1 = csqrt(zm);
   
   Real(w) = arcsinh(Real(zp1)*Imag(zm1) + Imag(zp1)*Real(zm1)); /* imag result */
   
   Imag(zm) = -Real(z);                      /* evaluate zm1 = csqrt(1.0-cconj(z)) */
   zm1 = csqrt(zm);
   
   ifoo = feholdexcept(&env);          /* save environ., clr flags/enables */
   Imag(w) = atan(Imag(z)/(Real(zp1)*Real(zm1) -Imag(zp1)*Imag(zm1))); /* real result */
   feclearexcept(FE_DIVBYZERO);        /* clear spurious DIVBYZERO exception */
   feupdateenv(&env);                  /* restore environment with new flags */
   
   return w;
}



/****************************************************************************
   double complex cacos(double complex z) returns the complex inverse cosine of its
   argument.  The algorithm is from Kahan's paper and is based on the
   formulae:
   
       real(cacos(z)) = 2.0*atan(real(csqrt(1.0-z)/real(csqrt(1.0+z))))
       imag(cacos(z)) = arcsinh(imag(csqrt(1.0-z)*csqrt(cconj(1.0+z))))
      
   Calls:  arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv.
****************************************************************************/

double complex cacos ( double complex z )
{
   double complex zp, zp1, zm1, w;
   fenv_t env;
   int ifoo;
   
   Real(zp) = 1.0 + Real(z);                /* zp1 = csqrt(1.0 + z) */
   Imag(zp) = Imag(z);
   zp1 = csqrt(zp);
   
   Real(zm1) = 1.0 - Real(z);               /* zm1 = csqrt(1.0 - z) */
   Imag(zm1) = -Imag(z);
   zm1 = csqrt(zm1);
   
   ifoo = feholdexcept(&env);         /* save environment, clr flags/enables */
   Real(w) = 2.0*atan(Real(zm1)/Real(zp1));    /* real result */
   feclearexcept(FE_DIVBYZERO);       /* clr possible spurious div-by-zero flag */
   feupdateenv(&env);                 /* update environment with new flags */
   
   Imag(zp) = -Imag(z);                     /* zp1 = csqrt(1.0 + cconj(z)) */
   zp1 = csqrt(zp);
   
   Imag(w) = arcsinh(Real(zp1)*Imag(zm1) + Imag(zp1)*Real(zm1)); /* imag result */
   
   return w;
}



/****************************************************************************
   double complex cacosh(double complex z) returns the complex inverse hyperbolic`cosine
   of its argument.  The algorithm is from Kahan's paper and is based on the
   formulae:
   
       real(cacosh(z)) = arcsinh(real(csqrt(cconj(z)-1.0)*csqrt(z+1.0)))
       imag(cacosh(z)) = 2.0*atan(imag(csqrt(z-1.0)/real(csqrt(z+1.0))))
      
   Calls:  arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv.
****************************************************************************/

double complex cacosh ( double complex z )
{
   double complex zp1, zm, zm1, w;
   fenv_t env;
   int ifoo;
   
   Real(zp1) = Real(z) + 1.0;              /* zp1 = csqrt(z + 1.0) */
   Imag(zp1) = Imag(z);
   zp1 = csqrt(zp1);
   
   Real(zm) = Real(z) - 1.0;               /* zm1 = csqrt(z - 1.0) */
   Imag(zm) = Imag(z);
   zm1 = csqrt(zm);
   
   ifoo = feholdexcept(&env);        /* save environment, clr flags/enables */
   Imag(w) = 2.0*atan(Imag(zm1)/Real(zp1));   /* imag result */
   feclearexcept(FE_DIVBYZERO);      /* clr possible spurious div-by-zero flag */
   feupdateenv(&env);                /* restore environment with new flags */
   
   Imag(zm) = -Imag(z);                    /* zm1 = csqrt(cconj(z) - 1.0) */
   zm1 = csqrt(zm);
   
   Real(w) = arcsinh(Real(zp1)*Real(zm1) - Imag(zp1)*Imag(zm1)); /* real result */
   
   return w;
}



/****************************************************************************
   double complex catan(double complex z) returns the complex inverse tangent
   of its argument.  The algorithm is from Kahan's paper and is based on
   the formula:
   
      catan(z) = i*(clog(1.0-i*z) - clog(1+i*z))/2.0.
   
   CONSTANTS:  FPKTHETA = sqrt(nextafterd(+INF,0.0))/4.0
               FPKRHO = 1.0/FPKTHETA
               FPKPI2 = pi/2.0
               FPKPI = pi
            
   Calls:  copysign, fabs, xdivc, sqrt, log, atan, log1p, and carg.
****************************************************************************/

double complex catan ( double complex z )
{
   double complex ctemp, w;
   double t1, t2, xi, eta, beta;
   
   xi = -Imag(z);
   beta = copysign(1.0,xi);               /* copes with unsigned zero */
   
   Imag(z) = -beta*Real(z);                     /* transform real & imag components */
   Real(z) = beta*xi;
   
   if ((Real(z) > FPKTHETA.dval) || (fabs(Imag(z)) > FPKTHETA.dval)) {
      xi = copysign(FPKPI2.dval,Imag(z));    /* avoid spurious overflow */
      ctemp = xdivc(1.0,z);
      eta = Real(ctemp);
   }
      
   else if (Real(z) == 1.0) {
      t1 = fabs(Imag(z)) + FPKRHO.dval;
      xi = log(sqrt(sqrt(4.0 + t1*t1))/sqrt(fabs(Imag(z))));
      eta = 0.5*copysign(FPKPI.dval-atan(2.0/(fabs(Imag(z))+FPKRHO.dval)),Imag(z));
   }
   
   else {                                 /* usual case */
      t2 = fabs(Imag(z)) + FPKRHO.dval;
      t1 = 1.0 - Real(z);
      t2 = t2*t2;
      xi = 0.25*log1p(4.0*Real(z)/(t1*t1 + t2));
      Real(ctemp) = (1.0 - Real(z))*(1.0 + Real(z)) - t2;
      Imag(ctemp) = Imag(z) + Imag(z);
      eta = 0.5*carg(ctemp);
   }
   
   Real(w) = -beta*eta;                      /* fix up signs of result */
   Imag(w) = -beta*xi;
   return w;
}



/****************************************************************************
   double complex catanh(double complex z) returns the complex inverse hyperbolic tangent
   of its argument.  The algorithm is from Kahan's paper and is based on
   the formula:
   
      catan(z) = (clog(1.0 + z) - clog(1 - z))/2.0.
   
   CONSTANTS:  FPKTHETA = sqrt(nextafterd(+INF,0.0))/4.0
               FPKRHO = 1.0/FPKTHETA
               FPKPI2 = pi/2.0
               FPKPI = pi
            
   Calls:  copysign, fabs, xdivc, sqrt, log, atan, log1p, and carg.
****************************************************************************/

double complex catanh( double complex z )   
{
   double complex ctemp, w;
   double t1, t2, xi, eta, beta;
   
   beta = copysign(1.0,Real(z));             /* copes with unsigned zero */
   
   Imag(z) = -beta*Imag(z);                     /* transform real & imag components */
   Real(z) = beta*Real(z);
   
   if ((Real(z) > FPKTHETA.dval) || (fabs(Imag(z)) > FPKTHETA.dval)) {
      eta = copysign(FPKPI2.dval,Imag(z));   /* avoid overflow */
      ctemp = xdivc(1.0,z);
      xi = Real(ctemp);
   }
      
   else if (Real(z) == 1.0) {
      t1 = fabs(Imag(z)) + FPKRHO.dval;
      xi = log(sqrt(sqrt(4.0 + t1*t1))/sqrt(fabs(Imag(z))));
      eta = 0.5*copysign(FPKPI.dval-atan(2.0/(fabs(Imag(z))+FPKRHO.dval)),Imag(z));
   }
   
   else {                                 /* usual case */
      t2 = fabs(Imag(z)) + FPKRHO.dval;
      t1 = 1.0 - Real(z);
      t2 = t2*t2;
      xi = 0.25*log1p(4.0*Real(z)/(t1*t1 + t2));
      Real(ctemp) = (1.0 - Real(z))*(1.0 + Real(z)) - t2;
      Imag(ctemp) = Imag(z) + Imag(z);
      eta = 0.5*carg(ctemp);
   }
   
   Real(w) = beta*xi;                        /* fix up signs of result */
   Imag(w) = -beta*eta;
   return w;
}

/* conj(), creal(), and cimag() are gcc built ins. */

double complex cproj( double complex z )
{
   double complex S;
   double zR = Real(z), zI = Imag(z), scale;

   if (isnan(zR) && isnan(zI))
      return z;
    
   if (isinf(zR) || isinf(zI)) {
      Real(S) = INFINITY;
      Imag(S) = copysign(0.0, zI);
      return S;
   }

   scale = 1.0 + zR * zR + zI * zI;
   Real(S) = (zR + zR)/scale;
   Imag(S) = (zI + zI)/scale;
   return S;
}

float complex cacosf( float complex z ){ return (float complex) cacos((double complex) z); }
float complex casinf( float complex z ){ return (float complex) casin((double complex) z); }
float complex catanf( float complex z ){ return (float complex) catan((double complex) z); }

float complex ccosf( float complex z ){ return (float complex) ccos((double complex) z); }
float complex csinf( float complex z ){ return (float complex) csin((double complex) z); }
float complex ctanf( float complex z ){ return (float complex) ctan((double complex) z); }

float complex cacoshf( float complex z ){ return (float complex) cacosh((double complex) z); }
float complex casinhf( float complex z ){ return (float complex) casinh((double complex) z); }
float complex catanhf( float complex z ){ return (float complex) catanh((double complex) z); }

float complex ccoshf( float complex z ){ return (float complex) ccosh((double complex) z); }
float complex csinhf( float complex z ){ return (float complex) csinh((double complex) z); }
float complex ctanhf( float complex z ){ return (float complex) ctanh((double complex) z); }

float complex cexpf( float complex z ){ return (float complex) cexp((double complex) z); }
float complex clogf( float complex z ){ return (float complex) clog((double complex) z); }

float complex cpowf( float complex x, float complex y ){ return (float complex) cpow((double complex) x, (double complex) y); }
float complex csqrtf( float complex z ){ return (float complex) csqrt((double complex) z); }

float cargf( float complex z ){ return (float) carg((double complex) z); }
float complex cprojf( float complex z ){ return (float complex) cproj((double complex) z); }