// complexf.cpp  implementation of  complexf.h
//               fuller set of complex functions

#include "complexf.h"

static double PI = 3.1415926535897932384626433832795028841971 ;
static double PI_2 = PI/2.0;
static double LOG_TWO = 0.6931471805599453094172321214581765680755 ;
static complex<double> COMPLEX_I = complex<double>(0.0, 1.0);

complex<double>  cot(complex<double> &z)
{
  return cos(z)/sin(z);
}

complex<double>  acot(complex<double> &z)
{
  complex<double> r;
  
  r = COMPLEX_I * log((z-COMPLEX_I)/(z+COMPLEX_I)) / 2.0 ;
  if(real(r) < 0.0) r = PI + r;  // insure principal value
  return r;
}

complex<double>  coth(complex<double> &z)
{
  return cosh(z)/sinh(z);
}

complex<double>  asinh(complex<double> &z)
{
  complex<double> r;
  r = log( z + sqrt(1.0 + z*z));
  if(imag(r) > PI_2)          // insure principal value
    r = complex<double>(real(r), PI-imag(r));
  else if(imag(r) < -PI_2)
    r = complex<double>(real(r), -(PI+imag(r)));
  return r;
}

complex<double>  acosh(complex<double> &z)
{
  complex<double> r;
  
  r= 2.0 * log( sqrt((1.0+z)/2.0) + sqrt((z-1.0)/2.0));
  if(real(r) <= 0.0) r = -r;     // insure principal value
  return r;
}

complex<double>  atanh(complex<double> &z)
{
  return (log(1.0 + z) - log(1.0 - z))/2.0 ;
}

complex<double>  acoth(complex<double> &z)
{
  complex<double> r;
  
  r = log (( 1.0 + z ) / ( z - 1.0 )) / 2.0;
  if(imag(r) < 0.0) r=r+complex<double>(0.0,PI); // insure principal value
  return r;
}

complex<double>  acos(complex<double> &z)
{
  complex<double> r;
  r = - COMPLEX_I * log(z + COMPLEX_I*sqrt(1.0-z*z));
  if(real(r) < 0.0)       // insure principal result
    r = complex<double>(-real(r), imag(r));
  return r;
}

complex<double>  asin(complex<double> &z)
{
  complex<double> r;
  
  r = -COMPLEX_I * log( COMPLEX_I*z + sqrt(1.0-z*z));
  if(real(r) > PI_2)        // insure principal value
    r = complex<double>(PI-real(r), imag(r));
  else if(real(r) < -PI_2)
    r = complex<double>(-(PI+real(r)), imag(r));
  return r;
}

complex<double>  atan(complex<double> &z)
{
  return -COMPLEX_I * (log(1.0+COMPLEX_I*z)
         -log(1.0-COMPLEX_I*z)) / 2.0 ;
}


complex<double>  tan(complex<double> &z)
{
  return sin(z)/cos(z);
}

complex<double>  tanh(complex<double> &z)
{
  return sinh(z)/cosh(z);
}