ESYS13  Revision_
UnaryFuncs.h
Go to the documentation of this file.
00001 
00002 /*******************************************************
00003 *
00004 * Copyright (c) 2003-2012 by University of Queensland
00005 * Earth Systems Science Computational Center (ESSCC)
00006 * http://www.uq.edu.au/esscc
00007 *
00008 * Primary Business: Queensland, Australia
00009 * Licensed under the Open Software License version 3.0
00010 * http://www.opensource.org/licenses/osl-3.0.php
00011 *
00012 *******************************************************/
00013 
00014 
00015 #if !defined escript_UnaryFuncs_20041124_H
00016 #define escript_UnaryFuncs_20041124_H
00017 #include "system_dep.h"
00018 
00019 namespace escript {
00020 
00021 #ifndef FP_NAN
00022 #define FP_NAN IEEE_NaN()
00023 #endif
00024 
00025 #ifndef INFINITY
00026 #define INFINITY IEEE_Infy()
00027 #endif
00028 
00029 //======================================================================
00030 
00031 inline
00032 double log1p (const double x)
00033 {
00034   volatile double y;
00035   y = 1 + x;
00036   return log(y) - ((y-1)-x)/y ;
00037 }
00038 
00039 //======================================================================
00040 
00041 inline
00042 float IEEE_NaN()
00043 {
00044    static unsigned char nan[4]={ 0, 0, 0xc0, 0x7f };
00045    return *( float *)nan;
00046 }
00047 
00048 //======================================================================
00049 
00050 inline
00051 double IEEE_Infy()
00052 {
00053    static unsigned char infy[8]={ 0, 0, 0, 0, 0, 0, 0xf0, 0x7f };
00054    return *( double *)infy;
00055 }
00056 
00057 
00058 //======================================================================
00059 
00060 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
00061 inline
00062 double
00063 acosh_substitute (const double x)
00064 {
00065   if (x > 1.0 / SQRT_DBL_EPSILON)
00066     {
00067       return log (x) + M_LN2;
00068     }
00069   else if (x > 2)
00070     {
00071       return log (2 * x - 1 / (sqrt (x * x - 1) + x));
00072     }
00073   else if (x > 1)
00074     {
00075       double t = x - 1;
00076       return log1p (t + sqrt (2 * t + t * t));
00077     }
00078   else if (x == 1)
00079     {
00080       return 0;
00081     }
00082   else
00083     {
00084       return FP_NAN;
00085     }
00086 }
00087 
00088 
00089 //======================================================================
00090 
00091 inline
00092 double
00093 asinh_substitute (const double x)
00094 {
00095   double a = fabs (x);
00096   double s = (x < 0) ? -1 : 1;
00097 
00098   if (a > 1 / SQRT_DBL_EPSILON)
00099     {
00100       return s * (log (a) + M_LN2);
00101     }
00102   else if (a > 2)
00103     {
00104       return s * log (2 * a + 1 / (a + sqrt (a * a + 1)));
00105     }
00106   else if (a > SQRT_DBL_EPSILON)
00107     {
00108       double a2 = a * a;
00109       return s * log1p (a + a2 / (1 + sqrt (1 + a2)));
00110     }
00111   else
00112     {
00113       return x;
00114     }
00115 }
00116 
00117 
00118 //======================================================================
00119 
00120 inline
00121 double
00122 atanh_substitute (const double x)
00123 {
00124   double a = fabs (x);
00125   double s = (x < 0) ? -1 : 1;
00126 
00127   if (a > 1)
00128     {
00129       return FP_NAN;
00130     }
00131   else if (a == 1)
00132     {
00133       return (x < 0) ? -INFINITY : INFINITY;
00134     }
00135   else if (a >= 0.5)
00136     {
00137       return s * 0.5 * log1p (2 * a / (1 - a));
00138     }
00139   else if (a > DBL_EPSILON)
00140     {
00141       return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a));
00142     }
00143   else
00144     {
00145       return x;
00146     }
00147 }
00148 #endif  // windows substitutes for stupid microsoft compiler.
00149 
00150 
00151 inline
00152 double
00153 fsign(double x)
00154 {
00155   if (x == 0) {
00156     return 0;
00157   } else {
00158     return x/fabs(x);
00159   }
00160 }
00161 
00162 }
00163 
00164 #endif