/[escript]/branches/arrexp_2137_win_merge/escript/src/UnaryFuncs.h
ViewVC logotype

Annotation of /branches/arrexp_2137_win_merge/escript/src/UnaryFuncs.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2213 - (hide annotations)
Wed Jan 14 00:23:39 2009 UTC (10 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 2934 byte(s)
In preparation for merging to trunk

1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 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 ksteube 1312
14 ksteube 1811
15 jgs 102 #if !defined escript_UnaryFuncs_20041124_H
16     #define escript_UnaryFuncs_20041124_H
17 woo409 757 #include "system_dep.h"
18 jgs 102
19     namespace escript {
20    
21 phornby 1031 #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 jfenwick 2005 inline
32 phornby 1031 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 jfenwick 2005 inline
42 phornby 1031 float IEEE_NaN()
43     {
44     static unsigned char nan[4]={ 0, 0, 0xc0, 0x7f };
45     return *( float *)nan;
46     }
47    
48     //======================================================================
49    
50 jfenwick 2005 inline
51 phornby 1031 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 phornby 2049 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
61 jfenwick 2005 inline
62 phornby 1031 double
63 phornby 1032 acosh_substitute (const double x)
64 phornby 1031 {
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 jfenwick 2005 inline
92 phornby 1031 double
93 phornby 1032 asinh_substitute (const double x)
94 phornby 1031 {
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 jfenwick 2005 inline
121 phornby 1031 double
122 phornby 1032 atanh_substitute (const double x)
123 phornby 1031 {
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 phornby 2049 #endif // windows substitutes for stupid microsoft compiler.
149 phornby 1031
150 phornby 2049
151 jgs 102 inline
152     double
153     fsign(double x)
154     {
155 jgs 123 if (x == 0) {
156     return 0;
157     } else {
158     return x/fabs(x);
159     }
160 jgs 102 }
161    
162 gross 1028 }
163    
164 jgs 102 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26