36#ifndef VIGRA_MATHUTIL_HXX
37#define VIGRA_MATHUTIL_HXX
40# pragma warning (disable: 4996)
49#include "sized_int.hxx"
50#include "numerictraits.hxx"
51#include "algorithm.hxx"
90# define M_PI 3.14159265358979323846
94# define M_2_PI 0.63661977236758134308
98# define M_PI_2 1.57079632679489661923
102# define M_PI_4 0.78539816339744830962
106# define M_SQRT2 1.41421356237309504880
110# define M_E 2.71828182845904523536
114# define M_EULER_GAMMA 0.5772156649015329
126using VIGRA_CSTD::pow;
127using VIGRA_CSTD::floor;
128using VIGRA_CSTD::ceil;
129using VIGRA_CSTD::exp;
138#define VIGRA_DEFINE_UNSIGNED_ABS(T) \
139 inline T abs(T t) { return t; }
141VIGRA_DEFINE_UNSIGNED_ABS(
bool)
142VIGRA_DEFINE_UNSIGNED_ABS(
unsigned char)
143VIGRA_DEFINE_UNSIGNED_ABS(
unsigned short)
144VIGRA_DEFINE_UNSIGNED_ABS(
unsigned int)
145VIGRA_DEFINE_UNSIGNED_ABS(
unsigned long)
146VIGRA_DEFINE_UNSIGNED_ABS(
unsigned long long)
148#undef VIGRA_DEFINE_UNSIGNED_ABS
150#define VIGRA_DEFINE_MISSING_ABS(T) \
151 inline T abs(T t) { return t < 0 ? static_cast<T>(-t) : t; }
153VIGRA_DEFINE_MISSING_ABS(
signed char)
154VIGRA_DEFINE_MISSING_ABS(
signed short)
156#if defined(_MSC_VER) && _MSC_VER < 1600
157VIGRA_DEFINE_MISSING_ABS(
signed long long)
160#undef VIGRA_DEFINE_MISSING_ABS
171inline bool isinf(REAL v)
173 return _finite(v) == 0 && _isnan(v) == 0;
177inline bool isnan(REAL v)
179 return _isnan(v) != 0;
183inline bool isfinite(REAL v)
185 return _finite(v) != 0;
193#define VIGRA_DEFINE_SCALAR_DOT(T) \
194 inline NumericTraits<T>::Promote dot(T l, T r) { return l*r; }
196VIGRA_DEFINE_SCALAR_DOT(
unsigned char)
197VIGRA_DEFINE_SCALAR_DOT(
unsigned short)
198VIGRA_DEFINE_SCALAR_DOT(
unsigned int)
199VIGRA_DEFINE_SCALAR_DOT(
unsigned long)
200VIGRA_DEFINE_SCALAR_DOT(
unsigned long long)
201VIGRA_DEFINE_SCALAR_DOT(
signed char)
202VIGRA_DEFINE_SCALAR_DOT(
signed short)
203VIGRA_DEFINE_SCALAR_DOT(
signed int)
204VIGRA_DEFINE_SCALAR_DOT(
signed long)
205VIGRA_DEFINE_SCALAR_DOT(
signed long long)
206VIGRA_DEFINE_SCALAR_DOT(
float)
207VIGRA_DEFINE_SCALAR_DOT(
double)
208VIGRA_DEFINE_SCALAR_DOT(
long double)
210#undef VIGRA_DEFINE_SCALAR_DOT
216inline float pow(
float v,
double e)
218 return std::pow(v, (
float)e);
221inline long double pow(
long double v,
double e)
223 return std::pow(v, (
long double)e);
238inline float round(
float t)
245inline double round(
double t)
252inline long double round(
long double t)
271 ? (
long long)(t + 0.5)
272 : (
long long)(t - 0.5);
280 return x && !(x & (x - 1));
332 static Int32 table[64];
336Int32 IntLog2<T>::table[64] = {
337 -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21,
338 29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6,
339 -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9,
340 11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25,
341 -1, 7, 24, -1, 23, -1, 31, -1};
370 return detail::IntLog2<Int32>::table[x >> 26];
382typename NumericTraits<T>::Promote
sq(T t)
389template <
class V,
unsigned>
392 static V call(
const V & x,
const V & y) {
return x * y; }
395struct cond_mult<V, 0>
397 static V call(
const V &,
const V & y) {
return y; }
400template <
class V,
unsigned n>
403 static V call(
const V & x)
406 ? cond_mult<V, n & 1>::call(x, power_static<V, n / 2>::call(x * x))
411struct power_static<V, 0>
413 static V call(
const V & )
426template <
unsigned n,
class V>
429 return detail::power_static<V, n>::call(x);
438 static UInt32 sqq_table[];
443UInt32 IntSquareRoot<T>::sqq_table[] = {
444 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57,
445 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83,
446 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102,
447 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
448 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
449 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
450 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
451 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
452 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
453 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
454 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
455 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
456 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
457 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
458 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
459 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
460 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
461 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
472 if (x >= 0x40000000) {
475 xn = sqq_table[x>>24] << 8;
477 xn = sqq_table[x>>22] << 7;
480 xn = sqq_table[x>>20] << 6;
482 xn = sqq_table[x>>18] << 5;
486 xn = sqq_table[x>>16] << 4;
488 xn = sqq_table[x>>14] << 3;
491 xn = sqq_table[x>>12] << 2;
493 xn = sqq_table[x>>10] << 1;
501 xn = (sqq_table[x>>8] >> 0) + 1;
503 xn = (sqq_table[x>>6] >> 1) + 1;
506 xn = (sqq_table[x>>4] >> 2) + 1;
508 xn = (sqq_table[x>>2] >> 3) + 1;
512 return sqq_table[x] >> 4;
516 xn = (xn + 1 + x / xn) / 2;
518 xn = (xn + 1 + x / xn) / 2;
529using VIGRA_CSTD::sqrt;
541 throw std::domain_error(
"sqrti(Int32): negative argument.");
542 return (
Int32)detail::IntSquareRoot<UInt32>::exec((
UInt32)v);
554 return detail::IntSquareRoot<UInt32>::exec(v);
566inline double hypot(
double a,
double b)
568 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b);
570 return absa * VIGRA_CSTD::sqrt(1.0 +
sq(absb/absa));
574 : absb * VIGRA_CSTD::sqrt(1.0 +
sq(absa/absb));
593 return t > NumericTraits<T>::zero()
594 ? NumericTraits<T>::one()
595 : t < NumericTraits<T>::zero()
596 ? -NumericTraits<T>::one()
597 : NumericTraits<T>::zero();
610 return t > NumericTraits<T>::zero()
612 : t < NumericTraits<T>::zero()
624template <
class T1,
class T2>
627 return t2 >= NumericTraits<T2>::zero()
648#define VIGRA_DEFINE_ODD_EVEN(T) \
649 inline bool even(T t) { return (t&1) == 0; } \
650 inline bool odd(T t) { return (t&1) == 1; }
652VIGRA_DEFINE_ODD_EVEN(
char)
653VIGRA_DEFINE_ODD_EVEN(
short)
654VIGRA_DEFINE_ODD_EVEN(
int)
655VIGRA_DEFINE_ODD_EVEN(
long)
656VIGRA_DEFINE_ODD_EVEN(
long long)
657VIGRA_DEFINE_ODD_EVEN(
unsigned char)
658VIGRA_DEFINE_ODD_EVEN(
unsigned short)
659VIGRA_DEFINE_ODD_EVEN(
unsigned int)
660VIGRA_DEFINE_ODD_EVEN(
unsigned long)
661VIGRA_DEFINE_ODD_EVEN(
unsigned long long)
663#undef VIGRA_DEFINE_ODD_EVEN
665#define VIGRA_DEFINE_NORM(T) \
666 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
667 inline NormTraits<T>::NormType norm(T t) { return abs(t); }
669VIGRA_DEFINE_NORM(
bool)
670VIGRA_DEFINE_NORM(
signed char)
671VIGRA_DEFINE_NORM(
unsigned char)
672VIGRA_DEFINE_NORM(
short)
673VIGRA_DEFINE_NORM(
unsigned short)
674VIGRA_DEFINE_NORM(
int)
675VIGRA_DEFINE_NORM(
unsigned int)
676VIGRA_DEFINE_NORM(
long)
677VIGRA_DEFINE_NORM(
unsigned long)
678VIGRA_DEFINE_NORM(
long long)
679VIGRA_DEFINE_NORM(
unsigned long long)
680VIGRA_DEFINE_NORM(
float)
681VIGRA_DEFINE_NORM(
double)
682VIGRA_DEFINE_NORM(
long double)
684#undef VIGRA_DEFINE_NORM
687inline typename NormTraits<std::complex<T> >::SquaredNormType
690 return sq(t.real()) +
sq(t.imag());
716inline typename NormTraits<T>::NormType
719 typedef typename NormTraits<T>::SquaredNormType SNT;
720 return sqrt(
static_cast<typename SquareRootTraits<SNT>::SquareRootArgument
>(
squaredNorm(t)));
736 double d =
hypot(a00 - a11, 2.0*a01);
737 *r0 =
static_cast<T
>(0.5*(a00 + a11 + d));
738 *r1 =
static_cast<T
>(0.5*(a00 + a11 - d));
755 T * r0, T * r1, T * r2)
757 double inv3 = 1.0 / 3.0, root3 = std::sqrt(3.0);
759 double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
760 double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
761 double c2 = a00 + a11 + a22;
762 double c2Div3 = c2*inv3;
763 double aDiv3 = (c1 - c2*c2Div3)*inv3;
766 double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
767 double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
770 double magnitude = std::sqrt(-aDiv3);
771 double angle = std::atan2(std::sqrt(-q),mbDiv2)*inv3;
772 double cs = std::cos(angle);
773 double sn = std::sin(angle);
774 *r0 =
static_cast<T
>(c2Div3 + 2.0*magnitude*cs);
775 *r1 =
static_cast<T
>(c2Div3 - magnitude*(cs + root3*sn));
776 *r2 =
static_cast<T
>(c2Div3 - magnitude*(cs - root3*sn));
788T ellipticRD(T x, T y, T z)
790 double f = 1.0, s = 0.0, X, Y, Z, m;
793 m = (x + y + 3.0*z) / 5.0;
797 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
799 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
800 s += f / (VIGRA_CSTD::sqrt(z)*(z + l));
809 double d = a - 6.0*b;
810 double e = d + 2.0*c;
811 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0)
812 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5);
816T ellipticRF(T x, T y, T z)
821 m = (x + y + z) / 3.0;
825 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
827 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
832 double d = X*Y -
sq(Z);
834 return (1.0 - d/10.0 + p/14.0 +
sq(d)/24.0 - d*p*3.0/44.0) / VIGRA_CSTD::sqrt(m);
857 double c2 =
sq(VIGRA_CSTD::cos(x));
858 double s = VIGRA_CSTD::sin(x);
859 return s*detail::ellipticRF(c2, 1.0 -
sq(k*s), 1.0);
881 double c2 =
sq(VIGRA_CSTD::cos(x));
882 double s = VIGRA_CSTD::sin(x);
884 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0));
887#if defined(_MSC_VER) && _MSC_VER < 1800
894 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
895 double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
896 t*(0.09678418+t*(-0.18628806+t*(0.27886807+
897 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
898 t*0.17087277)))))))));
922inline double erf(
double x)
924 return detail::erfImpl(x);
938 double a = degreesOfFreedom + noncentrality;
939 double b = (a + noncentrality) /
sq(a);
940 double t = (VIGRA_CSTD::pow((
double)
arg / a, 1.0/3.0) - (1.0 - 2.0 / 9.0 * b)) / VIGRA_CSTD::sqrt(2.0 / 9.0 * b);
941 return 0.5*(1.0 + erf(t/VIGRA_CSTD::sqrt(2.0)));
945void noncentralChi2OneIteration(T
arg, T & lans, T & dans, T & pans,
unsigned int & j)
950 lans = lans + VIGRA_CSTD::log(
arg / j);
951 dans = VIGRA_CSTD::exp(lans);
955 dans = dans *
arg / j;
962std::pair<double, double>
noncentralChi2CDF(
unsigned int degreesOfFreedom, T noncentrality, T
arg, T eps)
964 vigra_precondition(noncentrality >= 0.0 &&
arg >= 0.0 && eps > 0.0,
965 "noncentralChi2P(): parameters must be positive.");
966 if (
arg == 0.0 && degreesOfFreedom > 0)
967 return std::make_pair(0.0, 0.0);
970 double b1 = 0.5 * noncentrality,
971 ao = VIGRA_CSTD::exp(-b1),
973 lnrtpi2 = 0.22579135264473,
974 probability, density, lans, dans, pans,
sum, am, hold;
975 unsigned int maxit = 500,
977 if(degreesOfFreedom % 2)
980 lans = -0.5 * (
arg + VIGRA_CSTD::log(
arg)) - lnrtpi2;
981 dans = VIGRA_CSTD::exp(lans);
982 pans = erf(VIGRA_CSTD::sqrt(
arg/2.0));
988 dans = VIGRA_CSTD::exp(lans);
993 if(degreesOfFreedom == 0)
996 degreesOfFreedom = 2;
998 sum = 1.0 / ao - 1.0 - am;
1000 probability = 1.0 + am * pans;
1005 degreesOfFreedom = degreesOfFreedom - 1;
1007 sum = 1.0 / ao - 1.0;
1008 while(i < degreesOfFreedom)
1009 detail::noncentralChi2OneIteration(
arg, lans, dans, pans, i);
1010 degreesOfFreedom = degreesOfFreedom + 1;
1015 for(++m; m<maxit; ++m)
1018 detail::noncentralChi2OneIteration(
arg, lans, dans, pans, degreesOfFreedom);
1020 density = density + am * dans;
1022 probability = probability + hold;
1023 if((pans *
sum < eps2) && (hold < eps2))
1027 vigra_fail(
"noncentralChi2P(): no convergence.");
1028 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
1042inline double chi2(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1044 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0,
arg, accuracy).first;
1057inline double chi2CDF(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1059 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0,
arg, accuracy).second;
1074 double noncentrality,
double arg,
double accuracy = 1e-7)
1076 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality,
arg, accuracy).first;
1092 double noncentrality,
double arg,
double accuracy = 1e-7)
1094 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality,
arg, accuracy).second;
1113 return detail::noncentralChi2CDFApprox(degreesOfFreedom, noncentrality,
arg);
1123 T tmp = NumericTraits<T>::one();
1124 for(T f = l-m+1; f <= l+m; ++f)
1141template <
class REAL>
1144 vigra_precondition(
abs(x) <= 1.0,
"legendre(): x must be in [-1.0, 1.0].");
1151 return legendre(l,m,x) * s / detail::facLM<REAL>(l,m);
1156 REAL r = std::sqrt( (1.0-x) * (1.0+x) );
1158 for (
int i=1; i<=m; i++)
1167 REAL result_1 = x * (2.0 * m + 1.0) * result;
1171 for(
unsigned int i = m+2; i <= l; ++i)
1173 other = ( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m);
1188template <
class REAL>
1203template <
class REAL>
1209 return std::sin(M_PI * x);
1211 bool invert =
false;
1218 REAL rem = std::floor(x);
1225 rem = NumericTraits<REAL>::one();
1227 rem = std::sin(M_PI * rem);
1241template <
class REAL>
1249template <
class REAL>
1252 static REAL
gamma(REAL x);
1265template <
class REAL>
1266double GammaImpl<REAL>::g[] = {
1269 -0.6558780715202538,
1270 -0.420026350340952e-1,
1272 -0.421977345555443e-1,
1273 -0.9621971527877e-2,
1275 -0.11651675918591e-2,
1276 -0.2152416741149e-3,
1294template <
class REAL>
1295double GammaImpl<REAL>::a[] = {
1296 7.72156649015328655494e-02,
1297 3.22467033424113591611e-01,
1298 6.73523010531292681824e-02,
1299 2.05808084325167332806e-02,
1300 7.38555086081402883957e-03,
1301 2.89051383673415629091e-03,
1302 1.19270763183362067845e-03,
1303 5.10069792153511336608e-04,
1304 2.20862790713908385557e-04,
1305 1.08011567247583939954e-04,
1306 2.52144565451257326939e-05,
1307 4.48640949618915160150e-05
1310template <
class REAL>
1311double GammaImpl<REAL>::t[] = {
1312 4.83836122723810047042e-01,
1313 -1.47587722994593911752e-01,
1314 6.46249402391333854778e-02,
1315 -3.27885410759859649565e-02,
1316 1.79706750811820387126e-02,
1317 -1.03142241298341437450e-02,
1318 6.10053870246291332635e-03,
1319 -3.68452016781138256760e-03,
1320 2.25964780900612472250e-03,
1321 -1.40346469989232843813e-03,
1322 8.81081882437654011382e-04,
1323 -5.38595305356740546715e-04,
1324 3.15632070903625950361e-04,
1325 -3.12754168375120860518e-04,
1326 3.35529192635519073543e-04
1329template <
class REAL>
1330double GammaImpl<REAL>::u[] = {
1331 -7.72156649015328655494e-02,
1332 6.32827064025093366517e-01,
1333 1.45492250137234768737e+00,
1334 9.77717527963372745603e-01,
1335 2.28963728064692451092e-01,
1336 1.33810918536787660377e-02
1339template <
class REAL>
1340double GammaImpl<REAL>::v[] = {
1342 2.45597793713041134822e+00,
1343 2.12848976379893395361e+00,
1344 7.69285150456672783825e-01,
1345 1.04222645593369134254e-01,
1346 3.21709242282423911810e-03
1349template <
class REAL>
1350double GammaImpl<REAL>::s[] = {
1351 -7.72156649015328655494e-02,
1352 2.14982415960608852501e-01,
1353 3.25778796408930981787e-01,
1354 1.46350472652464452805e-01,
1355 2.66422703033638609560e-02,
1356 1.84028451407337715652e-03,
1357 3.19475326584100867617e-05
1360template <
class REAL>
1361double GammaImpl<REAL>::r[] = {
1363 1.39200533467621045958e+00,
1364 7.21935547567138069525e-01,
1365 1.71933865632803078993e-01,
1366 1.86459191715652901344e-02,
1367 7.77942496381893596434e-04,
1368 7.32668430744625636189e-06
1371template <
class REAL>
1372double GammaImpl<REAL>::w[] = {
1373 4.18938533204672725052e-01,
1374 8.33333333333329678849e-02,
1375 -2.77777777728775536470e-03,
1376 7.93650558643019558500e-04,
1377 -5.95187557450339963135e-04,
1378 8.36339918996282139126e-04,
1379 -1.63092934096575273989e-03
1382template <
class REAL>
1383REAL GammaImpl<REAL>::gamma(REAL x)
1385 int i, k, m, ix = (int)x;
1386 double ga = 0.0, gr = 0.0, r = 0.0, z = 0.0;
1388 vigra_precondition(x <= 171.0,
1389 "gamma(): argument cannot exceed 171.0.");
1396 for (i=2; i<ix; ++i)
1403 vigra_precondition(
false,
1404 "gamma(): gamma function is undefined for 0 and negative integers.");
1414 for (k=1; k<=m; ++k)
1425 for (k=23; k>=0; --k)
1435 ga = -M_PI/(x*ga*
sin_pi(x));
1455template <
class REAL>
1456REAL GammaImpl<REAL>::loggamma(REAL x)
1458 vigra_precondition(x > 0.0,
1459 "loggamma(): argument must be positive.");
1461 vigra_precondition(x <= 1.0e307,
1462 "loggamma(): argument must not exceed 1e307.");
1466 if (x < 4.2351647362715017e-22)
1470 else if ((x == 2.0) || (x == 1.0))
1476 const double tc = 1.46163214496836224576e+00;
1477 const double tf = -1.21486290535849611461e-01;
1478 const double tt = -3.63867699703950536541e-18;
1486 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1487 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1491 else if (x >= 0.23164)
1493 double y = x-(tc-1.0);
1496 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1497 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1498 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1499 double p = z*p1-(tt-w*(p2+y*p3));
1505 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1506 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1507 res += (-0.5*y + p1/p2);
1517 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1518 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1522 else if(x >= 1.23164)
1527 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1528 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1529 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1530 double p = z*p1-(tt-w*(p2+y*p3));
1536 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1537 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1538 res += (-0.5*y + p1/p2);
1544 double i = std::floor(x);
1546 double p = y*(s[0]+y*(s[1]+y*(s[2]+y*(s[3]+y*(s[4]+y*(s[5]+y*s[6]))))));
1547 double q = 1.0+y*(r[1]+y*(r[2]+y*(r[3]+y*(r[4]+y*(r[5]+y*r[6])))));
1557 else if (x < 2.8823037615171174e+17)
1559 double t = std::log(x);
1562 double yy = w[0]+z*(w[1]+y*(w[2]+y*(w[3]+y*(w[4]+y*(w[5]+y*w[6])))));
1563 res = (x-0.5)*(t-1.0)+yy;
1567 res = x*(std::log(x) - 1.0);
1589 return detail::GammaImpl<double>::gamma(x);
1605 return detail::GammaImpl<double>::loggamma(x);
1614FPT safeFloatDivision( FPT f1, FPT f2 )
1616 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
1617 ? NumericTraits<FPT>::max()
1618 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) ||
1619 f1 == NumericTraits<FPT>::zero()
1620 ? NumericTraits<FPT>::zero()
1636template <
class T1,
class T2>
1640 typedef typename PromoteTraits<T1, T2>::Promote T;
1642 return VIGRA_CSTD::fabs(r) <= epsilon;
1644 return VIGRA_CSTD::fabs(l) <= epsilon;
1645 T diff = VIGRA_CSTD::fabs( l - r );
1646 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
1647 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
1649 return (d1 <= epsilon && d2 <= epsilon);
1652template <
class T1,
class T2>
1655 typedef typename PromoteTraits<T1, T2>::Promote T;
1670template <
class T1,
class T2>
1677template <
class T1,
class T2>
1680 typedef typename PromoteTraits<T1, T2>::Promote T;
1695template <
class T1,
class T2>
1702template <
class T1,
class T2>
1705 typedef typename PromoteTraits<T1, T2>::Promote T;
1711#define VIGRA_MATH_FUNC_HELPER(TYPE) \
1712 inline TYPE clipLower(const TYPE t){ \
1713 return t < static_cast<TYPE>(0.0) ? static_cast<TYPE>(0.0) : t; \
1715 inline TYPE clipLower(const TYPE t,const TYPE valLow){ \
1716 return t < static_cast<TYPE>(valLow) ? static_cast<TYPE>(valLow) : t; \
1718 inline TYPE clipUpper(const TYPE t,const TYPE valHigh){ \
1719 return t > static_cast<TYPE>(valHigh) ? static_cast<TYPE>(valHigh) : t; \
1721 inline TYPE clip(const TYPE t,const TYPE valLow, const TYPE valHigh){ \
1724 else if(t>valHigh) \
1729 inline TYPE sum(const TYPE t){ \
1732 inline NumericTraits<TYPE>::RealPromote mean(const TYPE t){ \
1735 inline TYPE isZero(const TYPE t){ \
1736 return t==static_cast<TYPE>(0); \
1738 inline NumericTraits<TYPE>::RealPromote sizeDividedSquaredNorm(const TYPE t){ \
1739 return squaredNorm(t); \
1741 inline NumericTraits<TYPE>::RealPromote sizeDividedNorm(const TYPE t){ \
1747#pragma GCC diagnostic push
1748#pragma GCC diagnostic ignored "-Wtype-limits"
1751VIGRA_MATH_FUNC_HELPER(
unsigned char)
1752VIGRA_MATH_FUNC_HELPER(
unsigned short)
1753VIGRA_MATH_FUNC_HELPER(
unsigned int)
1754VIGRA_MATH_FUNC_HELPER(
unsigned long)
1755VIGRA_MATH_FUNC_HELPER(
unsigned long long)
1756VIGRA_MATH_FUNC_HELPER(
signed char)
1757VIGRA_MATH_FUNC_HELPER(
signed short)
1758VIGRA_MATH_FUNC_HELPER(
signed int)
1759VIGRA_MATH_FUNC_HELPER(
signed long)
1760VIGRA_MATH_FUNC_HELPER(
signed long long)
1761VIGRA_MATH_FUNC_HELPER(
float)
1762VIGRA_MATH_FUNC_HELPER(
double)
1763VIGRA_MATH_FUNC_HELPER(
long double)
1766#pragma GCC diagnostic pop
1769#undef VIGRA_MATH_FUNC_HELPER
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
double chi2(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Chi square distribution.
Definition: mathutil.hxx:1042
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition: sized_int.hxx:183
int round(FixedPoint< IntBits, FracBits > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:683
V power(const V &x)
Exponentiation to a positive integer power by squaring.
Definition: mathutil.hxx:427
double gamma(double x)
The gamma function.
Definition: mathutil.hxx:1587
R arg(const FFTWComplex< R > &a)
phase
Definition: fftw3.hxx:1009
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
double ellipticIntegralF(double x, double k)
The incomplete elliptic integral of the first kind.
Definition: mathutil.hxx:855
Int32 sqrti(Int32 v)
Signed integer square root.
Definition: mathutil.hxx:538
double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Cumulative chi square distribution.
Definition: mathutil.hxx:1057
REAL cos_pi(REAL x)
cos(pi*x).
Definition: mathutil.hxx:1242
REAL sin_pi(REAL x)
sin(pi*x).
Definition: mathutil.hxx:1204
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
Int32 roundi(FixedPoint16< IntBits, OverflowHandling > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:1775
UInt32 ceilPower2(UInt32 x)
Round up to the nearest power of 2.
Definition: mathutil.hxx:294
REAL legendre(unsigned int l, int m, REAL x)
Associated Legendre polynomial.
Definition: mathutil.hxx:1142
bool lessEqualAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point less-or-equal.
Definition: mathutil.hxx:1672
bool even(int t)
Check if an integer is even.
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22, T *r0, T *r1, T *r2)
Compute the eigenvalues of a 3x3 real symmetric matrix.
Definition: mathutil.hxx:754
UInt32 floorPower2(UInt32 x)
Round down to the nearest power of 2.
Definition: mathutil.hxx:317
Int32 log2i(UInt32 x)
Compute the base-2 logarithm of an integer.
Definition: mathutil.hxx:361
double loggamma(double x)
The natural logarithm of the gamma function.
Definition: mathutil.hxx:1603
T sign(T t)
The sign function.
Definition: mathutil.hxx:591
double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg)
Cumulative non-central chi square distribution (approximate).
Definition: mathutil.hxx:1111
void symmetric2x2Eigenvalues(T a00, T a01, T a11, T *r0, T *r1)
Compute the eigenvalues of a 2x2 real symmetric matrix.
Definition: mathutil.hxx:734
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
double noncentralChi2(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Non-central chi square distribution.
Definition: mathutil.hxx:1073
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point equality.
Definition: mathutil.hxx:1638
double ellipticIntegralE(double x, double k)
The incomplete elliptic integral of the second kind.
Definition: mathutil.hxx:879
bool greaterEqualAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point greater-or-equal.
Definition: mathutil.hxx:1697
detail::SelectIntegerType< 32, detail::SignedIntTypes >::type Int32
32-bit signed int
Definition: sized_int.hxx:175
int signi(T t)
The integer sign function.
Definition: mathutil.hxx:608
double noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Cumulative non-central chi square distribution.
Definition: mathutil.hxx:1091
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition: fixedpoint.hxx:1640
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
bool isPower2(UInt32 x)
Determine whether x is a power of 2 Bit twiddle from https://graphics.stanford.edu/~seander/bithacks....
Definition: mathutil.hxx:278
bool odd(int t)
Check if an integer is odd.