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