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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years ago) by ksteube
File MIME type: text/plain
File size: 2779 byte(s)
Copyright updated in all files

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     double log1p (const double x)
32     {
33     volatile double y;
34     y = 1 + x;
35     return log(y) - ((y-1)-x)/y ;
36     }
37    
38     //======================================================================
39    
40     float IEEE_NaN()
41     {
42     static unsigned char nan[4]={ 0, 0, 0xc0, 0x7f };
43     return *( float *)nan;
44     }
45    
46     //======================================================================
47    
48     double IEEE_Infy()
49     {
50     static unsigned char infy[8]={ 0, 0, 0, 0, 0, 0, 0xf0, 0x7f };
51     return *( double *)infy;
52     }
53    
54    
55     //======================================================================
56    
57    
58     double
59 phornby 1032 acosh_substitute (const double x)
60 phornby 1031 {
61     if (x > 1.0 / SQRT_DBL_EPSILON)
62     {
63     return log (x) + M_LN2;
64     }
65     else if (x > 2)
66     {
67     return log (2 * x - 1 / (sqrt (x * x - 1) + x));
68     }
69     else if (x > 1)
70     {
71     double t = x - 1;
72     return log1p (t + sqrt (2 * t + t * t));
73     }
74     else if (x == 1)
75     {
76     return 0;
77     }
78     else
79     {
80     return FP_NAN;
81     }
82     }
83    
84    
85     //======================================================================
86    
87     double
88 phornby 1032 asinh_substitute (const double x)
89 phornby 1031 {
90     double a = fabs (x);
91     double s = (x < 0) ? -1 : 1;
92    
93     if (a > 1 / SQRT_DBL_EPSILON)
94     {
95     return s * (log (a) + M_LN2);
96     }
97     else if (a > 2)
98     {
99     return s * log (2 * a + 1 / (a + sqrt (a * a + 1)));
100     }
101     else if (a > SQRT_DBL_EPSILON)
102     {
103     double a2 = a * a;
104     return s * log1p (a + a2 / (1 + sqrt (1 + a2)));
105     }
106     else
107     {
108     return x;
109     }
110     }
111    
112    
113     //======================================================================
114    
115     double
116 phornby 1032 atanh_substitute (const double x)
117 phornby 1031 {
118     double a = fabs (x);
119     double s = (x < 0) ? -1 : 1;
120    
121     if (a > 1)
122     {
123     return FP_NAN;
124     }
125     else if (a == 1)
126     {
127     return (x < 0) ? -INFINITY : INFINITY;
128     }
129     else if (a >= 0.5)
130     {
131     return s * 0.5 * log1p (2 * a / (1 - a));
132     }
133     else if (a > DBL_EPSILON)
134     {
135     return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a));
136     }
137     else
138     {
139     return x;
140     }
141     }
142    
143 jgs 102 inline
144     double
145     fsign(double x)
146     {
147 jgs 123 if (x == 0) {
148     return 0;
149     } else {
150     return x/fabs(x);
151     }
152 jgs 102 }
153    
154 gross 1028 }
155    
156 jgs 102 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26