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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1
2 /*******************************************************
3 *
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
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 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 acosh_substitute (const double x)
60 {
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 asinh_substitute (const double x)
89 {
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 atanh_substitute (const double x)
117 {
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 inline
144 double
145 fsign(double x)
146 {
147 if (x == 0) {
148 return 0;
149 } else {
150 return x/fabs(x);
151 }
152 }
153
154 }
155
156 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26