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()
596 ? -NumericTraits<T>::one()
597 : NumericTraits<T>::zero();
610 return 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
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)));
737 *
r0 =
static_cast<T
>(0.5*(
a00 +
a11 + d));
738 *
r1 =
static_cast<T
>(0.5*(
a00 +
a11 - d));
757 double inv3 = 1.0 / 3.0,
root3 = std::sqrt(3.0);
770 double magnitude = std::sqrt(-
aDiv3);
774 *
r0 =
static_cast<T
>(
c2Div3 + 2.0*magnitude*
cs);
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);
936double noncentralChi2CDFApprox(
unsigned int degreesOfFreedom, T noncentrality, T
arg)
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)));
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].");
1156 REAL r = std::sqrt( (1.0-x) * (1.0+x) );
1158 for (
int i=1;
i<=
m;
i++)
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);
1225 rem = NumericTraits<REAL>::one();
1227 rem = std::sin(M_PI *
rem);
1241template <
class REAL>
1244 return sin_pi(x+0.5);
1249template <
class REAL>
1252 static REAL gamma(REAL x);
1253 static REAL loggamma(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
Class for a single RGB value.
Definition rgbvalue.hxx:128
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
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
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.