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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1032 - (hide 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 jgs 102 // $Id$
2     /*
3 elspeth 615 ************************************************************
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 jgs 102 */
13    
14     #if !defined escript_UnaryFuncs_20041124_H
15     #define escript_UnaryFuncs_20041124_H
16 woo409 757 #include "system_dep.h"
17 jgs 102
18     namespace escript {
19    
20 phornby 1031 #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 phornby 1032 acosh_substitute (const double x)
59 phornby 1031 {
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 phornby 1032 asinh_substitute (const double x)
88 phornby 1031 {
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 phornby 1032 atanh_substitute (const double x)
116 phornby 1031 {
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 jgs 102 inline
143     double
144     fsign(double x)
145     {
146 jgs 123 if (x == 0) {
147     return 0;
148     } else {
149     return x/fabs(x);
150     }
151 jgs 102 }
152    
153 gross 1028 }
154    
155 jgs 102 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26