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

Revision 1032 - (show annotations)
Wed Mar 14 06:32:09 2007 UTC (12 years, 6 months ago) by phornby
File MIME type: text/plain
File size: 3023 byte(s)
Implement inverse hyp. functions.

 1 // \$Id\$ 2 /* 3 ************************************************************ 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 */ 13 14 #if !defined escript_UnaryFuncs_20041124_H 15 #define escript_UnaryFuncs_20041124_H 16 #include "system_dep.h" 17 18 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_substitute (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_substitute (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_substitute (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 143 double 144 fsign(double x) 145 { 146 if (x == 0) { 147 return 0; 148 } else { 149 return x/fabs(x); 150 } 151 } 152 153 } 154 155 #endif

## Properties

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