1 // -*- C++ -*- 2 /** 3 * @brief OpenTURNS wrapper to a library of special functions 4 * 5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca 6 * 7 * This library is free software: you can redistribute it and/or modify 8 * it under the terms of the GNU Lesser General Public License as published by 9 * the Free Software Foundation, either version 3 of the License, or 10 * (at your option) any later version. 11 * 12 * This library is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 * GNU Lesser General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General Public License 18 * along with this library. If not, see <http://www.gnu.org/licenses/>. 19 * 20 */ 21 #ifndef OPENTURNS_SPECFUNC_HXX 22 #define OPENTURNS_SPECFUNC_HXX 23 24 #include "openturns/OTprivate.hxx" 25 26 BEGIN_NAMESPACE_OPENTURNS 27 28 class OT_API SpecFunc 29 { 30 31 public: 32 33 // 0.39894228040143267 = 1 / sqrt(2.pi) 34 static const Scalar ISQRT2PI; 35 // 2.5066282746310005024 = sqrt(2.pi) 36 static const Scalar SQRT2PI; 37 // 0.91893853320467274177 = log(sqrt(2.pi)) 38 static const Scalar LOGSQRT2PI; 39 // 0.57721566490153286 = Euler constant gamma 40 static const Scalar EulerConstant; 41 // 1.64493406684822643 = pi^2 / 6 42 static const Scalar PI2_6; 43 // 1.28254983016118640 = pi / sqrt(6) 44 static const Scalar PI_SQRT6; 45 // 0.45005320754569466 = gamma * sqrt(6) / pi 46 static const Scalar EULERSQRT6_PI; 47 // 3.28986813369645287 = pi^2 / 3 48 static const Scalar PI2_3; 49 // 0.55132889542179204 = sqrt(3) / pi 50 static const Scalar SQRT3_PI; 51 // 1.81379936423421785 = pi / sqrt(3) 52 static const Scalar PI_SQRT3; 53 // 6.283185307179586476925286 = 2*pi 54 static const Scalar TWOPI; 55 // 1.20205690315959429 = Zeta(3) 56 static const Scalar ZETA3; 57 // Maximum number of iterations for algorithms 58 static const UnsignedInteger MaximumIteration; 59 // Maximum precision for algorithms 60 static const Scalar Precision; 61 // Minimum positive real number 62 static const Scalar MinScalar; 63 static const Scalar LogMinScalar; 64 // Maximum positive real number 65 static const Scalar MaxScalar; 66 static const Scalar LogMaxScalar; 67 // Minimum negative real number; 68 static const Scalar LowestScalar; 69 // Real number accuracy 70 static const Scalar ScalarEpsilon; 71 72 // Information about capabilities 73 static Bool IsBoostAvailable(); 74 static Bool IsMPFRAvailable(); 75 static Bool IsMPCAvailable(); 76 77 // Some facilities for NaN and inf 78 static Bool IsNaN(const Scalar value); 79 static Bool IsInf(const Scalar value); 80 static Bool IsNormal(const Scalar value); 81 82 // Modified first kind Bessel function of order 0: BesselI0(x) = \sum_{m=0}\infty\frac{1}{m!^2}\left(\frac{x}{2}\right)^{2m} 83 private: 84 static Scalar SmallCaseBesselI0(const Scalar x); 85 static Scalar LargeCaseLogBesselI0(const Scalar x); 86 public: 87 static Scalar BesselI0(const Scalar x); 88 static Scalar LogBesselI0(const Scalar x); 89 // Modified first kind Bessel function of order 1: BesselI1(x) = \sum_{m=0}\infty\frac{1}{m!(m+1)!}\left(\frac{x}{2}\right)^{2m+1} 90 private: 91 static Scalar SmallCaseBesselI1(const Scalar x); 92 static Scalar LargeCaseLogBesselI1(const Scalar x); 93 public: 94 static Scalar BesselI1(const Scalar x); 95 static Scalar LogBesselI1(const Scalar x); 96 // Difference between the logarithms of BesselI1 and BesselI0: 97 // DeltaLogBesselI10(x) = log(BesselI1(x)) - log(BesselI0(x)) 98 private: 99 static Scalar LargeCaseDeltaLogBesselI10(const Scalar x); 100 public: 101 static Scalar DeltaLogBesselI10(const Scalar x); 102 // Modified second kind Bessel function of order nu: BesselK(nu, x)=\frac{\pi}{2}\frac{I_{-\nu}(x)-I_[\nu}(x)}{\sin{\nu\pi}} 103 static Scalar LogBesselK(const Scalar nu, 104 const Scalar x); 105 static Scalar BesselK(const Scalar nu, 106 const Scalar x); 107 static Scalar BesselKDerivative(const Scalar nu, 108 const Scalar x); 109 // Beta function: beta(a, b) = \int_0^1 t^{a-1}.(1-t)^{b-1} dt 110 static Scalar Beta(const Scalar a, 111 const Scalar b); 112 // Incomplete beta function: betaInc(a, b, x) = \int_0^x t^{a-1}.(1-t)^{b-1} dt 113 static Scalar IncompleteBeta(const Scalar a, 114 const Scalar b, 115 const Scalar x, 116 const Bool tail = false); 117 // Incomplete beta function inverse with respect to x 118 static Scalar IncompleteBetaInverse(const Scalar a, 119 const Scalar b, 120 const Scalar x, 121 const Bool tail = false); 122 // Incomplete beta ratio function: betaRatioInc(a, b, x) = \int_0^x t^{a-1}.(1-t)^{b-1} dt / beta(a, b) 123 static Scalar RegularizedIncompleteBeta(const Scalar a, 124 const Scalar b, 125 const Scalar x, 126 const Bool tail = false); 127 // Incomplete beta ratio function inverse with respect to x 128 static Scalar RegularizedIncompleteBetaInverse(const Scalar a, 129 const Scalar b, 130 const Scalar x, 131 const Bool tail = false); 132 // Natural logarithm of the beta function 133 static Scalar LnBeta(const Scalar a, 134 const Scalar b); 135 static Scalar LogBeta(const Scalar a, 136 const Scalar b); 137 // Dawson function: Dawson(x) = \exp(-x^2) * \int_0^x \exp(t^2) dt 138 static Scalar Dawson(const Scalar x); 139 static Complex Dawson(const Complex & z); 140 // Debye function of order n: DebyeN(x, n) = n / x^n \int_0^x t^n/(\exp(t)-1) dt 141 static Scalar Debye(const Scalar x, 142 const UnsignedInteger n); 143 // DiLog function: Dilog(x) = -\int_0^x \log(1-t)/t dt 144 static Scalar DiLog(const Scalar x); 145 // Exponential integral function: Ei(x) = -\int_{-x}^{\infty}exp(-t)/t dt 146 static Scalar Ei(const Scalar x); 147 // Complex exponential integral function: Ei(z) = -\int_{-z}^{\infty}exp(-t)/t dt 148 static Complex Ei(const Complex & z); 149 // Complex Faddeeva function: Faddeeva(z) = exp(-z^2)\erfc(-I*z) 150 static Complex Faddeeva(const Complex & z); 151 // Imaginary part of the Faddeeva function: FaddeevaIm(z) = Im(Faddeeva(x)) 152 static Scalar FaddeevaIm(const Scalar x); 153 // Factorial and log-Factorial functions 154 static Scalar Factorial(UnsignedInteger n); 155 static Scalar LogFactorial(UnsignedInteger n); 156 // Gamma function: gamma(a) = \int_0^{\infty} t^{a-1}\exp(-t) dt 157 static Scalar Gamma(const Scalar a); 158 // igamma1pm1(a) = 1 / gamma(1 + a) - 1 159 static Scalar IGamma1pm1(const Scalar a); 160 // GammaCorrection(a) = LogGamma(a) - log(sqrt(2.Pi)) + a - (a - 1/2) log(a) 161 static Scalar GammaCorrection(const Scalar a); 162 // Complex gamma function: gamma(a) = \int_0^{\infty} t^{a-1}\exp(-t) dt 163 static Complex Gamma(const Complex & a); 164 // Natural logarithm of the gamma function 165 static Scalar LnGamma(const Scalar a); 166 static Scalar LogGamma(const Scalar a); 167 static Scalar LogGamma1p(const Scalar a); 168 static Complex LogGamma(const Complex & a); 169 // Incomplete gamma function: gamma(a, x) = \int_0^x t^{a-1}\exp(-t) dt 170 static Scalar IncompleteGamma(const Scalar a, 171 const Scalar x, 172 const Bool tail = false); 173 // Incomplete gamma function inverse with respect to x 174 static Scalar IncompleteGammaInverse(const Scalar a, 175 const Scalar x, 176 const Bool tail = false); 177 // Regularized incomplete gamma function: gamma(a, x) = \int_0^x t^{a-1}\exp(-t) dt / \Gamma(a) 178 static Scalar RegularizedIncompleteGamma(const Scalar a, 179 const Scalar x, 180 const Bool tail = false); 181 // Regularized incomplete gamma function inverse with respect to x 182 static Scalar RegularizedIncompleteGammaInverse(const Scalar a, 183 const Scalar x, 184 const Bool tail = false); 185 // Digamma function: psi(x) = ((dgamma/dx) / gamma)(x) 186 static Scalar DiGamma(const Scalar x); 187 static Scalar Psi(const Scalar x); 188 // Inverse of the DiGamma function 189 static Scalar DiGammaInv(const Scalar a); 190 // Trigamma function: TriGamma(x) = ((d^2gamma/dx^2) / gamma)(x) 191 static Scalar TriGamma(const Scalar x); 192 // Stirling error: Stirlerr(n) = log(n!) - log( sqrt(2*pi*n)*(n/e)^n ) 193 static Scalar Stirlerr(const UnsignedInteger n); 194 // Hypergeometric function of type (1,1): hyperGeom_1_1(p1, q1, x) = \sum_{n=0}^{\infty} [\prod_{k=0}^{n-1} (p1 + k) / (q1 + k)] * x^n / n! 195 static Scalar HyperGeom_1_1(const Scalar p1, 196 const Scalar q1, 197 const Scalar x); 198 // Complex hypergeometric function of type (1,1): hyperGeom_1_1(p1, q1, x) = \sum_{n=0}^{\infty} [\prod_{k=0}^{n-1} (p1 + k) / (q1 + k)] * x^n / n! 199 static Complex HyperGeom_1_1(const Scalar p1, 200 const Scalar q1, 201 const Complex & x); 202 // Hypergeometric function of type (2,1): hyperGeom_2_1(p1, p2, q1, x) = \sum_{n=0}^{\infty} [\prod_{k=0}^{n-1} (p1 + k) . (p2 + k) / (q1 + k)] * x^n / n! 203 static Scalar HyperGeom_2_1(const Scalar p1, 204 const Scalar p2, 205 const Scalar q1, 206 const Scalar x); 207 // Hypergeometric function of type (2,2): hyperGeom_2_1(p1, p2, q1, q2, x) = \sum_{n=0}^{\infty} [\prod_{k=0}^{n-1} (p1 + k) . (p2 + k) / (q1 + k) / (q2 + k)] * x^n / n! 208 static Scalar HyperGeom_2_2(const Scalar p1, 209 const Scalar p2, 210 const Scalar q1, 211 const Scalar q2, 212 const Scalar x); 213 // Erf function erf(x) = 2 / \sqrt(\pi) . \int_0^x \exp(-t^2) dt 214 static Scalar Erf(const Scalar x); 215 static Complex Erf(const Complex & z); 216 // Erf function erfi(x) = -i.erf(iz) 217 static Scalar ErfI(const Scalar x); 218 static Complex ErfI(const Complex & z); 219 // Erf function erfc(x) = 1 - erf(x) 220 static Scalar ErfC(const Scalar x); 221 static Complex ErfC(const Complex & z); 222 // Erf function erfcx(x) = exp(x^2).erfc(x) 223 static Scalar ErfCX(const Scalar x); 224 static Complex ErfCX(const Complex & z); 225 // Inverse of the erf function 226 static Scalar ErfInverse(const Scalar x); 227 // Real branch of Lambert W function (principal or secndary) 228 static Scalar LambertW(const Scalar x, 229 const Bool principal = true); 230 // Accurate value of log(1+z) for |z|<<1 231 static Complex Log1p(const Complex & z); 232 // Accurate value of exp(z)-1 for |z|<<1 233 static Complex Expm1(const Complex & z); 234 // Accurate value of log(1-exp(-x)) for all x 235 static Complex Log1MExp(const Scalar x); 236 // MarcumQ- function 237 // static Scalar MarcumQFunction(const Scalar a,const Scalar b); 238 239 // Next power of two 240 static UnsignedInteger NextPowerOfTwo(const UnsignedInteger n); 241 242 // Integer power 243 static Scalar IPow(const Scalar x, const SignedInteger n); 244 245 // Integer root 246 static Scalar IRoot(const Scalar x, const SignedInteger n); 247 248 // Integer log2 249 static UnsignedInteger Log2(const Unsigned64BitsInteger n); 250 251 // Compute the number of bits sets to 1 in n 252 // Best known algorithm for 64 bits n and fast multiply 253 static UnsignedInteger BitCount(const Unsigned64BitsInteger n); 254 255 // Missing functions in cmath wrt math.h as of C++98 256 static Scalar Acosh(const Scalar x); 257 static Scalar Asinh(const Scalar x); 258 static Scalar Atanh(const Scalar x); 259 static Scalar Cbrt(const Scalar x); 260 261 // binomial coefficient C(n, k) 262 static UnsignedInteger BinomialCoefficient(const UnsignedInteger n, 263 const UnsignedInteger k); 264 }; /* class SpecFunc */ 265 266 END_NAMESPACE_OPENTURNS 267 268 #endif /* OPENTURNS_SPECFUNC_HXX */ 269