/* 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 * 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 #include #include #include #include /*Used for macros like GSL_MAX*/ #include /*Floating point enviroment*/ #include "complex_math_eh.h" /*This header uses flags*/ /*#include "complex_math_eh_a.h"*//*This alternative header uses GSL*/ /*error handling macros*/ /**************************************************************************** *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_GSL_COMPLEX_OVERFLOW; } } } } 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) { 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; } /* 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 */ { 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; } /* 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*/ { 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; } /* 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*/ { 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; }