Contents of /branches/subworld2/escriptcore/src/UnaryFuncs.h

Revision 5504 - (show annotations)
Wed Mar 4 22:58:13 2015 UTC (4 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 3121 byte(s)
```Again with a more up to date copy

```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2015 by University of Queensland 5 * http://www.uq.edu.au 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 * http://www.opensource.org/licenses/osl-3.0.php 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development 2012-2013 by School of Earth Sciences 13 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 * 15 *****************************************************************************/ 16 17 18 #if !defined escript_UnaryFuncs_20041124_H 19 #define escript_UnaryFuncs_20041124_H 20 #include "system_dep.h" 21 22 namespace escript { 23 24 #ifndef FP_NAN 25 #define FP_NAN IEEE_NaN() 26 #endif 27 28 #ifndef INFINITY 29 #define INFINITY IEEE_Infy() 30 #endif 31 32 //====================================================================== 33 34 inline 35 double log1p (const double x) 36 { 37 volatile double y; 38 y = 1 + x; 39 return log(y) - ((y-1)-x)/y ; 40 } 41 42 //====================================================================== 43 44 inline 45 float IEEE_NaN() 46 { 47 static unsigned char nan[4]={ 0, 0, 0xc0, 0x7f }; 48 return *( float *)nan; 49 } 50 51 //====================================================================== 52 53 inline 54 double IEEE_Infy() 55 { 56 static unsigned char infy[8]={ 0, 0, 0, 0, 0, 0, 0xf0, 0x7f }; 57 return *( double *)infy; 58 } 59 60 61 //====================================================================== 62 63 #if defined (_WIN32) && !defined(__INTEL_COMPILER) 64 inline 65 double 66 acosh_substitute (const double x) 67 { 68 if (x > 1.0 / SQRT_DBL_EPSILON) 69 { 70 return log (x) + M_LN2; 71 } 72 else if (x > 2) 73 { 74 return log (2 * x - 1 / (sqrt (x * x - 1) + x)); 75 } 76 else if (x > 1) 77 { 78 double t = x - 1; 79 return log1p (t + sqrt (2 * t + t * t)); 80 } 81 else if (x == 1) 82 { 83 return 0; 84 } 85 else 86 { 87 return FP_NAN; 88 } 89 } 90 91 92 //====================================================================== 93 94 inline 95 double 96 asinh_substitute (const double x) 97 { 98 double a = fabs (x); 99 double s = (x < 0) ? -1 : 1; 100 101 if (a > 1 / SQRT_DBL_EPSILON) 102 { 103 return s * (log (a) + M_LN2); 104 } 105 else if (a > 2) 106 { 107 return s * log (2 * a + 1 / (a + sqrt (a * a + 1))); 108 } 109 else if (a > SQRT_DBL_EPSILON) 110 { 111 double a2 = a * a; 112 return s * log1p (a + a2 / (1 + sqrt (1 + a2))); 113 } 114 else 115 { 116 return x; 117 } 118 } 119 120 121 //====================================================================== 122 123 inline 124 double 125 atanh_substitute (const double x) 126 { 127 double a = fabs (x); 128 double s = (x < 0) ? -1 : 1; 129 130 if (a > 1) 131 { 132 return FP_NAN; 133 } 134 else if (a == 1) 135 { 136 return (x < 0) ? -INFINITY : INFINITY; 137 } 138 else if (a >= 0.5) 139 { 140 return s * 0.5 * log1p (2 * a / (1 - a)); 141 } 142 else if (a > DBL_EPSILON) 143 { 144 return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a)); 145 } 146 else 147 { 148 return x; 149 } 150 } 151 #endif // windows substitutes for stupid microsoft compiler. 152 153 154 inline 155 double 156 fsign(double x) 157 { 158 if (x == 0) { 159 return 0; 160 } else { 161 return x/fabs(x); 162 } 163 } 164 165 } 166 167 #endif

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision