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

Contents of /branches/schroedinger/escript/src/UnaryFuncs.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1886 - (show annotations)
Wed Oct 15 01:34:18 2008 UTC (13 years ago) by jfenwick
File MIME type: text/plain
File size: 2820 byte(s)
Branch commit.
Added unary ops up to pos.
toString now prints expression.
Added inlines to UnaryFuncs.h.

Still only supporting DataExpanded.

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