/[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

trunk/escript/src/Data/UnaryFuncs.h revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC trunk/escript/src/UnaryFuncs.h revision 2005 by jfenwick, Mon Nov 10 01:21:39 2008 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ******************************************************************************  *
4   *                                                                            *  * Copyright (c) 2003-2008 by University of Queensland
5   *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  * Earth Systems Science Computational Center (ESSCC)
6   *                                                                            *  * http://www.uq.edu.au/esscc
7   * This software is the property of ACcESS. No part of this code              *  *
8   * may be copied in any form or by any means without the expressed written    *  * Primary Business: Queensland, Australia
9   * consent of ACcESS.  Copying, use or modification of this software          *  * Licensed under the Open Software License version 3.0
10   * by any unauthorised person is illegal unless that person has a software    *  * http://www.opensource.org/licenses/osl-3.0.php
11   * license agreement with ACcESS.                                             *  *
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"
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
61    double
62    acosh_substitute (const double x)
63    {
64      if (x > 1.0 / SQRT_DBL_EPSILON)
65        {
66          return log (x) + M_LN2;
67        }
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    
88    //======================================================================
89    
90    inline
91    double
92    asinh_substitute (const double x)
93    {
94      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
120    double
121    atanh_substitute (const double x)
122    {
123      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  fsign(double x)  fsign(double x)
# Line 29  fsign(double x) Line 156  fsign(double x)
156    }    }
157  }  }
158    
159  } // end of namespace  }
160    
161  #endif  #endif

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

  ViewVC Help
Powered by ViewVC 1.1.26