/[escript]/trunk/escript/src/UnaryFuncs.h
ViewVC logotype

Contents of /trunk/escript/src/UnaryFuncs.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 3053 byte(s)
Round 1 of copyright fixes
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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26