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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC revision 1031 by phornby, Wed Mar 14 06:03:21 2007 UTC
# Line 17  Line 17 
17    
18  namespace escript {  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 (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 (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 (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  inline
143  double  double
144  fsign(double x)  fsign(double x)

Legend:
Removed from v.1028  
changed lines
  Added in v.1031

  ViewVC Help
Powered by ViewVC 1.1.26