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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1032 - (show annotations)
Wed Mar 14 06:32:09 2007 UTC (12 years, 7 months ago) by phornby
File MIME type: text/plain
File size: 3023 byte(s)
Implement inverse hyp. functions.


1 // $Id$
2 /*
3 ************************************************************
4 * Copyright 2006 by ACcESS MNRF *
5 * *
6 * http://www.access.edu.au *
7 * Primary Business: Queensland, Australia *
8 * Licensed under the Open Software License version 3.0 *
9 * http://www.opensource.org/licenses/osl-3.0.php *
10 * *
11 ************************************************************
12 */
13
14 #if !defined escript_UnaryFuncs_20041124_H
15 #define escript_UnaryFuncs_20041124_H
16 #include "system_dep.h"
17
18 namespace escript {
19
20 #ifndef FP_NAN
21 #define FP_NAN IEEE_NaN()
22 #endif
23
24 #ifndef INFINITY
25 #define INFINITY IEEE_Infy()
26 #endif
27
28 //======================================================================
29
30 double log1p (const double x)
31 {
32 volatile double y;
33 y = 1 + x;
34 return log(y) - ((y-1)-x)/y ;
35 }
36
37 //======================================================================
38
39 float IEEE_NaN()
40 {
41 static unsigned char nan[4]={ 0, 0, 0xc0, 0x7f };
42 return *( float *)nan;
43 }
44
45 //======================================================================
46
47 double IEEE_Infy()
48 {
49 static unsigned char infy[8]={ 0, 0, 0, 0, 0, 0, 0xf0, 0x7f };
50 return *( double *)infy;
51 }
52
53
54 //======================================================================
55
56
57 double
58 acosh_substitute (const double x)
59 {
60 if (x > 1.0 / SQRT_DBL_EPSILON)
61 {
62 return log (x) + M_LN2;
63 }
64 else if (x > 2)
65 {
66 return log (2 * x - 1 / (sqrt (x * x - 1) + x));
67 }
68 else if (x > 1)
69 {
70 double t = x - 1;
71 return log1p (t + sqrt (2 * t + t * t));
72 }
73 else if (x == 1)
74 {
75 return 0;
76 }
77 else
78 {
79 return FP_NAN;
80 }
81 }
82
83
84 //======================================================================
85
86 double
87 asinh_substitute (const double x)
88 {
89 double a = fabs (x);
90 double s = (x < 0) ? -1 : 1;
91
92 if (a > 1 / SQRT_DBL_EPSILON)
93 {
94 return s * (log (a) + M_LN2);
95 }
96 else if (a > 2)
97 {
98 return s * log (2 * a + 1 / (a + sqrt (a * a + 1)));
99 }
100 else if (a > SQRT_DBL_EPSILON)
101 {
102 double a2 = a * a;
103 return s * log1p (a + a2 / (1 + sqrt (1 + a2)));
104 }
105 else
106 {
107 return x;
108 }
109 }
110
111
112 //======================================================================
113
114 double
115 atanh_substitute (const double x)
116 {
117 double a = fabs (x);
118 double s = (x < 0) ? -1 : 1;
119
120 if (a > 1)
121 {
122 return FP_NAN;
123 }
124 else if (a == 1)
125 {
126 return (x < 0) ? -INFINITY : INFINITY;
127 }
128 else if (a >= 0.5)
129 {
130 return s * 0.5 * log1p (2 * a / (1 - a));
131 }
132 else if (a > DBL_EPSILON)
133 {
134 return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a));
135 }
136 else
137 {
138 return x;
139 }
140 }
141
142 inline
143 double
144 fsign(double x)
145 {
146 if (x == 0) {
147 return 0;
148 } else {
149 return x/fabs(x);
150 }
151 }
152
153 }
154
155 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26