/[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 2005 by jfenwick, Mon Nov 10 01:21:39 2008 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2008 by University of Queensland
5   *                                                          *  * Earth Systems Science Computational Center (ESSCC)
6   *              http://www.access.edu.au                    *  * http://www.uq.edu.au/esscc
7   *       Primary Business: Queensland, Australia            *  *
8   *  Licensed under the Open Software License version 3.0    *  * Primary Business: Queensland, Australia
9   *     http://www.opensource.org/licenses/osl-3.0.php       *  * 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  #if !defined escript_UnaryFuncs_20041124_H
16  #define escript_UnaryFuncs_20041124_H  #define escript_UnaryFuncs_20041124_H
17  #include "system_dep.h"  #include "system_dep.h"
18    
19  namespace escript {  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  inline
61  double  double
62  fsign(double x)  acosh_substitute (const double x)
63  {  {
64    if (x == 0) {    if (x > 1.0 / SQRT_DBL_EPSILON)
65      return 0;      {
66    } else {        return log (x) + M_LN2;
67      return x/fabs(x);      }
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  /* substitute functions for _WIN32 */  
88    //======================================================================
89    
90  inline  inline
91  double  double
92  asinh_substitute(double x)  asinh_substitute (const double x)
93  {  {
94      return 0;    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  inline
120  double  double
121  acosh_substitute(double x)  atanh_substitute (const double x)
122  {  {
123      return 0;    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  inline
149  double  double
150  atanh_substitute(double x)  fsign(double x)
151  {  {
152      if (x == 0) {
153      return 0;      return 0;
154      } else {
155        return x/fabs(x);
156      }
157  }  }
158    
159    }
160    
 } // end of namespace  
161  #endif  #endif

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

  ViewVC Help
Powered by ViewVC 1.1.26