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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (show annotations)
Mon Nov 10 01:21:39 2008 UTC (10 years, 10 months ago) by jfenwick
File MIME type: text/plain
File size: 2820 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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 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 inline
61 double
62 acosh_substitute (const double x)
63 {
64 if (x > 1.0 / SQRT_DBL_EPSILON)
65 {
66 return log (x) + M_LN2;
67 }
68 else if (x > 2)
69 {
70 return log (2 * x - 1 / (sqrt (x * x - 1) + x));
71 }
72 else if (x > 1)
73 {
74 double t = x - 1;
75 return log1p (t + sqrt (2 * t + t * t));
76 }
77 else if (x == 1)
78 {
79 return 0;
80 }
81 else
82 {
83 return FP_NAN;
84 }
85 }
86
87
88 //======================================================================
89
90 inline
91 double
92 asinh_substitute (const double x)
93 {
94 double a = fabs (x);
95 double s = (x < 0) ? -1 : 1;
96
97 if (a > 1 / SQRT_DBL_EPSILON)
98 {
99 return s * (log (a) + M_LN2);
100 }
101 else if (a > 2)
102 {
103 return s * log (2 * a + 1 / (a + sqrt (a * a + 1)));
104 }
105 else if (a > SQRT_DBL_EPSILON)
106 {
107 double a2 = a * a;
108 return s * log1p (a + a2 / (1 + sqrt (1 + a2)));
109 }
110 else
111 {
112 return x;
113 }
114 }
115
116
117 //======================================================================
118
119 inline
120 double
121 atanh_substitute (const double x)
122 {
123 double a = fabs (x);
124 double s = (x < 0) ? -1 : 1;
125
126 if (a > 1)
127 {
128 return FP_NAN;
129 }
130 else if (a == 1)
131 {
132 return (x < 0) ? -INFINITY : INFINITY;
133 }
134 else if (a >= 0.5)
135 {
136 return s * 0.5 * log1p (2 * a / (1 - a));
137 }
138 else if (a > DBL_EPSILON)
139 {
140 return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a));
141 }
142 else
143 {
144 return x;
145 }
146 }
147
148 inline
149 double
150 fsign(double x)
151 {
152 if (x == 0) {
153 return 0;
154 } else {
155 return x/fabs(x);
156 }
157 }
158
159 }
160
161 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26