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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26