00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
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