exception handling complex elementary functions (fwd)

Andrew Binkley abinkley@cs.toronto.edu
Fri Jan 9 05:10:00 GMT 2004


Hi,
I've done as requested.   I decided not to use exact same format as
specfunc/error.h so the returned result is always similar to that of the
original functions.  Attached are complex\math.c (with some #ifndef
additions) and complex\math_eh.c.  Please let me know if you have any
other requests.

Thanks,
Andrew Binkley

On Wed, 7 Jan 2004, Brian Gough Wrote:
>Could you
>
> -- put the new functions in a separate file math_eh.c, so I don't need
> to put them directly in the existing file -- it is easier for
> maintenance that way.  Run the new file through indent to bring it
> into line with the GNU coding conventions.
>
> -- move the definitions in complex_math_eh.h into this file math_eh.c
> (there is no need for a separate header file in this case).
>
> -- It would be good practice to write these macros with arguments
> rather than using the various variable names implicitly in the
> definition.
>
> -- The macros are not exported so they do not need a GSL_COMPLEX
> prefix.  See the first definition in specfunc/error.h for an example
> of how it is done in the special functions.
>
> I'll deal with the configure script issues, such as testing whether
> the floating-point exception handling functions are available.
>
> Thanks,
>
> --
> Brian Gough
>
-------------- next part --------------
/* complex/math.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jorma Olavi T?htinen, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Basic complex arithmetic functions

 * Original version by Jorma Olavi T?htinen <jotahtin@cc.hut.fi>
 *
 * Modified for GSL by Brian Gough, 3/2000
 */

/* The following references describe the methods used in these
 * functions,
 *
 *   T. E. Hull and Thomas F. Fairgrieve and Ping Tak Peter Tang,
 *   "Implementing Complex Elementary Functions Using Exception
 *   Handling", ACM Transactions on Mathematical Software, Volume 20
 *   (1994), pp 215-244, Corrigenda, p553
 *
 *   Hull et al, "Implementing the complex arcsin and arccosine
 *   functions using exception handling", ACM Transactions on
 *   Mathematical Software, Volume 23 (1997) pp 299-335
 *
 *   Abramowitz and Stegun, Handbook of Mathematical Functions, "Inverse
 *   Circular Functions in Terms of Real and Imaginary Parts", Formulas
 *   4.4.37, 4.4.38, 4.4.39
 */

 /*  If COMPLEX_EXCEPTION_HANDLE is defined, alternative functions in
  *  complex/math_eh.c are used
  */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

/**********************************************************************
 * Complex numbers
 **********************************************************************/

#ifndef HIDE_INLINE_STATIC
gsl_complex
gsl_complex_rect (double x, double y)
{				/* return z = x + i y */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, x, y);
  return z;
}
#endif

gsl_complex
gsl_complex_polar (double r, double theta)
{				/* return z = r exp(i theta) */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, r * cos (theta), r * sin (theta));
  return z;
}

/**********************************************************************
 * Properties of complex numbers
 **********************************************************************/

double
gsl_complex_arg (gsl_complex z)
{				/* return arg(z),  -pi < arg(z) <= +pi */
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);

  if (x == 0.0 && y == 0.0)
    {
      return 0;
    }

  return atan2 (y, x);
}

#ifndef COMPLEX_EXCEPTION_HANDLE
double
gsl_complex_abs (gsl_complex z)
{				/* return |z| */
  return hypot (GSL_REAL (z), GSL_IMAG (z));
}
#endif

double
gsl_complex_abs2 (gsl_complex z)
{				/* return |z|^2 */
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);

  return (x * x + y * y);
}

double
gsl_complex_logabs (gsl_complex z)
{				/* return log|z| */
  double xabs = fabs (GSL_REAL (z));
  double yabs = fabs (GSL_IMAG (z));
  double max, u;

  if (xabs >= yabs)
    {
      max = xabs;
      u = yabs / xabs;
    }
  else
    {
      max = yabs;
      u = xabs / yabs;
    }

  /* Handle underflow when u is close to 0 */

  return log (max) + 0.5 * log1p (u * u);
}


/***********************************************************************
 * Complex arithmetic operators
 ***********************************************************************/

gsl_complex
gsl_complex_add (gsl_complex a, gsl_complex b)
{				/* z=a+b */
  double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  double br = GSL_REAL (b), bi = GSL_IMAG (b);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, ar + br, ai + bi);
  return z;
}

gsl_complex
gsl_complex_add_real (gsl_complex a, double x)
{				/* z=a+x */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_REAL (a) + x, GSL_IMAG (a));
  return z;
}

gsl_complex
gsl_complex_add_imag (gsl_complex a, double y)
{				/* z=a+iy */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_REAL (a), GSL_IMAG (a) + y);
  return z;
}


gsl_complex
gsl_complex_sub (gsl_complex a, gsl_complex b)
{				/* z=a-b */
  double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  double br = GSL_REAL (b), bi = GSL_IMAG (b);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, ar - br, ai - bi);
  return z;
}

gsl_complex
gsl_complex_sub_real (gsl_complex a, double x)
{				/* z=a-x */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_REAL (a) - x, GSL_IMAG (a));
  return z;
}

gsl_complex
gsl_complex_sub_imag (gsl_complex a, double y)
{				/* z=a-iy */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_REAL (a), GSL_IMAG (a) - y);
  return z;
}

gsl_complex
gsl_complex_mul (gsl_complex a, gsl_complex b)
{				/* z=a*b */
  double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  double br = GSL_REAL (b), bi = GSL_IMAG (b);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, ar * br - ai * bi, ar * bi + ai * br);
  return z;
}

gsl_complex
gsl_complex_mul_real (gsl_complex a, double x)
{				/* z=a*x */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, x * GSL_REAL (a), x * GSL_IMAG (a));
  return z;
}

gsl_complex
gsl_complex_mul_imag (gsl_complex a, double y)
{				/* z=a*iy */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, -y * GSL_IMAG (a), y * GSL_REAL (a));
  return z;
}

gsl_complex
gsl_complex_div (gsl_complex a, gsl_complex b)
{				/* z=a/b */
  double ar = GSL_REAL (a), ai = GSL_IMAG (a);
  double br = GSL_REAL (b), bi = GSL_IMAG (b);

  double s = 1.0 / gsl_complex_abs (b);

  double sbr = s * br;
  double sbi = s * bi;

  double zr = (ar * sbr + ai * sbi) * s;
  double zi = (ai * sbr - ar * sbi) * s;

  gsl_complex z;
  GSL_SET_COMPLEX (&z, zr, zi);
  return z;
}

gsl_complex
gsl_complex_div_real (gsl_complex a, double x)
{				/* z=a/x */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_REAL (a) / x, GSL_IMAG (a) / x);
  return z;
}

gsl_complex
gsl_complex_div_imag (gsl_complex a, double y)
{				/* z=a/(iy) */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_IMAG (a) / y, -GSL_REAL (a) / y);
  return z;
}

gsl_complex
gsl_complex_conjugate (gsl_complex a)
{				/* z=conj(a) */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_REAL (a), -GSL_IMAG (a));
  return z;
}

gsl_complex
gsl_complex_negative (gsl_complex a)
{				/* z=-a */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, -GSL_REAL (a), -GSL_IMAG (a));
  return z;
}

gsl_complex
gsl_complex_inverse (gsl_complex a)
{				/* z=1/a */
  double s = 1.0 / gsl_complex_abs (a);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, (GSL_REAL (a) * s) * s, -(GSL_IMAG (a) * s) * s);
  return z;
}

/**********************************************************************
 * Elementary complex functions
 **********************************************************************/

#ifndef COMPLEX_EXCEPTION_HANDLE
gsl_complex
gsl_complex_sqrt (gsl_complex a)
{				/* z=sqrt(a) */
  gsl_complex z;

  if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0)
    {
      GSL_SET_COMPLEX (&z, 0, 0);
    }
  else
    {
      double x = fabs (GSL_REAL (a));
      double y = fabs (GSL_IMAG (a));
      double w;

      if (x >= y)
	{
	  double t = y / x;
	  w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + t * t)));
	}
      else
	{
	  double t = x / y;
	  w = sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 + t * t)));
	}

      if (GSL_REAL (a) >= 0.0)
	{
	  double ai = GSL_IMAG (a);
	  GSL_SET_COMPLEX (&z, w, ai / (2.0 * w));
	}
      else
	{
	  double ai = GSL_IMAG (a);
	  double vi = (ai >= 0) ? w : -w;
	  GSL_SET_COMPLEX (&z, ai / (2.0 * vi), vi);
	}
    }

  return z;
}
#endif

gsl_complex
gsl_complex_sqrt_real (double x)
{				/* z=sqrt(x) */
  gsl_complex z;

  if (x >= 0)
    {
      GSL_SET_COMPLEX (&z, sqrt (x), 0.0);
    }
  else
    {
      GSL_SET_COMPLEX (&z, 0.0, sqrt (-x));
    }

  return z;
}

#ifndef COMPLEX_EXCEPTION_HANDLE
gsl_complex
gsl_complex_exp (gsl_complex a)
{				/* z=exp(a) */
  double rho = exp (GSL_REAL (a));
  double theta = GSL_IMAG (a);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, rho * cos (theta), rho * sin (theta));
  return z;
}
#endif

gsl_complex
gsl_complex_pow (gsl_complex a, gsl_complex b)
{				/* z=a^b */
  gsl_complex z;

  if (GSL_REAL (a) == 0 && GSL_IMAG (a) == 0.0)
    {
      GSL_SET_COMPLEX (&z, 0.0, 0.0);
    }
  else
    {
      double logr = gsl_complex_logabs (a);
      double theta = gsl_complex_arg (a);

      double br = GSL_REAL (b), bi = GSL_IMAG (b);

      double rho = exp (logr * br - bi * theta);
      double beta = theta * br + bi * logr;

      GSL_SET_COMPLEX (&z, rho * cos (beta), rho * sin (beta));
    }

  return z;
}

gsl_complex
gsl_complex_pow_real (gsl_complex a, double b)
{				/* z=a^b */
  gsl_complex z;

  if (GSL_REAL (a) == 0 && GSL_IMAG (a) == 0)
    {
      GSL_SET_COMPLEX (&z, 0, 0);
    }
  else
    {
      double logr = gsl_complex_logabs (a);
      double theta = gsl_complex_arg (a);
      double rho = exp (logr * b);
      double beta = theta * b;
      GSL_SET_COMPLEX (&z, rho * cos (beta), rho * sin (beta));
    }

  return z;
}

#ifndef COMPLEX_EXCEPTION_HANDLE
gsl_complex
gsl_complex_log (gsl_complex a)
{				/* z=log(a) */
  double logr = gsl_complex_logabs (a);
  double theta = gsl_complex_arg (a);

  gsl_complex z;
  GSL_SET_COMPLEX (&z, logr, theta);
  return z;
}
#endif

gsl_complex
gsl_complex_log10 (gsl_complex a)
{				/* z = log10(a) */
  return gsl_complex_mul_real (gsl_complex_log (a), 1 / log (10.));
}

gsl_complex
gsl_complex_log_b (gsl_complex a, gsl_complex b)
{
  return gsl_complex_div (gsl_complex_log (a), gsl_complex_log (b));
}

/***********************************************************************
 * Complex trigonometric functions
 ***********************************************************************/

#ifndef COMPLEX_EXCEPTION_HANDLE
gsl_complex
gsl_complex_sin (gsl_complex a)
{				/* z = sin(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);

  gsl_complex z;

  if (I == 0.0)
    {
      /* avoid returing negative zero (-0.0) for the imaginary part  */

      GSL_SET_COMPLEX (&z, sin (R), 0.0);
    }
  else
    {
      GSL_SET_COMPLEX (&z, sin (R) * cosh (I), cos (R) * sinh (I));
    }

  return z;
}
#endif

#ifndef COMPLEX_EXCEPTION_HANDLE
gsl_complex
gsl_complex_cos (gsl_complex a)
{				/* z = cos(a) */
  double R = GSL_REAL (a), I = GSL_IMAG (a);

  gsl_complex z;

  if (I == 0.0)
    {
      /* avoid returing negative zero (-0.0) for the imaginary part  */

      GSL_SET_COMPLEX (&z, cos (R), 0.0);
    }
  else
    {
      GSL_SET_COMPLEX (&z, cos (R) * cosh (I), sin (R) * sinh (-I));
    }

  return z;
}
#endif

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

  gsl_complex z;

  if (fabs (I) < 1)
    {
      double D = pow (cos (R), 2.0) + pow (sinh (I), 2.0);

      GSL_SET_COMPLEX (&z, 0.5 * sin (2 * R) / D, 0.5 * sinh (2 * I) / D);
    }
  else
    {
      double u = exp (-I);
      double C = 2 * u / (1 - pow (u, 2.0));
      double D = 1 + pow (cos (R), 2.0) * pow (C, 2.0);

      double S = pow (C, 2.0);
      double T = 1.0 / tanh (I);

      GSL_SET_COMPLEX (&z, 0.5 * sin (2 * R) * S / D, T / D);
    }

  return z;
}

gsl_complex
gsl_complex_sec (gsl_complex a)
{				/* z = sec(a) */
  gsl_complex z = gsl_complex_cos (a);
  return gsl_complex_inverse (z);
}

gsl_complex
gsl_complex_csc (gsl_complex a)
{				/* z = csc(a) */
  gsl_complex z = gsl_complex_sin (a);
  return gsl_complex_inverse (z);
}


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);
}


-------------- next part --------------

/* complex/math_eh.c
 * 
 * Copyright (C) 2003 Andrew Binkley
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Complex elementary functions using exception handling
 *
 * Original version by Andrew Binkley.  Take out the removebfmail
 * before sending <abinkley at removebfmail cs dot toronto dot edu>
 * Supervised by Professor Thomas Fairgrieve
 * Written for the GSL 6/2003
 */
/* The following reference describes the methods used in these
 * functions:
 *
 *   T. E. Hull and Thomas F. Fairgrieve and Ping Tak Peter Tang,
 *   "Implementing Complex Elementary Functions Using Exception
 *   Handling", ACM Transactions on Mathematical Software, Volume 20
 *   (1994), pp 215-244, Corrigenda, p553
 */

#ifdef COMPLEX_EXCEPTION_HANDLE	/*Used these methods if requested at compile time */

#include <math.h>
#include <float.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_math.h>	/*Used for macros like GSL_MAX */
#include <fenv.h>		/*Floating point enviroment */


#define CABS_COMPLEX_OVERFLOW( answer )                                      \
                          do{                                                \
                          feraiseexcept (FE_OVERFLOW);                       \
                          return answer;                                     \
			  } while (0)

#define COMPLEX_OVERFLOW( answer, answerreal, answerimag)                    \
                          do{                                                \
                          GSL_SET_COMPLEX (&answer, answerreal, answerimag); \
                          feraiseexcept (FE_OVERFLOW);                       \
                          return answer;                                     \
                          } while (0)

#define COMPLEX_UNDERFLOW( answer, answerreal, answerimag)                   \
                          do{                                                \
                          GSL_SET_COMPLEX (&answer, answerreal, answerimag); \
                          feraiseexcept (FE_UNDERFLOW);                      \
			  return answer;                                     \
			  } while (0)


#define COMPLEX_DOMAINERROR( answer, answerreal, answerimag)                 \
                          do{                                                \
                          GSL_SET_COMPLEX (&answer, answerreal, answerimag); \
                          feraiseexcept (FE_INVALID);                        \
			  return answer;                                     \
			  } while (0)

/****************************************************************************
 *Helper Functions
 ***************************************************************************/

/*int fhandle() 
 *This function is a close approximation to the enable, handle construct used by 
 *Hull, Fairgrieve and Tang.  By testing against values in fenv.h the type of
 *exception can be determined.
 *It is important to note that this function will only detect (and clear) overflows
 *and underflows.  NaN etc... are unaffected by this function
 */

int
fhandle ()
{
  int result = 0;		/*Changed if underflow and/or overflow */
  /*Get and clear status */
  result = fetestexcept (FE_UNDERFLOW | FE_OVERFLOW);
  feclearexcept (FE_UNDERFLOW | FE_OVERFLOW);
  return result;
}

/*********************************************************************************
 *Complex Elementary functions
 ********************************************************************************/

/* This is the cabs function that closely corresponds to Hull, Fairgrieve and
 * Tang.  First it tries to compute the result directly.  In the case of 
 * overflow or underflow,they are trapped and alternative methods executed.
 * Even then, an overflow is still possible.  This is handled with the standard 
 * gsl error handling mechanisms.
 */
double
gsl_complex_abs_eh (gsl_complex z)
{
/*Declarations and initializations*/
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);
  double answer, scaledx, scaledy, scaledanswer;

/*DBL_MAN_DIG is in float.h*/
  const int precision = DBL_MANT_DIG;
  int logbx, logby;

  /*Enable */
  answer = sqrt (x * x + y * y);
  /*Handle */
  if (fhandle ())		/*Overflow or underflow has occurred */
    {
      if (x == 0 || y == 0)
	{
	  answer = fabs (x) + fabs (y);
	}

      else
	{
	  logbx = ilogb (x);
	  logby = ilogb (y);

	  if (2 * fabs (logbx - logby) > precision + 1)
	    {
	      /*This case is invoked if exponents are so different
	       *that one may safely be ignored */
	      answer = GSL_MAX (fabs (x), fabs (y));
	    }

	  else			/*Scale so that fabs(x) is near 1 */
	    {
	      scaledx = scalbn (x, -logbx);
	      scaledy = scalbn (y, -logbx);
	      scaledanswer = sqrt ((scaledx) * (scaledx)
				   + (scaledy) * (scaledy));

	      /*The unscaled expression may overflow */
	      /*Enable */
	      answer = scalbn (scaledanswer, logbx);
	      /*Handle */
	      if (fhandle ())	/*must be overflow in scalb */
		{
		  CABS_COMPLEX_OVERFLOW (answer);
		}
	    }
	}
    }

  return answer;
}

/* This is the csqrt function that closely corresponds to Hull, Fairgrieve and 
 * Tang.  First it tries to compute the result directly.  In the case of 
 * overflow or underflow, they are trapped and alternative methods executed.  
 */
gsl_complex
gsl_complex_sqrt_eh (gsl_complex z)
{
/*Declarations and initializations*/
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);
  const double sqrt2 = M_SQRT2;
  double t, scaledx, scaledy, scaledt, temp;
  double answerreal, answerimag;

  const int precision = DBL_MANT_DIG;
  int logbx, logby, evennearlogbx;

  gsl_complex answer;

  /*Enable */
  t = sqrt (2.0 * (sqrt (x * x + y * y) + fabs (x)));

  if (x > 0)
    {
      answerreal = t / 2;
      answerimag = y / t;
    }

  else if (x < 0)
    {
      answerreal = fabs (y) / t;
      answerimag = copysign (t / 2, y);
      /*Returns t/2 with sign of y */
    }

  else				/*This is the x = 0 case */
    {
      temp = sqrt (fabs (y)) / sqrt2;
      answerreal = temp;
      answerimag = copysign (temp, y);
    }

  /*Handle */
  if (fhandle ())
    {
      if (x == 0)
	{
	  temp = sqrt (fabs (y)) / sqrt2;
	  answerreal = temp;
	  answerimag = copysign (temp, y);
	}

      else if (y == 0)		/*Trivial case */
	{
	  if (x > 0)
	    {
	      answerreal = sqrt (x);
	      answerimag = 0;
	    }

	  else
	    {
/*Declarations and initializations*/
	      answerreal = 0;
	      answerimag = sqrt (-x);
	    }
	}

      else			/*Determine t */
	{
	  logbx = ilogb (x);
	  logby = ilogb (y);

	  if (logby - logbx > precision)	/*x may be ignored */
	    {
	      t = sqrt2 * sqrt (fabs (y));
	    }

	  else if ((2 * (logbx - logby)) > (precision + 1))
	    {			/*y may be ignored */
	      t = 2 * sqrt (fabs (x));
	    }

	  else			/*Scale and unscale so that abs(x) is near 1 with */
	    {			/*even exponent */
	      evennearlogbx = (logbx + logbx % 2);
	      scaledx = scalbn (x, -evennearlogbx);
	      scaledy = scalbn (y, -evennearlogbx);
	      scaledt = sqrt (2 *
			      (sqrt
			       ((scaledx) * (scaledx) +
				(scaledy) * (scaledy)) + fabs (scaledx)));
	      t = scalbn (scaledt, evennearlogbx / 2);
	    }

	  /*t has now been determined, must generate approp. answer using t */
	  if (x > 0)
	    {
	      answerreal = t / 2;

	      /*The original pseudo code for this function calls for an
	       *additional enable-handle construct to handle underflow of 
	       *answerimag here. Instead we allow a gradual underflow to zero.
	       */
	      answerimag = y / t;
	    }

	  else			/*(x < 0) */
	    {
	      /*The original pseudo code for this function calls for an
	       *additional enable-handle construct to handle underflow of
	       *answerreal here.  Instead we allow a gradual underflow to zero.
	       */
	      answerreal = fabs (y) / t;
	      answerimag = copysign (t / 2, y);
	    }
	}
    }

  GSL_SET_COMPLEX (&answer, answerreal, answerimag);
  return answer;
}

/* This is the cexp function that closely corresponds to Hull, Fairgrieve and
 * Tang.  First it tries to compute the result directly.  In the case of 
 * overflow or underflow,they are trapped and alternative methods executed.
 */
gsl_complex
gsl_complex_exp_eh (gsl_complex z)
{
/*Declarations and initializations*/
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);
  double expx;
  double answerreal, answerimag;

  gsl_complex answer;

  /*Enable */
  expx = exp (x);
  answerreal = expx * cos (y);
  answerimag = expx * sin (y);

  /*Handle */
  if (fhandle ())		/*Underflow or overflow has occurred */
    {
      /*Enable */
      expx = exp (x);		/*Test if range error in expx */
      /*Handle */
      if (fhandle ())		/*If error is exp, unavoidable */
	{
	  if (x > 0)
	    {
	      COMPLEX_OVERFLOW (answer, answerreal, answerimag);
	    }

	  else			/*Allow gradual underflow */
	    {
	      COMPLEX_UNDERFLOW (answer, answerreal, answerimag);
	    }
	}
      else			/*The underflow is not in exp */
	{
	  /*The original pseudo code for this function calls for several
	   *additional enable-handle constructs to handle underflow of
	   *answerreal and/or answerimag here.  Instead we allow a gradual
	   *underflow to zero, and reset the underflow flag.
	   */

	  COMPLEX_UNDERFLOW (answer, answerreal, answerimag);
	}
    }

  GSL_SET_COMPLEX (&answer, answerreal, answerimag);
  return answer;
}

/* This is the clog function that closely corresponds to Hull, Fairgrieve and
 * Tang.  Three special cases are considered first.  The real part is then 
 * calculated.  The calculation carried out is determined by testing whether max 
 * (|x| , |y|) is outside the interval (1/2, 2^0.5) or not.  In the case that 
 * log1p must be used it is preferable to have doubled (i.e. long double) precision
 * available.  Since the precision of long double isn't standardized the error
 * bounds are invalid on certain machines.
 */
gsl_complex
gsl_complex_log_eh (gsl_complex z)
{
/*Declarations and initializations*/
  const double sqrt2 = M_SQRT2;
  const double logradix = 1.0 / M_LOG2E;
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);
  double M, m, scaledM, scaledm, scaledr;
  double answerreal, answerimag;

  const int precision = DBL_MANT_DIG;
  int scale;

  gsl_complex answer;

  M = GSL_MAX (fabs (x), fabs (y));
  m = GSL_MIN (fabs (x), fabs (y));

  if (M == 0)			/*Log isn't defined in this case */
    {
      COMPLEX_DOMAINERROR (answer, answerreal, answerimag);
    }

  /*Enable */
  answerimag = gsl_complex_arg (z);	/*Approximation to arg(z) */
  /*Handle */
  if (fhandle ())		/*Underflow of arg has occurred */
    {				/*It can be allowed to be gradual */
      answerreal = log (M);
      COMPLEX_UNDERFLOW (answer, answerreal, answerimag);
    }
/*Now determine real part*/
  if (m == 0)
    {
      answerreal = log (M);
    }

  else if (M <= 0.5 || M >= sqrt2)
    {
      /*Enable */
      answerreal = 0.5 * log (M * M + m * m);
      /*Handle */
      if (fhandle ())		/*Must be overflow in M^2 or */
	{			/*underflow in m^2 */

	  /* Some log() functions set FE_DIVBYZERO in error when 
	     both M^2 and m^2 have underflowed.  Reset this flag. */
	  feclearexcept (FE_DIVBYZERO);

	  if ((2 * (ilogb (M) - ilogb (m))) > (precision + 1))
	    {
	      /*m may be safely ignored */
	      answerreal = log (M);
	    }

	  else			/*Scale to avoid exceptions and keep components */
	    {			/*of real answer with the same sign */
	      /*Find the approp. scale value */
	      if (M > 1)	/*Must have been an overflow */
		{
		  scale = ilogb (M);
		}

	      else		/*Must have been an underflow */
		{
		  scale = ilogb (M) + 2;
		}

	      /*Use the scaled value */
	      scaledM = scalbn (M, -scale);
	      scaledm = scalbn (m, -scale);
	      scaledr = scaledM * scaledM + scaledm * scaledm;
	      answerreal = scale * logradix + 0.5 * log (scaledr);
	    }
	}
    }
  else				/* 1/2 < M < sqrt2 */
    {
      /*Enable */
      /*If in this block we will want to use long double precision
       *to evaluate the argument to log1p
       */
      /*Block Declarations and initializations */
      long double ldm = m;
      long double ldM = M;
      double arglog1p;

      arglog1p = (ldM - 1) * (ldM + 1) + ldm * ldm;
      answerreal = 0.5 * log1p (arglog1p);

      /*Handle */
      if (fhandle ())		/*Must be underflow in m*m */
	{
	  answerreal = log (M);
	}
    }

  GSL_SET_COMPLEX (&answer, answerreal, answerimag);
  return answer;
}

/* This is the csin function that closely corresponds to Hull, Fairgrieve and 
 * Tang.  It is fairly straightforward, however there are some fringe overflows
 * that could be avoided
 */

gsl_complex
gsl_complex_sin_eh (gsl_complex z)
{
/*Declare*/
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);
  double coshy, sinhy;
  double answerreal, answerimag;
  gsl_complex answer;

  /*Enable */
  answerreal = sin (x) * cosh (y);
  answerimag = cos (x) * sinh (y);

  /*Handle */
  if (fhandle ())		/*Underflow or overflow */
    {
      /*Enable */
      coshy = cosh (y);
      sinhy = sinh (y);
      /*Handle */
      if (fhandle ())		/*Must be overflow */
	{
	  COMPLEX_OVERFLOW (answer, answerreal, answerimag);
	}
      else			/*Must be underflow */
	{
	  /*The original pseudo code for this function calls for an
	   *additional enable-handle construct to handle underflow here. 
	   *Instead we allow a gradual underflow to zero.
	   */
	  COMPLEX_UNDERFLOW (answer, answerreal, answerimag);
	}
    }

  GSL_SET_COMPLEX (&answer, answerreal, answerimag);
  return answer;
}


/* This is the ccos function that closely corresponds to Hull, Fairgrieve and
 * Tang.  It is fairly straightforward, however there are some fringe 
 * overflows that could be avoided.
 */
gsl_complex
gsl_complex_cos_eh (gsl_complex z)
{
/*Declare*/
  double x = GSL_REAL (z);
  double y = GSL_IMAG (z);
  double coshy, sinhy;
  double answerreal, answerimag;
  gsl_complex answer;

  /*Enable */
  answerreal = cos (x) * cosh (y);
  answerimag = -sin (x) * sinh (y);
  /*Handle */
  if (fhandle ())		/*Overflow or underflow */
    {
      /*Enable */
      coshy = cosh (y);
      sinhy = sinh (y);
      /*Handle */
      if (fhandle ())		/*Must be overflow */
	{
	  COMPLEX_OVERFLOW (answer, answerreal, answerimag);
	}
      else			/*Must be underflow */
	{
	  /*The original pseudo code for this function calls for an
	   *additional enable-handle construct to handle underflow here. 
	   *Instead we allow a gradual underflow to zero.
	   */
	  COMPLEX_UNDERFLOW (answer, answerreal, answerimag);
	}
    }

  GSL_SET_COMPLEX (&answer, answerreal, answerimag);
  return answer;
}

#endif	/*End of conditional #ifdef COMPLEX_EXCEPTION_HANDLE*/


More information about the Gsl-discuss mailing list