exception handling complex elementary functions (fwd)

Andrew Binkley abinkley@cs.toronto.edu
Tue Jan 6 04:12:00 GMT 2004



---------- Forwarded message ----------
Date: Mon, 5 Jan 2004 22:56:07 -0500
From: Andrew Binkley <andrew.binkley@utoronto.ca>
To: abinkley@cs.toronto.edu
Subject: exception handling complex elementary functions

Hello,
I want to get some complex elementary functions that use exception handling
included in the main gsl release.  These functions have the advantage of
being within strict error bounds.  One file below is to replace
complex\math.c.  Another is to be placed in the gsl include directory
(complex_math_eh.h).  The third is some documentation.

The strict error bound functions are used if COMPLEX_EXCEPTION_HANDLE is
defined (see math.c).  I assume that this would be done at compile time, but
am not sure how.  The user can still use the existing gsl complex elementary
functions if COMPLEX_EXCEPTION_HANDLE is not defined.

I have compiled this addition with the gsl and was successful.  I have also
tested them for error bounds in a separate set-up.

Any pointers/suggestions about what I can do to get this included in the
main release would be appreciated.  Or if it looks ready, could you please
check it in.

Thanks,
Andrew Binkley

-------------- 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
 */

 /*  1/2004
  *  Added strict error bound functions if COMPLEX_EXCEPTION_HANDLE is 
  *  defined.  Otherwise the original methods are used.
  */

 /* 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
  *
  */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_errno.h>
#ifdef COMPLEX_EXCEPTION_HANDLE	/*Additional includes for strict*/
#include <float.h>		/*error bound methods*/
#include <fenv.h>               /*Floating point enviroment*/
#include <gsl/complex_math_eh.h>/*This header specifies flags*/
#endif				/*to be used for over+undeflow*/


#ifdef COMPLEX_EXCEPTION_HANDLE	/*Additional function for strict*/
/****************************************************************************
 *Helper Function for strict error bound methods
 ***************************************************************************/

/*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;
}
#endif

/**********************************************************************
 * 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

#ifdef COMPLEX_EXCEPTION_HANDLE
/* 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 (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_GSL_COMPLEX_OVERFLOW;
                }
            }
        }
    }

  return answer;
}
#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

#ifdef COMPLEX_EXCEPTION_HANDLE
/* 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 (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;
}
#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;
}

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;
}
#endif

#ifdef COMPLEX_EXCEPTION_HANDLE
/* 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 (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)
            {
              GSL_COMPLEX_OVERFLOW;
            }

          else                   /*Allow gradual underflow*/
            {
              GSL_COMPLEX_UNDERFLOW;
            }
         }
       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.
           */

              GSL_COMPLEX_UNDERFLOW;
         }
     }

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

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

#ifdef COMPLEX_EXCEPTION_HANDLE

/* 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 (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 */
    {
      GSL_COMPLEX_DOMAINERROR;
    }

  /*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);
      GSL_COMPLEX_UNDERFLOW;  
    }
/*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;
}
#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

#ifdef COMPLEX_EXCEPTION_HANDLE
/* 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 (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*/
        {
          GSL_COMPLEX_OVERFLOW;
        }
      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.
           */
          GSL_COMPLEX_UNDERFLOW;
        }
    }

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

#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

#ifdef COMPLEX_EXCEPTION_HANDLE
/* 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 (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*/
        {
          GSL_COMPLEX_OVERFLOW;
        }
      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.
           */
          GSL_COMPLEX_UNDERFLOW;
        }
    }

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

#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/complex_math_eh.h
 * 
 * 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.
 */

/*
 *Define error handling macros.
 *These macros are invoked when we cannot avoid an error.  Because
 *of different return types cabs has its own macro. These macros
 *set floating point register flags as opposed to aborting the 
 *functions for overflow and underflow.  DOMAINERROR still aborts.
 */

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

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

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


#define GSL_COMPLEX_DOMAINERROR                                              \
                          do{                                                \
                          GSL_SET_COMPLEX (&answer, answerreal, answerimag); \
                          feraiseexcept (FE_INVALID);                        \
			  return answer;                                     \
			} while (0)
-------------- next part --------------
Complex Elementary Functions Using Exception Handling


    All these functions are alternative functions that are only used if 
COMPLEX_EXCEPTION_HANDLE is defined at compile time.  These functions 
are offered as an alternative to other GSL complex elementary functions 
since while slightly slower to execute they have better error bounds.  E 
is the machine epsilon.  For more details on the error bounds see: 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.


.......


Properties of complex numbers


Function: double gsl_complex_abs (gsl_complex z)

    This function returns the magnitude of the complex number z, |z|.  
Any overflowed components are set to inf.  The appropriate 
floating-point register flags are set.  Theoretical error bound: E + 
E(sqrt).


.......


Elementary Complex Functions


Function: gsl_complex gsl_complex_sqrt (gsl_complex z)

    This function returns the square root of the complex number z, \sqrt 
z.  Any underflowed components are allowed to gradually underflow.  The 
appropriate floating point register flags are set.  There is a branch 
cut along the negative real axis from -oo to 0.  Theoretical error 
bound: sqrt( 2.5E^2 + 4.5E*E(sqrt) + 2.25E(sqrt)^2 ).


Function: gsl_complex gsl_complex_exp (gsl_complex z)

    This function returns the complex exponential of the complex number 
z, \exp(z).  Any overflowed components are set to inf.  Any underflowed 
components are allowed to gradually underflow.  The appropriate floating 
point register flags are set.  Theoretical error bound: sqrt( E^2 + ( E 
+ E(exp) + E(sincos) )^2 ).


Function: gsl_complex gsl_complex_log (gsl_complex z)

    This function returns the complex natural logarithm (base e) of the 
complex number z, \log(z).  Any underflowed components are allowed to 
gradually underflow.  The appropriate floating point register flags are 
set.  There is a branch cut along the negative real axis from -oo to 0.  
Depending on the precision of long double the theoretical error bound is 
between max( 2.886E + E(log), 6.495E + E(log1p), 4.457E + E(log1p) 
+E(arg) ) and max( 2.886E + E(log), 2.165E + E(log1p), E(arg) ).

.......


Complex Trigonometric Functions


Function: gsl_complex gsl_complex_sin (gsl_complex z)

    This function returns the complex sine of the complex number z, 
\sin(z).  Any overflowed components are set to inf.  Any underflowed 
components are allowed to gradually underflow.  The appropriate floating 
point register flags are set.  Theoretical error bound: E + E(sincos) + 
E(sinhcosh).


Function: gsl_complex gsl_complex_cos (gsl_complex z)

    This function returns the complex cosine of the complex number z, 
\cos(z).  Any overflowed components are set to inf.  Any underflowed 
components are allowed to gradually underflow.  The appropriate floating 
point register flags are set.  Theoretical error bound: E + E(sincos) + 
E(sinhcosh).






More information about the Gsl-discuss mailing list