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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 2934 byte(s)
Don't panic.
Updating copyright stamps

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26