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

Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 6 months ago) by ksteube
File MIME type: text/plain
File size: 2779 byte(s)
```Copyright updated in all files

```
 1 ksteube 1312 2 /******************************************************* 3 ksteube 1811 * 4 * Copyright (c) 2003-2008 by University of Queensland 5 * Earth Systems Science Computational Center (ESSCC) 6 * http://www.uq.edu.au/esscc 7 * 8 * Primary Business: Queensland, Australia 9 * Licensed under the Open Software License version 3.0 10 * http://www.opensource.org/licenses/osl-3.0.php 11 * 12 *******************************************************/ 13 ksteube 1312 14 ksteube 1811 15 jgs 102 #if !defined escript_UnaryFuncs_20041124_H 16 #define escript_UnaryFuncs_20041124_H 17 woo409 757 #include "system_dep.h" 18 jgs 102 19 namespace escript { 20 21 phornby 1031 #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 double log1p (const double x) 32 { 33 volatile double y; 34 y = 1 + x; 35 return log(y) - ((y-1)-x)/y ; 36 } 37 38 //====================================================================== 39 40 float IEEE_NaN() 41 { 42 static unsigned char nan[4]={ 0, 0, 0xc0, 0x7f }; 43 return *( float *)nan; 44 } 45 46 //====================================================================== 47 48 double IEEE_Infy() 49 { 50 static unsigned char infy[8]={ 0, 0, 0, 0, 0, 0, 0xf0, 0x7f }; 51 return *( double *)infy; 52 } 53 54 55 //====================================================================== 56 57 58 double 59 phornby 1032 acosh_substitute (const double x) 60 phornby 1031 { 61 if (x > 1.0 / SQRT_DBL_EPSILON) 62 { 63 return log (x) + M_LN2; 64 } 65 else if (x > 2) 66 { 67 return log (2 * x - 1 / (sqrt (x * x - 1) + x)); 68 } 69 else if (x > 1) 70 { 71 double t = x - 1; 72 return log1p (t + sqrt (2 * t + t * t)); 73 } 74 else if (x == 1) 75 { 76 return 0; 77 } 78 else 79 { 80 return FP_NAN; 81 } 82 } 83 84 85 //====================================================================== 86 87 double 88 phornby 1032 asinh_substitute (const double x) 89 phornby 1031 { 90 double a = fabs (x); 91 double s = (x < 0) ? -1 : 1; 92 93 if (a > 1 / SQRT_DBL_EPSILON) 94 { 95 return s * (log (a) + M_LN2); 96 } 97 else if (a > 2) 98 { 99 return s * log (2 * a + 1 / (a + sqrt (a * a + 1))); 100 } 101 else if (a > SQRT_DBL_EPSILON) 102 { 103 double a2 = a * a; 104 return s * log1p (a + a2 / (1 + sqrt (1 + a2))); 105 } 106 else 107 { 108 return x; 109 } 110 } 111 112 113 //====================================================================== 114 115 double 116 phornby 1032 atanh_substitute (const double x) 117 phornby 1031 { 118 double a = fabs (x); 119 double s = (x < 0) ? -1 : 1; 120 121 if (a > 1) 122 { 123 return FP_NAN; 124 } 125 else if (a == 1) 126 { 127 return (x < 0) ? -INFINITY : INFINITY; 128 } 129 else if (a >= 0.5) 130 { 131 return s * 0.5 * log1p (2 * a / (1 - a)); 132 } 133 else if (a > DBL_EPSILON) 134 { 135 return s * 0.5 * log1p (2 * a + 2 * a * a / (1 - a)); 136 } 137 else 138 { 139 return x; 140 } 141 } 142 143 jgs 102 inline 144 double 145 fsign(double x) 146 { 147 jgs 123 if (x == 0) { 148 return 0; 149 } else { 150 return x/fabs(x); 151 } 152 jgs 102 } 153 154 gross 1028 } 155 156 jgs 102 #endif

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26