escript  Revision_
UnaryFuncs.h
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
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
float IEEE_NaN()
Definition: UnaryFuncs.h:45
double log1p(const double x)
Definition: UnaryFuncs.h:35
#define M_LN2
Definition: escriptcore/src/system_dep.h:51
Definition: AbstractContinuousDomain.cpp:24
#define INFINITY
Definition: UnaryFuncs.h:29
double IEEE_Infy()
Definition: UnaryFuncs.h:54
#define FP_NAN
Definition: UnaryFuncs.h:25
#define SQRT_DBL_EPSILON
Definition: escriptcore/src/system_dep.h:47
double fsign(double x)
Definition: UnaryFuncs.h:156