escript
Revision_
|
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