/[escript]/branches/subworld2/escriptcore/src/UnaryFuncs.h
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26