complex.h

Go to the documentation of this file.

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */
// define the range of used complex (and imaginary) numbers 
//#pragma STDC CX_LIMITED_RANGE OFF
 
              // _Imaginary_I, _Complex_I, _Imaginary, _Complex, 
              // I are macros expanding to: 
//complex       // to _Complex.
//_Complex_I    // to a constant expression of type const float _Complex, 
              // with the value of the imaginary unit */
              // (that is, a number i such that i.i=-1).
//imaginary     // to _Imaginary
//_Imaginary_I  // to a constant expression of type const float _Imaginary 
              // with the value of the imaginary unit, or expands to 
              // either _Imaginary_I or _Complex_I.  
              // If _Imaginary_I is not defined, I expands to _Complex_I.
// I have commented out this. Otherwise complex.h is not includeable.

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#ifndef _COMPLEX_H
#define _COMPLEX_H

#include <nlibc.h>
#include <math.h>

#define _math_setdouble(U,D) { \
  register union { double d; unsigned i; } u; \
  asm("\tatr0h %0 <"#U">" : "=r" (u.i)); D=u.d; \
}

#define _Complex_I {0., 1.}                    /* imaginary unit */
#define I _Complex_I                           /* imaginary unit */ 

#define M_PI           3.14159265358979323846  /* pi */
#define M_PI_2         1.57079632679489661923  /* pi/2 */
#define M_PI_4         0.78539816339744830962  /* pi/4 */

#define conj(x)  (~(x))
#define creal(x) ((x).re)
#define cimag(x) ((x).im)


/* ------------------------------------------------------------------------------------- */
/*  nan(char *) not implemented, using workaround:            */
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double cnan(void);
#else
#if defined(_uses_cnan_complex_h) || !defined(__cflow_processed)
math_inline double cnan(void)
{
  //  unsigned int nan=0x7FF0000000000001ULL;
  //  return *((double*)&nan);
  register double nan; 
  _math_setdouble( 0x7FF0000000000001 , nan );
 return nan;
}
#endif  // cnan
#endif // Has Main

#define NAN                 (cnan())  // (nan("")    

/*  using HUGE_VAL (=INFINITY) comes with warnings:           */
#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double cinfty(void);
#else
#if defined(_uses_cinfty_complex_h) || !defined(__cflow_processed)
math_inline double cinfty(void)
{
  register double infty; 
  _math_setdouble( 0x7FF0000000000000 , infty );
 return infty;
}
#endif  // cinfty
#endif // Has Main

#define HUGE_VALD           (cinfty())    

/*
                    glibc     ape                         
#define FP_NAN        1   //   3
#define FP_INFINITE   2   //   2
#define FP_ZERO       3   //   4
#define FP_NORMAL     4   //   1
#define FP_SUBNORMAL  5   //   1

#if defined(_uses_cfpclassify_complex_h) || !defined(__cflow_processed)
math_inline int cfpclassify(double x) {
  register unsigned int _emask,_mmask;
  register int ie,ii,ij;
  register union {
    double xr;
    unsigned int xint;
  } xu;
 
  _emask=0x7FF0000000000000ULL;
  _mmask=0x000fffffffffffffULL;
  xu.xr = x;
  ii = (_emask&xu.xint);
  ij = (_mmask&xu.xint);
  where ((ii != _emask) && (ij != 0))  ie = FP_NORMAL;
  where ((ii == _emask) && (ij == 0) ) ie = FP_INFINITE;
  where ((ii == _emask) && (ij != 0) ) ie = FP_NAN;
  where ((ii == 0) && (ij == 0))       ie = FP_ZERO;
  return ie;
}
#endif
*/
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex sincos(double x);
#else
#if defined(_uses_sincos_complex_h) || !defined(__cflow_processed)
math_inline complex sincos(double x)
{
  register double twooverpi, pihalflead, pihalftrail1, pihalftrail2, magic2, f, y, siny, cosy;
  register complex sc;
  register int k;
  register union {
    double z;
    int zint;
  } zu;
  _math_setdouble(  0x3fe45f306dc9c883 , twooverpi   );
  _math_setdouble(  0x4337ffffffffffff , magic2      );
  _math_setdouble(  0x3ff921fb50000000 , pihalflead  );
  _math_setdouble(  0x3e5110b460000000 , pihalftrail1 );
  _math_setdouble(  0x3c91a62633145c07 , pihalftrail2 );
  zu.z = magic2 + twooverpi * x;
  f = zu.z - magic2;
  y = x - f * pihalflead;
  y = y - f * pihalftrail1;
  y = y - f * pihalftrail2;
  siny = sin_1oct(y);
  cosy = cos_1oct(y);
  k = 0x3;
  k = zu.zint & k ;
  where ( k == 3 ) {
    sc.re =  siny;
    sc.im =  cosy;
  }
  where ( k == 0 ) {
    sc.re =  cosy;
    sc.im = -siny;
  }
  where ( k == 1 ) {
    sc.re = -siny;
    sc.im = -cosy;
  }
  where ( k == 2 ) {
    sc.re = -cosy;
    sc.im =  siny;
  }
  return sc;                    // (sin(x), cos(x))
}
#endif  // sincos
#endif // Has Main

/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex sinhcosh(double x);
#else
#if defined(_uses_sinhcosh_complex_h) || !defined(__cflow_processed)
math_inline complex sinhcosh(double x)
{
  register double expx, expxm;
  register complex shch;
  expx  = 0.5 * expm1( x);
  expxm = 0.5 * expm1(-x);
  shch.re = (expx - expxm);
  shch.im = (expx + expxm + 1.0);
  return shch;                   // (sinh(x), cosh(x))
} 
#endif  // sinhcosh
#endif // Has Main

/* ------------------------------------------------------------------------------------- */



/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double cabs(complex z) ;
#else
#if defined(_uses_cabs_complex_h) || !defined(__cflow_processed)
math_inline double cabs(complex z) 
{
  register complex res;
  
  res = hypot(z.re,z.im);
  return res;
}
#endif // cabs()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline double carg(complex z) ;
#else
#if !defined(__cflow_processed) || defined(_uses_carg_complex_h)
math_inline double carg(complex z) 
{
  register complex res;

  res = atan2(z.im,z.re);
  return res;
}
#endif // carg()
#endif // Has Main


/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex cproj(complex z) ;
#else
#if defined(_uses_cproj_complex_h) || !defined(__cflow_processed)
math_inline complex cproj(complex z) 
{
  register complex res;
  register double den;

#if defined(__FAST_MATH__)
  den = z.re * z.re + z.im * z.im +1.0;       //test:  den = z*(~z) + 1.0;
  res.re = (2.0 * z.re)/den;                  //test:  res = (2.0 * z) / den;
  res.im = (2.0 * z.im)/den;

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where (rcls == FP_NAN  && icls == FP_NAN) {
    res = z;
  } 
  else where ((!(rcls == FP_ZERO) && !(rcls == FP_NORMAL)) || (!(icls == FP_ZERO) && !(icls == FP_NORMAL))) {
    res.re = HUGE_VALD;
    res.im = copysign(0.0, z.im);
  }
  else {
    den = z.re * z.re + z.im * z.im + 1.0; 
    res.re = (2.0 * z.re)/den;             
    res.im = (2.0 * z.im)/den;
  }
#endif
  return res;
}
#endif // cproj()
#endif // Has Main


/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex cexp(complex z) ;
#else
#if defined(_uses_cexp_complex_h) || !defined(__cflow_processed) 
math_inline complex cexp(complex z) 
{
  register complex res, sc;
  register double exp_val;

#if defined(__FAST_MATH__)
  sc = sincos(z.im);
  exp_val = exp(z.re);
  res.re = exp_val * sc.im;
  res.im = exp_val * sc.re;


#else
  register double value;
  register int rcls, icls, ecls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((rcls == FP_ZERO) || (rcls == FP_NORMAL)) {              // Real part is finite.  
    where ((icls == FP_ZERO) || (icls == FP_NORMAL)) {              // Imaginary part is finite.
      exp_val = exp(z.re);
      ecls = fpclassify(exp_val);
      sc = sincos(z.im);
      where (ecls == FP_ZERO || ecls == FP_NORMAL) {
      //      where (isfinite(exp_val)) {
        res.re = exp_val * sc.im;
        res.im = exp_val * sc.re;
      }
      else {
        res.re = copysign(exp_val, sc.im);
        res.im = copysign(exp_val, sc.re);
      }
    }
    else { // If the imaginary part is +-inf or NaN and the real part is not +-inf 
           //   the result is NaN + iNaN. 
      res.re = NAN;
      res.im = NAN;
      //#ifdef FE_INVALID
      //feraiseexcept(FE_INVALID);
      //#endif
    }
  }
  else where (rcls == FP_INFINITE) {           // Real part is infinite. 
    where ((icls == FP_ZERO) || (icls == FP_NORMAL)) {  // Imaginary part is finite.
      value = signbit(z.re) ? 0.0 : HUGE_VALD;
      where (icls == FP_ZERO) {                // Imaginary part is 0.0.  
        res.re = value;
        res.im = z.im;
      }
      else {
        sc = sincos(z.im);
        res.re = copysign(value, sc.im);
        res.im = copysign(value, sc.re);
      }
    }
    else where (signbit(z.re) == 0) {
          res.re = HUGE_VALD;
          res.im = NAN;
          //#ifdef FE_INVALID
          //where (icls == FP_INFINITE) feraiseexcept(FE_INVALID);
          //#endif
    }
    else {
      res.re = 0.0;
      res.im = copysign(0.0, z.im);
    }
  }
  else {                       // If the real part is NaN the result is NaN + iNaN.
     res.re = NAN;
     res.im = NAN;
    //#ifdef FE_INVALID
    //where (rcls != FP_NAN || icls != FP_NAN) feraiseexcept(FE_INVALID);
    //#endif
  }

#endif
  return res;
}
#endif // cexp()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex clog(complex z);
#else
#if defined(_uses_clog_complex_h) || !defined(__cflow_processed)
math_inline complex clog(complex z)
{
  register complex res;

#if defined(__FAST_MATH__)
  res.re = log(hypot(z.re, z.im));  //test: res.re = 0.5*log(~z*z);
  res.im = atan2(z.im,z.re);

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where (rcls == FP_ZERO && icls == FP_ZERO) {
    res.im = signbit(z.re) ? M_PI : 0.0;
    res.im = copysign(res.im, z.im);
    res.re = -HUGE_VALD;  //orig.:  res.re = -1.0 / fabs(z.re);  -> exeption,  why ?
  }
  else where (rcls != FP_NAN && icls != FP_NAN) {
    res.re = log(hypot(z.re, z.im));
    res.im = atan2(z.im, z.re);

  } else {
    res.im = NAN;
    where (rcls == FP_INFINITE || icls == FP_INFINITE)  res.re = HUGE_VALD;
    else res.re = NAN;
  }
#endif
  return res;
}
#endif // clog()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex cpow(complex x, complex y) ;
#else
#if defined(_uses_cpow_complex_h) || !defined(__cflow_processed)
math_inline complex cpow(complex x, complex y) 
{
  register complex res;

  res = cexp(y*clog(x));
  return res;
}
#endif // cpow()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex csqrt(complex z)  ;
#else
#if defined(_uses_csqrt_complex_h) || !defined(__cflow_processed) 
math_inline complex csqrt(complex z)  
{
  register complex res;
  register double d, r, s;

#if defined(__FAST_MATH__)
  d = hypot(z.re, z.im);
  where(z.re > 0) {
    r = sqrt(0.5 * d + 0.5 * z.re); //test: r = sqrt(0.5*(d+z.re));  
    s = (0.5 * z.im) / r;
  } else {
    s = sqrt(0.5 * d - 0.5 * z.re); //test: s = sqrt(0.5*(d-z.re));
    r = fabs((0.5 * z.im) / s);
  }
  res.re = r;
  res.im = copysign(s, z.im);

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((rcls == FP_INFINITE || rcls == FP_NAN) || (icls == FP_INFINITE) || (icls == FP_NAN)) {
    where (icls == FP_INFINITE) {
      res.re = HUGE_VALD;
      res.im = z.im;
    }
    else where (rcls == FP_INFINITE) {
      where (z.re < 0.0) {
        res.re = (icls == FP_NAN ? NAN:0);
        res.im = copysign(HUGE_VALD, z.im);
      }
      else {
        res.re = z.re;
        res.im = (icls == FP_NAN ? NAN:copysign(0.0, z.im));
      }
    }
    else {
      res.re = NAN;
      res.im = NAN;
    }
  }
  else {
    where (icls == FP_ZERO) {
      where (z.re < 0.0) {
        res.re = 0.0;
        res.im = copysign(sqrt(-z.re), z.im);
      }
      else {
        res.re = fabs(sqrt(z.re));
        res.im = copysign(0.0, z.im);
      }
    }
    else where (rcls == FP_ZERO) {
      r = sqrt(0.5 * fabs(z.im));
      res.re = copysign(r, z.im);
      res.im = r;
    }
    else {
      d = hypot(z.re, z.im);
      where (z.re > 0) {
        r = sqrt(0.5 * d + 0.5 * z.re);
        s = (0.5 * z.im) / r;
      }
      else {
        s = sqrt(0.5 * d - 0.5 * z.re);
        r = fabs((0.5 * z.im) / s);
      }
      res.re = r;
      res.im = copysign (s, z.im);
    }
  }
#endif
  return res;
}
#endif // csqrt()
#endif // Has Main



/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex cacosh(complex z) ;
#else
#if defined(_uses_cacosh_complex_h) || !defined(__cflow_processed)
math_inline complex cacosh(complex z) 
{
  register complex y,res;

#if defined(__FAST_MATH__)
  y.re = (z.re - z.im)*(z.re + z.im) - 1.0;   //test: y.re = ~z*z-1.0;
  y.im = 2.0 * z.re * z.im;
  y    = csqrt(y);
  y.re += z.re;                               //test: y += z;
  y.im += z.im;
  res  = clog(y);

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((rcls == FP_INFINITE || rcls == FP_NAN) || (icls == FP_INFINITE || icls == FP_NAN)) {
    where (icls == FP_INFINITE) {
      res.re = HUGE_VALD;
      where (rcls == FP_NAN) res.im = NAN;
      else res.im=copysign((rcls==FP_INFINITE ? (z.re<0.0 ? M_PI-M_PI_4:M_PI_4):M_PI_2), z.im);
    }
    else where (rcls == FP_INFINITE) {
      res.re = HUGE_VALD;
      where ((icls == FP_ZERO) || (icls == FP_NORMAL)) res.im = copysign(signbit(z.re) ? M_PI : 0.0, z.im);
      else res.im = NAN;
    }
    else {
      res.re = NAN;
      res.im = NAN;
    }
  }
  else where (rcls == FP_ZERO && icls == FP_ZERO) {
    res.re = 0.0;
    res.im = copysign(M_PI_2, z.im);
  }
  else {
    y.re = (z.re - z.im) * (z.re + z.im) - 1.0;
    y.im = 2.0 * z.re * z.im;
    y = csqrt(y);
    y.re += z.re;
    y.im += z.im;
    res = clog(y);
  }
#endif
  return res;
}
#endif // cacosh()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex casinh(complex z) ;
#else
#if defined(_uses_casinh_complex_h) || !defined(__cflow_processed) 
math_inline complex casinh(complex z) 
{
  register complex y,res;

#if defined(__FAST_MATH__)
  y.re = (z.re - z.im)*(z.re + z.im) + 1.0;  //test:  y.re = ~z*z+1.0 
  y.im = 2.0 * z.re * z.im;
  y    = csqrt(y);
  y.re += z.re;                              //test:  y   += z;
  y.im += z.im;
  res  = clog(y);

#else
  register int rcls, icls;
  rcls=fpclassify(z.re);
  icls=fpclassify(z.im);
  where ((rcls == FP_INFINITE || rcls == FP_NAN) || (icls == FP_INFINITE || icls == FP_NAN)) {
    where (icls == FP_INFINITE) {
      res.re = copysign (HUGE_VALD, z.re);
      where (rcls == FP_NAN) 
        res.im = NAN;
      else 
        res.im = copysign((rcls>=FP_ZERO ? M_PI_2:M_PI_4), z.im);
    } 
    else where (rcls == FP_INFINITE || rcls == FP_NAN) {
      res.re = z.re;
      where (((icls == FP_ZERO || icls == FP_NORMAL) && rcls == FP_INFINITE) || (rcls == FP_NAN && icls == FP_ZERO))
        res.im = copysign(0.0, z.im);
      else
        res.im = NAN;
    }
    else {
      res.re = NAN;
      res.im = NAN;
    }
  }
  else where (rcls == FP_ZERO && icls == FP_ZERO) {
    res = z;
  }
  else {
    y.re = (z.re - z.im)*(z.re + z.im) + 1.0;  
    y.im  = 2.0 * z.re * z.im;
    y     = csqrt(y);
    y.re += z.re;
    y.im += z.im;
    res   = clog(y);
  }
#endif
  return res;
}
#endif // casinh()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex catanh(complex z) ;
#else
#if defined(_uses_catanh_complex_h) || !defined(__cflow_processed) 
math_inline complex catanh(complex z) 
{
  register complex res;
  register double i2, num, den;

#if defined(__FAST_MATH__)
  i2 = z.im * z.im;
  num = 1.0 + z.re;
  num = i2 + num * num;
  den = 1.0 - z.re;
  den = i2 + den * den;
  res.re = 0.25 * (log(num) - log(den)); //test: res.re = 0.25*log(num/den);
  den = 1 - z.re * z.re - i2;
  res.im = 0.5 * atan2(2.0 * z.im, den);

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((rcls == FP_INFINITE || rcls == FP_NAN) || (icls == FP_INFINITE || icls == FP_NAN)) {
    where (icls == FP_INFINITE) {
        res.re = copysign(0.0, z.re);
        res.im = copysign(M_PI_2, z.im);
    }
    else where (rcls == FP_INFINITE || rcls == FP_ZERO) {
      res.re = copysign(0.0, z.re);
      where (icls == FP_ZERO || icls == FP_NORMAL) res.im = copysign(M_PI_2, z.im);
      else res.im = NAN;
    }
    else {
      res.re = NAN;
      res.im = NAN;
    }
  }
  else where (rcls == FP_ZERO && icls == FP_ZERO) res = z;
  else {
    i2 = z.im * z.im;
    num = 1.0 + z.re;
    num = i2 + num * num;
    den = 1.0 - z.re;
    den = i2 + den * den;
    res.re = 0.25 * (log(num) - log(den));
    den = 1 - z.re * z.re - i2;
    res.im = 0.5 * atan2(2.0 * z.im, den);
  }
#endif
  return res;
}
#endif // catanh()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */


#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex cacos(complex z);
#else
#if defined(_uses_cacos_complex_h) || !defined(__cflow_processed) 
math_inline complex cacos(complex z)
{
/* old version:
  register complex res, y;
  y = casin(z);
  res.re = (double) M_PI_2 - y.re;
  res.im = -y.im;
  return res;
*/
  register complex y, res;

#if defined(__FAST_MATH__)
  y = cacosh(z);
  res.re =  y.im;
  res.im = -y.re; // cacos(z) = -I*cacosh(z)

#else 
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where (rcls == FP_NAN || icls == FP_NAN) {
    //    where (z.re == 0.0) res.re = M_PI_2, res.im = z.im;
    where (rcls == FP_ZERO) res.re = M_PI_2, res.im = z.im;
    else where (rcls == FP_INFINITE || icls == FP_INFINITE) {
      res.re = copysign(HUGE_VALD, z.im);
      res.im = NAN;
    } 
    else {
      res.re = NAN;
      res.im = NAN;
    }
  } 
  else {
    y = cacosh(z);
    res.re =  y.im;
    res.im = -y.re;
  }
#endif
  return res;
}
#endif // cacos()
#endif // Has Main


/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex casin(complex z) ;
#else
#if defined(_uses_casin_complex_h) || !defined(__cflow_processed) 
math_inline complex casin(complex z) 
{
  register complex y, res;

#if defined(__FAST_MATH__)
  y.re = -z.im;
  y.im =  z.re;
  y = casinh(y);     //test: res = -I * casinh(I*z);
  res.re =  y.im;
  res.im = -y.re;

#else 
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where (rcls == FP_NAN || icls == FP_NAN) {
    //    where (z.re == 0.0) res = z;
    where (rcls == FP_ZERO) res = z;
    else where (rcls == FP_INFINITE || icls == FP_INFINITE) {
      res.re = NAN;
      res.im = copysign(HUGE_VALD, z.im);
    } 
    else {
      res.re = NAN;
      res.im = NAN;
    }
  } 
  else {
    y.re = -z.im;
    y.im =  z.re;
    y = casinh(y);
    res.re =  y.im;
    res.im = -y.re;
  }
#endif
  return res;
}
#endif // casin()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex catan(complex z);
#else
#if defined(_uses_catan_complex_h) || !defined(__cflow_processed) 
math_inline complex catan(complex z)
{
  register complex res;
  register double r2, num ,den;

#if defined(__FAST_MATH__)
  r2 = z.re * z.re;
  den = 1 - r2 - z.im * z.im;
  res.re = 0.5 * atan2(2.0 * z.re, den);
  num = z.im + 1.0;
  num = r2 + num * num;
  den = z.im - 1.0;
  den = r2 + den * den;
  res.im = 0.25*(log(num)-log(den));  //test:  res.im = 0.25 * log(num/den); 

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((rcls == FP_INFINITE || rcls == FP_NAN) || (icls == FP_INFINITE || icls == FP_NAN)) {
    where (rcls == FP_INFINITE) {
      res.re = copysign(M_PI_2, z.re);
      res.im = copysign(0.0, z.im);
    }
    else where (icls == FP_INFINITE) {
      res.im = copysign(0.0, z.im);
      where (rcls == FP_ZERO || rcls == FP_NORMAL) res.re = copysign(M_PI_2, z.re);
      else res.re = NAN;
    }
    else where (icls == FP_ZERO || icls == FP_INFINITE) {
      res.re = NAN;
      res.im = copysign(0.0, z.im);
    }
    else {
      res.re = NAN;
      res.im = NAN;
    }
  }
  else where (rcls == FP_ZERO && icls == FP_ZERO) res = z;
  else {
    r2 = z.re * z.re;
    den = 1 - r2 - z.im * z.im;
    res.re = 0.5 * atan2(2.0 * z.re, den);
    num = z.im + 1.0;
    num = r2 + num * num;
    den = z.im - 1.0;
    den = r2 + den * den;
    res.im = 0.25*(log(num)-log(den));  //test:  res.im = 0.25 * log(num/den); 
  }
#endif
  return res;
}
#endif // catan()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex ccosh(complex z);
#else
#if defined(_uses_ccosh_complex_h) || !defined(__cflow_processed) 
math_inline complex ccosh(complex z)
{
  register complex res, sc, shch;

#if defined(__FAST_MATH__) 
  sc = sincos(z.im);
  shch = sinhcosh(z.re);
  res.re = shch.im * sc.im;
  res.im = shch.re * sc.re;

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where (rcls == FP_ZERO || rcls == FP_NORMAL) {         /* Real part is finite.  */
    where (icls == FP_ZERO || icls == FP_NORMAL) {          /* Imaginary part is finite.  */
      sc = sincos(z.im);
      shch = sinhcosh(z.re);
      res.re = shch.im * sc.im;
      res.im = shch.re * sc.re;
    }
    else {
      res.im = z.re == 0.0 ? 0.0 : NAN;
      res.re = NAN;                     // res.re = NAN + NAN; ?? (s.u.)
      //#ifdef FE_INVALID
      //where (icls == FP_INFINITE) feraiseexcept(FE_INVALID);
      //#endif
    }
  }
  else where (rcls == FP_INFINITE) {      // Real part is infinite. 
    where (icls == FP_ZERO) {          // Imaginary part is 0.0. 
      res.re = HUGE_VALD;
      res.im = z.im * copysign (1.0, z.re);
    }
    else where (icls ==  FP_NORMAL) {          // Imaginary part is finite. 
      sc = sincos(z.im);
      res.re = copysign (HUGE_VALD, sc.im);
      res.im = (copysign (HUGE_VALD, sc.re) * copysign (1.0, z.re));
    }
    else {
      res.re = HUGE_VALD;
      res.im = NAN;         // res.im = NAN + NAN; ?? The addition raises the invalid exception.
      //#ifdef FE_INVALID
      //where (icls == FP_INFINITE) feraiseexcept(FE_INVALID);
      //#endif
    }
  }
  else {
    res.re = NAN;
    res.im = z.im == 0.0 ? z.im : NAN;
  }
#endif
  return res;
}
#endif // ccosh()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex csinh(complex z) ;
#else
#if defined(_uses_csinh_complex_h) || !defined(__cflow_processed) 
math_inline complex csinh(complex z) 
{
  register complex res, sc, shch;

#if defined(__FAST_MATH__)
  sc = sincos(z.im);
  shch = sinhcosh(z.re);
  res.re = shch.re * sc.im;
  res.im = shch.im * sc.re;

#else
  register int rcls, icls, negate;
  negate = signbit(z.re);
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  z.re = fabs(z.re);
  where (rcls == FP_ZERO || rcls == FP_NORMAL) {      // Real part is finite. 
    where (icls == FP_ZERO || icls == FP_NORMAL) {          // Imaginary part is finite.
      sc = sincos(z.im);
      shch = sinhcosh(z.re);
      res.re = shch.re * sc.im;
      res.im = shch.im * sc.re;
      where (negate) res.re = -res.re;
    }
    else {
      where (rcls == FP_ZERO) {      // Real part is 0.0.
        res.re = copysign(0.0, negate ? -1.0 : 1.0);
        res.im = NAN;           //    res.im = NAN + NAN;
        //#ifdef FE_INVALID
        //where (icls == FP_INFINITE) feraiseexcept(FE_INVALID);
        //#endif
      }
      else {
        res.re = NAN;
        res.im = NAN;
        //#ifdef FE_INVALID
        //feraiseexcept (FE_INVALID);
        //#endif
      }
    }
  }
  else where (rcls == FP_INFINITE) {       // Real part is infinite.
    where (icls == FP_ZERO) {          // Imaginary part is 0.0.
      res.re = negate ? -HUGE_VALD : HUGE_VALD;
      res.im = z.im;
    }
    else where (icls == FP_NORMAL) {          // Imaginary part is finite.
      sc = sincos(z.im);
      res.re = copysign(HUGE_VALD, sc.im);
      res.im = copysign(HUGE_VALD, sc.re);
      where (negate) res.re = -res.re;
    }
    else {       
      res.re = HUGE_VALD;
      res.im = NAN;          // res.im = NAN + NAN; The addition raises the invalid exception.
      //#ifdef FE_INVALID
      //where (icls == FP_INFINITE) feraiseexcept(FE_INVALID);
      //#endif
    }
  }
  else {
    res.re = NAN;
    res.im = z.im == 0.0 ? z.im : NAN;
  }
#endif
  return res;
}
#endif // csinh()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex ctanh(complex z) ;
#else
#if defined(_uses_ctanh_complex_h) || !defined(__cflow_processed) 
math_inline complex ctanh(complex z) 
{
  register complex res, sc, shch;
  register double den;

#if defined(__FAST_MATH__)
  sc = sincos(2.0 * z.im);
  shch = sinhcosh(2.0 * z.re);
  den = shch.im + sc.im;            //test: den = 1/(shch.im + sc.im); 
  res.re = shch.re / den;           //test: res.re = shch.re * den
  res.im = sc.re   / den;           //test: res.im = shch.re * den

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((!(rcls == FP_ZERO) && !(rcls == FP_NORMAL)) || (!(icls == FP_ZERO) && !(icls == FP_NORMAL))) {
  //  where (!isfinite(z.re) || !isfinite(z.im)) {
    where (rcls == FP_INFINITE) {
    //    where (isinf(z.re)) {
      res.re = copysign (1.0, z.re);
      res.im = copysign (0.0, z.im);
    }
    else where (icls == FP_ZERO) res = z; 
    //    else where (z.im == 0.0) res = z;
    else {
      res.re = NAN;
      res.im = NAN;
      //#ifdef FE_INVALID
      //where (isinf (z.im)) feraiseexcept (FE_INVALID);
      //#endif
    }
  }
  else {
    sc = sincos(2.0 * z.im);
    shch = sinhcosh(2.0 * z.re);
    den = shch.im + sc.im;
    res.re = shch.re / den;
    res.im = sc.re   / den;
  }
#endif
  return res;
}
#endif // ctanh()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex ccos(complex z) ;
#else
#if defined(_uses_ccos_complex_h) || !defined(__cflow_processed) 
math_inline complex ccos(complex z) 
{
/* old version:
  register complex y,res;
  y.re = -z.im;
  y.im =  z.re;
  res = ccosh(y); 
  return res;
*/
  register complex res, sc, shch;

#if defined(__FAST_MATH__)
  sc = sincos(z.re);
  shch = sinhcosh(z.im);
  res.re =  shch.im * sc.im;
  res.im = -shch.re * sc.re;      //res = ccosh(I*z);

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((!(rcls == FP_ZERO) && !(rcls == FP_NORMAL)) || (icls == FP_NAN)) {
  //  where (!isfinite(z.re) || isnan(z.im)) {
    where (z.re == 0.0 || z.im == 0.0) {
      res.re = NAN;
      res.im = 0.0;
      //#ifdef FE_INVALID
      //where (isinf(z.re)) feraiseexcept(FE_INVALID);
      //#endif
    }
    else where (icls == FP_INFINITE) {
    //    else where (isinf (z.im)) {
      res.re = HUGE_VALD;
      res.im = NAN;
      //#ifdef FE_INVALID
      //where (isinf(z.re)) feraiseexcept(FE_INVALID);
      //#endif
    }
    else {
      res.re = NAN;
      res.im = NAN;
      //#ifdef FE_INVALID
      //where (isfinite(z.im)) feraiseexcept(FE_INVALID);
      //#endif
    }
  }
  else {
    sc = sincos(z.re);
    shch = sinhcosh(z.im);
    res.re =  shch.im * sc.im;
    res.im = -shch.re * sc.re;
  }
#endif
  return res;
}
#endif // ccos()
#endif // Has Main

/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex csin(complex z) ;
#else
#if defined(_uses_csin_complex_h) || !defined(__cflow_processed) 
math_inline complex csin(complex z) 
{
  register complex res, sc, shch;

#if defined(__FAST_MATH__)
  sc = sincos(z.re);
  shch = sinhcosh(z.im);
  res.re = shch.im * sc.re;
  res.im = shch.re * sc.im;
  where (signbit(z.re)) res.re = -res.re;   //res = -I*csinh(I*z);

#else
  register int rcls, icls, negate;
  negate = signbit(z.re);
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  z.re = fabs(z.re);
  where ((icls == FP_ZERO) || (icls == FP_NORMAL)) {      // Imaginary part is finite. 
    where ((rcls == FP_ZERO) || (rcls == FP_NORMAL)) {          // Real part is finite.
      sc = sincos(z.re);
      shch = sinhcosh(z.im);
      res.re = shch.im * sc.re;
      res.im = shch.re * sc.im;
      where (negate)  res.re = -res.re;
    }
    else {
      where (icls == FP_ZERO) {         // Imaginary part is 0.0. 
        res.re = NAN;
        res.im = z.im;
        //#ifdef FE_INVALID
        //where (rcls == FP_INFINITE) feraiseexcept(FE_INVALID);
        //#endif
      }
      else {
        res.re = NAN;
        res.im = NAN;
        //#ifdef FE_INVALID
        //feraiseexcept(FE_INVALID);
        //#endif
      }
    }
  }
  else where (icls == FP_INFINITE) {      // Imaginary part is infinite.
    where (rcls == FP_ZERO) {      // Real part is 0.0.
      res.re = copysign(0.0, negate ? -1.0 : 1.0);
      res.im = z.im;
    }
    else where (rcls > FP_ZERO) {   // Real part is finite.
      sc = sincos(z.re);
      res.re = copysign(HUGE_VALD, sc.re);
      res.im = copysign(HUGE_VALD, sc.im);
      where (negate)         res.re = -res.re;
      where (signbit(z.im)) res.im = -res.im;
    }
    else {     
      res.re = NAN;
      res.im = HUGE_VALD;
      //#ifdef FE_INVALID
      //where (rcls == FP_INFINITE) feraiseexcept(FE_INVALID);
      //#endif
    }
  }
  else {
    where (rcls == FP_ZERO) res.re = copysign(0.0, negate ? -1.0 : 1.0);
    else res.re = NAN;
    res.im = NAN;
  }
#endif
  return res;
}
#endif // csin()
#endif // Has Main


/* ------------------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------------------- */

#if !defined(__HAS_MAIN) && !defined(_INLINE_MATH_)
extern math_inline complex ctan(complex z) ;
#else
#if defined(_uses_ctan_complex_h) || !defined(__cflow_processed) 
math_inline complex ctan(complex z) 
{
  register complex res, sc, shch;
  register double den;

#if defined(__FAST_MATH__)
  sc = sincos(2.0 * z.re);          // res = -I*ctanh(I*z);
  shch = sinhcosh(2.0 * z.im);
  den = sc.im + shch.im;          //test: den = 1/(sc.im + shch.im)
  res.re = sc.re   / den;         //test: res.re = sc.re   * den;
  res.im = shch.re / den;         //test: res.im = shch.re * den;

#else
  register int rcls, icls;
  rcls = fpclassify(z.re);
  icls = fpclassify(z.im);
  where ((!(rcls == FP_ZERO) && !(rcls == FP_NORMAL)) || (!(icls == FP_ZERO) && !(icls == FP_NORMAL))) {
    //  where (!isfinite(z.re) || !isfinite(z.im)) {
    where (icls == FP_INFINITE) {
      //    where (isinf(z.im)) {
      res.re = copysign(0.0, z.re);
      res.im = copysign(1.0, z.im);
    }
    else where (rcls == FP_ZERO) res = z;
    //    else where (z.re == 0.0) res = z;
    else {
      res.re = NAN;
      res.im = NAN;
      //#ifdef FE_INVALID
      //      where (isinf(z.re)) feraiseexcept(FE_INVALID);
      //#endif
    }
  }
  else {
    sc = sincos(2.0 * z.re);
    shch = sinhcosh(2.0 * z.im);
    den = sc.im + shch.im;  
    res.re = sc.re   / den; 
    res.im = shch.re / den;
  }
#endif
  return res;
}
#endif // ctan()
#endif // Has Main


#undef _math_setdouble
#undef HUGE_VALD
#undef NAN

#endif   /* ifndef _COMPLEX_H  */

---