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