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