math.c

一个复数运算的例子

}


gsl_complex
gsl_complex_cot (gsl_complex a)
{				/* z = cot(a) */
  gsl_complex z = gsl_complex_tan (a);
  return gsl_complex_inverse (z);
}

/**********************************************************************
 * Inverse Complex Trigonometric Functions
 **********************************************************************/

gsl_complex
gsl_complex_arcsin (gsl_complex a)
{				/* z = arcsin(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);
  gsl_complex z;

  if (I == 0)
    {
      z = gsl_complex_arcsin_real (R);
    }
  else
    {
      double x = fabs (R), y = fabs (I);
      double r = hypot (x + 1, y), s = hypot (x - 1, y);
      double A = 0.5 * (r + s);
      double B = x / A;
      double y2 = y * y;

      double real, imag;

      const double A_crossover = 1.5, B_crossover = 0.6417;

      if (B <= B_crossover)
	{
	  real = asin (B);
	}
      else
	{
	  if (x <= 1)
	    {
	      double D = 0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));
	      real = atan (x / sqrt (D));
	    }
	  else
	    {
	      double Apx = A + x;
	      double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));
	      real = atan (x / (y * sqrt (D)));
	    }
	}

      if (A <= A_crossover)
	{
	  double Am1;

	  if (x < 1)
	    {
	      Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));
	    }
	  else
	    {
	      Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));
	    }

	  imag = log1p (Am1 + sqrt (Am1 * (A + 1)));
	}
      else
	{
	  imag = log (A + sqrt (A * A - 1));
	}

      GSL_SET_COMPLEX (&z, (R >= 0) ? real : -real, (I >= 0) ? imag : -imag);
    }

  return z;
}

gsl_complex
gsl_complex_arcsin_real (double a)
{				/* z = arcsin(a) */
  gsl_complex z;

  if (fabs (a) <= 1.0)
    {
      GSL_SET_COMPLEX (&z, asin (a), 0.0);
    }
  else
    {
      if (a < 0.0)
	{
	  GSL_SET_COMPLEX (&z, -M_PI_2, acosh (-a));
	}
      else
	{
	  GSL_SET_COMPLEX (&z, M_PI_2, -acosh (a));
	}
    }

  return z;
}

gsl_complex
gsl_complex_arccos (gsl_complex a)
{				/* z = arccos(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);
  gsl_complex z;

  if (I == 0)
    {
      z = gsl_complex_arccos_real (R);
    }
  else
    {
      double x = fabs (R), y = fabs (I);
      double r = hypot (x + 1, y), s = hypot (x - 1, y);
      double A = 0.5 * (r + s);
      double B = x / A;
      double y2 = y * y;

      double real, imag;

      const double A_crossover = 1.5, B_crossover = 0.6417;

      if (B <= B_crossover)
	{
	  real = acos (B);
	}
      else
	{
	  if (x <= 1)
	    {
	      double D = 0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));
	      real = atan (sqrt (D) / x);
	    }
	  else
	    {
	      double Apx = A + x;
	      double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));
	      real = atan ((y * sqrt (D)) / x);
	    }
	}

      if (A <= A_crossover)
	{
	  double Am1;

	  if (x < 1)
	    {
	      Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));
	    }
	  else
	    {
	      Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));
	    }

	  imag = log1p (Am1 + sqrt (Am1 * (A + 1)));
	}
      else
	{
	  imag = log (A + sqrt (A * A - 1));
	}

      GSL_SET_COMPLEX (&z, (R >= 0) ? real : M_PI - real, (I >= 0) ? -imag : imag);
    }

  return z;
}

gsl_complex
gsl_complex_arccos_real (double a)
{				/* z = arccos(a) */
  gsl_complex z;

  if (fabs (a) <= 1.0)
    {
      GSL_SET_COMPLEX (&z, acos (a), 0);
    }
  else
    {
      if (a < 0.0)
	{
	  GSL_SET_COMPLEX (&z, M_PI, -acosh (-a));
	}
      else
	{
	  GSL_SET_COMPLEX (&z, 0, acosh (a));
	}
    }

  return z;
}

gsl_complex
gsl_complex_arctan (gsl_complex a)
{				/* z = arctan(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);
  gsl_complex z;

  if (I == 0)
    {
      GSL_SET_COMPLEX (&z, atan (R), 0);
    }
  else
    {
      /* FIXME: This is a naive implementation which does not fully
         take into account cancellation errors, overflow, underflow
         etc.  It would benefit from the Hull et al treatment. */

      double r = hypot (R, I);

      double imag;

      double u = 2 * I / (1 + r * r);

      /* FIXME: the following cross-over should be optimized but 0.1
         seems to work ok */

      if (fabs (u) < 0.1)
	{
	  imag = 0.25 * (log1p (u) - log1p (-u));
	}
      else
	{
	  double A = hypot (R, I + 1);
	  double B = hypot (R, I - 1);
	  imag = 0.5 * log (A / B);
	}

      if (R == 0)
	{
	  if (I > 1)
	    {
	      GSL_SET_COMPLEX (&z, M_PI_2, imag);
	    }
	  else if (I < -1)
	    {
	      GSL_SET_COMPLEX (&z, -M_PI_2, imag);
	    }
	  else
	    {
	      GSL_SET_COMPLEX (&z, 0, imag);
	    };
	}
      else
	{
	  GSL_SET_COMPLEX (&z, 0.5 * atan2 (2 * R, ((1 + r) * (1 - r))), imag);
	}
    }

  return z;
}

gsl_complex
gsl_complex_arcsec (gsl_complex a)
{				/* z = arcsec(a) */
  gsl_complex z = gsl_complex_inverse (a);
  return gsl_complex_arccos (z);
}

gsl_complex
gsl_complex_arcsec_real (double a)
{				/* z = arcsec(a) */
  gsl_complex z;

  if (a <= -1.0 || a >= 1.0)
    {
      GSL_SET_COMPLEX (&z, acos (1 / a), 0.0);
    }
  else
    {
      if (a >= 0.0)
	{
	  GSL_SET_COMPLEX (&z, 0, acosh (1 / a));
	}
      else
	{
	  GSL_SET_COMPLEX (&z, M_PI, -acosh (-1 / a));
	}
    }

  return z;
}

gsl_complex
gsl_complex_arccsc (gsl_complex a)
{				/* z = arccsc(a) */
  gsl_complex z = gsl_complex_inverse (a);
  return gsl_complex_arcsin (z);
}

gsl_complex
gsl_complex_arccsc_real (double a)
{				/* z = arccsc(a) */
  gsl_complex z;

  if (a <= -1.0 || a >= 1.0)
    {
      GSL_SET_COMPLEX (&z, asin (1 / a), 0.0);
    }
  else
    {
      if (a >= 0.0)
	{
	  GSL_SET_COMPLEX (&z, M_PI_2, -acosh (1 / a));
	}
      else
	{
	  GSL_SET_COMPLEX (&z, -M_PI_2, -acosh (-1 / a));
	}
    }

  return z;
}

gsl_complex
gsl_complex_arccot (gsl_complex a)
{				/* z = arccot(a) */
  gsl_complex z;

  if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0)
    {
      GSL_SET_COMPLEX (&z, M_PI_2, 0);
    }
  else
    {
      z = gsl_complex_inverse (a);
      z = gsl_complex_arctan (z);
    }

  return z;
}

/**********************************************************************
 * Complex Hyperbolic Functions
 **********************************************************************/

gsl_complex
gsl_complex_sinh (gsl_complex a)
{				/* z = sinh(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, sinh (R) * cos (I), cosh (R) * sin (I));
  return z;
}

gsl_complex
gsl_complex_cosh (gsl_complex a)
{				/* z = cosh(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, cosh (R) * cos (I), sinh (R) * sin (I));
  return z;
}

gsl_complex
gsl_complex_tanh (gsl_complex a)
{				/* z = tanh(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);

  gsl_complex z;

  if (fabs(R) < 1.0) 
    {
      double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);
      
      GSL_SET_COMPLEX (&z, sinh (R) * cosh (R) / D, 0.5 * sin (2 * I) / D);
    }
  else
    {
      double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);
      double F = 1 + pow (cos (I) / sinh (R), 2.0);

      GSL_SET_COMPLEX (&z, 1.0 / (tanh (R) * F), 0.5 * sin (2 * I) / D);
    }

  return z;
}

gsl_complex
gsl_complex_sech (gsl_complex a)
{				/* z = sech(a) */
  gsl_complex z = gsl_complex_cosh (a);
  return gsl_complex_inverse (z);
}

gsl_complex
gsl_complex_csch (gsl_complex a)
{				/* z = csch(a) */
  gsl_complex z = gsl_complex_sinh (a);
  return gsl_complex_inverse (z);
}

gsl_complex
gsl_complex_coth (gsl_complex a)
{				/* z = coth(a) */
  gsl_complex z = gsl_complex_tanh (a);
  return gsl_complex_inverse (z);
}

/**********************************************************************
 * Inverse Complex Hyperbolic Functions
 **********************************************************************/

gsl_complex
gsl_complex_arcsinh (gsl_complex a)
{				/* z = arcsinh(a) */
  gsl_complex z = gsl_complex_mul_imag(a, 1.0);
  z = gsl_complex_arcsin (z);
  z = gsl_complex_mul_imag (z, -1.0);
  return z;
}

gsl_complex
gsl_complex_arccosh (gsl_complex a)
{				/* z = arccosh(a) */
  gsl_complex z = gsl_complex_arccos (a);
  z = gsl_complex_mul_imag (z, GSL_IMAG(z) > 0 ? -1.0 : 1.0);
  return z;
}

gsl_complex
gsl_complex_arccosh_real (double a)
{				/* z = arccosh(a) */
  gsl_complex z;

  if (a >= 1)
    {
      GSL_SET_COMPLEX (&z, acosh (a), 0);
    }
  else
    {
      if (a >= -1.0)
	{
	  GSL_SET_COMPLEX (&z, 0, acos (a));
	}
      else
	{
	  GSL_SET_COMPLEX (&z, acosh (-a), M_PI);
	}
    }

  return z;
}

gsl_complex
gsl_complex_arctanh (gsl_complex a)
{				/* z = arctanh(a) */
  if (GSL_IMAG (a) == 0.0)
    {
      return gsl_complex_arctanh_real (GSL_REAL (a));
    }
  else
    {
      gsl_complex z = gsl_complex_mul_imag(a, 1.0);
      z = gsl_complex_arctan (z);
      z = gsl_complex_mul_imag (z, -1.0);
      return z;
    }
}

gsl_complex
gsl_complex_arctanh_real (double a)
{				/* z = arctanh(a) */
  gsl_complex z;

  if (a > -1.0 && a < 1.0)
    {
      GSL_SET_COMPLEX (&z, atanh (a), 0);
    }
  else
    {
      GSL_SET_COMPLEX (&z, atanh (1 / a), (a < 0) ? M_PI_2 : -M_PI_2);
    }

  return z;
}

gsl_complex
gsl_complex_arcsech (gsl_complex a)
{				/* z = arcsech(a); */
  gsl_complex t = gsl_complex_inverse (a);
  return gsl_complex_arccosh (t);
}

gsl_complex
gsl_complex_arccsch (gsl_complex a)
{				/* z = arccsch(a) */
  gsl_complex t = gsl_complex_inverse (a);
  return gsl_complex_arcsinh (t);
}

gsl_complex
gsl_complex_arccoth (gsl_complex a)
{				/* z = arccoth(a) */
  gsl_complex t = gsl_complex_inverse (a);
  return gsl_complex_arctanh (t);
}