IT++ Logo
itcompat.cpp
Go to the documentation of this file.
1
29#include <itpp/base/itcompat.h>
30#include <itpp/base/itassert.h>
31#include <limits>
32
34
35#ifndef HAVE_TGAMMA
36// "True" gamma function
37double tgamma(double x)
38{
39 double s = (2.50662827510730 + 190.9551718930764 / (x + 1)
40 - 216.8366818437280 / (x + 2) + 60.19441764023333
41 / (x + 3) - 3.08751323928546 / (x + 4) + 0.00302963870525
42 / (x + 5) - 0.00001352385959072596 / (x + 6)) / x;
43 if (s < 0)
44 return -exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(-s));
45 else
46 return exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(s));
47}
48#endif
49
50#if !defined(HAVE_LGAMMA) || (HAVE_DECL_SIGNGAM != 1)
51// The sign of the Gamma function is returned in the external integer
52// signgam declared in <math.h>. It is 1 when the Gamma function is positive
53// or zero, -1 when it is negative. However, MinGW definition of lgamma()
54// function does not use the global signgam variable.
55int signgam;
56// Logarithm of an absolute value of gamma function
57double lgamma(double x)
58{
59 double gam = tgamma(x);
60 signgam = (gam < 0) ? -1 : 1;
61 return log(fabs(gam));
62}
63#endif
64
65#ifndef HAVE_CBRT
66// Cubic root
67double cbrt(double x) { return std::pow(x, 1.0 / 3.0); }
68#endif
69
70
71#ifndef HAVE_EXPM1
72#ifndef M_LN2
73#define M_LN2 0.69314718055994530941723212146 // ln(2)
74#endif
75// Implementation taken from GSL (http://www.gnu.org/software/gsl/)
76double expm1(double x)
77{
78 /* FIXME: this should be improved */
79 if (std::fabs(x) < M_LN2) {
80 /* Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ... */
81 double i = 1.0;
82 double sum = x;
83 double term = x / 1.0;
84 do {
85 i++;
86 term *= x / i;
87 sum += term;
88 }
89 while (std::fabs(term) > (std::fabs(sum)
90 * std::numeric_limits<double>::epsilon()));
91 return sum;
92 }
93 else {
94 return std::exp(x) - 1.0;
95 }
96}
97#endif // HAVE_EXPM1
98
99
100#ifndef HAVE_ERFC
101// Complementary error function
102double erfc(double Y)
103{
104 int ISW, I;
105 double P[4], Q[3], P1[6], Q1[5], P2[4], Q2[3];
106 double XMIN, XLARGE, SQRPI;
107 double X, RES, XSQ, XNUM, XDEN, XI, XBIG, ERFCret;
108 P[1] = 0.3166529;
109 P[2] = 1.722276;
110 P[3] = 21.38533;
111 Q[1] = 7.843746;
112 Q[2] = 18.95226;
113 P1[1] = 0.5631696;
114 P1[2] = 3.031799;
115 P1[3] = 6.865018;
116 P1[4] = 7.373888;
117 P1[5] = 4.318779e-5;
118 Q1[1] = 5.354217;
119 Q1[2] = 12.79553;
120 Q1[3] = 15.18491;
121 Q1[4] = 7.373961;
122 P2[1] = 5.168823e-2;
123 P2[2] = 0.1960690;
124 P2[3] = 4.257996e-2;
125 Q2[1] = 0.9214524;
126 Q2[2] = 0.1509421;
127 XMIN = 1.0E-5;
128 XLARGE = 4.1875E0;
129 XBIG = 9.0;
130 SQRPI = 0.5641896;
131 X = Y;
132 ISW = 1;
133 if (X < 0) {
134 ISW = -1;
135 X = -X;
136 }
137 if (X < 0.477) {
138 if (X >= XMIN) {
139 XSQ = X * X;
140 XNUM = (P[1] * XSQ + P[2]) * XSQ + P[3];
141 XDEN = (XSQ + Q[1]) * XSQ + Q[2];
142 RES = X * XNUM / XDEN;
143 }
144 else RES = X * P[3] / Q[2];
145 if (ISW == -1) RES = -RES;
146 RES = 1.0 - RES;
147 goto slut;
148 }
149 if (X > 4.0) {
150 if (ISW > 0) goto ulf;
151 if (X < XLARGE) goto eva;
152 RES = 2.0;
153 goto slut;
154 }
155 XSQ = X * X;
156 XNUM = P1[5] * X + P1[1];
157 XDEN = X + Q1[1];
158 for (I = 2;I <= 4;I++) {
159 XNUM = XNUM * X + P1[I];
160 XDEN = XDEN * X + Q1[I];
161 }
162 RES = XNUM / XDEN;
163 goto elin;
164ulf:
165 if (X > XBIG) goto fred;
166eva:
167 XSQ = X * X;
168 XI = 1.0 / XSQ;
169 XNUM = (P2[1] * XI + P2[2]) * XI + P2[3];
170 XDEN = XI + Q2[1] * XI + Q2[2];
171 RES = (SQRPI + XI * XNUM / XDEN) / X;
172elin:
173 RES = RES * exp(-XSQ);
174 if (ISW == -1) RES = 2.0 - RES;
175 goto slut;
176fred:
177 RES = 0.0;
178slut:
179 ERFCret = RES;
180 return ERFCret;
181}
182#endif
183
184
185#ifndef HAVE_ASINH
186// Arcus sinhyp
187double asinh(double x)
188{
189 return ((x >= 0) ? std::log(x + std::sqrt(x * x + 1))
190 : -std::log(-x + std::sqrt(x * x + 1)));
191}
192#endif
193
194#ifndef HAVE_ACOSH
195// Arcus coshyp
196double acosh(double x)
197{
198 it_error_if(x < 1, "acosh(): Argument out of range");
199 return std::log(x + std::sqrt(x * x - 1.0));
200}
201#endif
202
203#ifndef HAVE_ATANH
204// Arcus tanhyp
205double atanh(double x)
206{
207 it_error_if(std::fabs(x) >= 1, "atanh(): Argument out of range");
208 return 0.5 * std::log((x + 1) / (x - 1));
209}
210#endif
211
212
213#ifndef HAVE_RINT
214double rint(double x)
215{
216 // zero or NaN case
217 if ((x == 0.0) || (x != x))
218 return x;
219
220 // negative case
221 bool neg = false;
222 if (x < 0.0) {
223 x = -x;
224 neg = true;
225 }
226
227 double y = std::floor(x + 0.5);
228 int i = static_cast<int>(y);
229 if ((y - x >= 0.5) && (i & 1))
230 --y;
231
232 return neg ? -y : y;
233}
234#endif // HAVE_RINT
235
vec erfc(const vec &x)
Complementary error function.
Definition: error.cpp:240
#define it_error_if(t, s)
Abort if t is true.
Definition: itassert.h:117
vec asinh(const vec &x)
Inverse sine hyperbolic function.
Definition: trig_hyp.cpp:36
vec acosh(const vec &x)
Inverse cosine hyperbolic function.
Definition: trig_hyp.cpp:40
vec atanh(const vec &x)
Inverse tan hyperbolic function.
Definition: trig_hyp.cpp:44
vec log(const vec &x)
The natural logarithm of the elements.
Definition: log_exp.h:241
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition: log_exp.h:155
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
Error handling functions - header file.
IT++ compatibility types and functions.
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.5