36 #ifndef VIGRA_MATHUTIL_HXX
37 #define VIGRA_MATHUTIL_HXX
40 # pragma warning (disable: 4996) // hypot/_hypot confusion
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
109 #ifndef M_EULER_GAMMA
110 # define M_EULER_GAMMA 0.5772156649015329
122 using VIGRA_CSTD::pow;
134 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \
135 inline T abs(T t) { return t; }
137 VIGRA_DEFINE_UNSIGNED_ABS(
bool)
138 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned char)
139 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned short)
140 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned int)
141 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned long)
142 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned long long)
144 #undef VIGRA_DEFINE_UNSIGNED_ABS
146 #define VIGRA_DEFINE_MISSING_ABS(T) \
147 inline T abs(T t) { return t < 0 ? static_cast<T>(-t) : t; }
149 VIGRA_DEFINE_MISSING_ABS(
signed char)
150 VIGRA_DEFINE_MISSING_ABS(
signed short)
152 #if defined(_MSC_VER) && _MSC_VER < 1600
153 VIGRA_DEFINE_MISSING_ABS(
signed long long)
156 #undef VIGRA_DEFINE_MISSING_ABS
161 #define VIGRA_DEFINE_SCALAR_DOT(T) \
162 inline NumericTraits<T>::Promote dot(T l, T r) { return l*r; }
164 VIGRA_DEFINE_SCALAR_DOT(
unsigned char)
165 VIGRA_DEFINE_SCALAR_DOT(
unsigned short)
166 VIGRA_DEFINE_SCALAR_DOT(
unsigned int)
167 VIGRA_DEFINE_SCALAR_DOT(
unsigned long)
168 VIGRA_DEFINE_SCALAR_DOT(
unsigned long long)
169 VIGRA_DEFINE_SCALAR_DOT(
signed char)
170 VIGRA_DEFINE_SCALAR_DOT(
signed short)
171 VIGRA_DEFINE_SCALAR_DOT(
signed int)
172 VIGRA_DEFINE_SCALAR_DOT(
signed long)
173 VIGRA_DEFINE_SCALAR_DOT(
signed long long)
174 VIGRA_DEFINE_SCALAR_DOT(
float)
175 VIGRA_DEFINE_SCALAR_DOT(
double)
176 VIGRA_DEFINE_SCALAR_DOT(
long double)
178 #undef VIGRA_DEFINE_SCALAR_DOT
184 inline float pow(
float v,
double e)
186 return std::pow(v, (
float)e);
189 inline long double pow(
long double v,
double e)
191 return std::pow(v, (
long double)e);
202 #ifdef DOXYGEN // only for documentation
206 inline float round(
float t)
213 inline double round(
double t)
220 inline long double round(
long double t)
291 static Int32 table[64];
295 Int32 IntLog2<T>::table[64] = {
296 -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21,
297 29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6,
298 -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9,
299 11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25,
300 -1, 7, 24, -1, 23, -1, 31, -1};
329 return detail::IntLog2<Int32>::table[x >> 26];
341 typename NumericTraits<T>::Promote
sq(T t)
348 template <
class V,
unsigned>
351 static V call(
const V & x,
const V & y) {
return x * y; }
354 struct cond_mult<V, 0>
356 static V call(
const V &,
const V & y) {
return y; }
359 template <
class V,
unsigned n>
362 static V call(
const V & x)
365 ? cond_mult<V, n & 1>::call(x, power_static<V, n / 2>::call(x * x))
370 struct power_static<V, 0>
372 static V call(
const V & x)
385 template <
unsigned n,
class V>
388 return detail::power_static<V, n>::call(x);
398 static UInt32 sqq_table[];
403 UInt32 IntSquareRoot<T>::sqq_table[] = {
404 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57,
405 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83,
406 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102,
407 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
408 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
409 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
410 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
411 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
412 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
413 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
414 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
415 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
416 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
417 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
418 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
419 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
420 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
421 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
432 if (x >= 0x40000000) {
435 xn = sqq_table[x>>24] << 8;
437 xn = sqq_table[x>>22] << 7;
440 xn = sqq_table[x>>20] << 6;
442 xn = sqq_table[x>>18] << 5;
446 xn = sqq_table[x>>16] << 4;
448 xn = sqq_table[x>>14] << 3;
451 xn = sqq_table[x>>12] << 2;
453 xn = sqq_table[x>>10] << 1;
461 xn = (sqq_table[x>>8] >> 0) + 1;
463 xn = (sqq_table[x>>6] >> 1) + 1;
466 xn = (sqq_table[x>>4] >> 2) + 1;
468 xn = (sqq_table[x>>2] >> 3) + 1;
472 return sqq_table[x] >> 4;
476 xn = (xn + 1 + x / xn) / 2;
478 xn = (xn + 1 + x / xn) / 2;
501 throw std::domain_error(
"sqrti(Int32): negative argument.");
502 return (
Int32)detail::IntSquareRoot<UInt32>::exec((
UInt32)v);
514 return detail::IntSquareRoot<UInt32>::exec(v);
517 #ifdef VIGRA_NO_HYPOT
526 inline double hypot(
double a,
double b)
528 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b);
553 return t > NumericTraits<T>::zero()
554 ? NumericTraits<T>::one()
555 : t < NumericTraits<T>::zero()
556 ? -NumericTraits<T>::one()
557 : NumericTraits<T>::zero();
570 return t > NumericTraits<T>::zero()
572 : t < NumericTraits<T>::zero()
584 template <
class T1,
class T2>
587 return t2 >= NumericTraits<T2>::zero()
593 #ifdef DOXYGEN // only for documentation
608 #define VIGRA_DEFINE_ODD_EVEN(T) \
609 inline bool even(T t) { return (t&1) == 0; } \
610 inline bool odd(T t) { return (t&1) == 1; }
612 VIGRA_DEFINE_ODD_EVEN(
char)
613 VIGRA_DEFINE_ODD_EVEN(
short)
614 VIGRA_DEFINE_ODD_EVEN(
int)
615 VIGRA_DEFINE_ODD_EVEN(
long)
616 VIGRA_DEFINE_ODD_EVEN(
long long)
617 VIGRA_DEFINE_ODD_EVEN(
unsigned char)
618 VIGRA_DEFINE_ODD_EVEN(
unsigned short)
619 VIGRA_DEFINE_ODD_EVEN(
unsigned int)
620 VIGRA_DEFINE_ODD_EVEN(
unsigned long)
621 VIGRA_DEFINE_ODD_EVEN(
unsigned long long)
623 #undef VIGRA_DEFINE_ODD_EVEN
625 #define VIGRA_DEFINE_NORM(T) \
626 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
627 inline NormTraits<T>::NormType norm(T t) { return abs(t); }
629 VIGRA_DEFINE_NORM(
bool)
630 VIGRA_DEFINE_NORM(
signed char)
631 VIGRA_DEFINE_NORM(
unsigned char)
632 VIGRA_DEFINE_NORM(
short)
633 VIGRA_DEFINE_NORM(
unsigned short)
634 VIGRA_DEFINE_NORM(
int)
635 VIGRA_DEFINE_NORM(
unsigned int)
636 VIGRA_DEFINE_NORM(
long)
637 VIGRA_DEFINE_NORM(
unsigned long)
638 VIGRA_DEFINE_NORM(
long long)
639 VIGRA_DEFINE_NORM(
unsigned long long)
640 VIGRA_DEFINE_NORM(
float)
641 VIGRA_DEFINE_NORM(
double)
642 VIGRA_DEFINE_NORM(
long double)
644 #undef VIGRA_DEFINE_NORM
647 inline typename NormTraits<std::complex<T> >::SquaredNormType
650 return sq(t.real()) +
sq(t.imag());
653 #ifdef DOXYGEN // only for documentation
661 NormTraits<T>::SquaredNormType
squaredNorm(T
const & t);
674 inline typename NormTraits<T>::NormType
677 typedef typename NormTraits<T>::SquaredNormType SNT;
678 return sqrt(
static_cast<typename SquareRootTraits<SNT>::SquareRootArgument
>(
squaredNorm(t)));
694 double d =
hypot(a00 - a11, 2.0*a01);
695 *r0 =
static_cast<T
>(0.5*(a00 + a11 + d));
696 *r1 =
static_cast<T
>(0.5*(a00 + a11 - d));
714 T * r0, T * r1, T * r2)
716 static double inv3 = 1.0 / 3.0, root3 =
std::sqrt(3.0);
718 double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
719 double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
720 double c2 = a00 + a11 + a22;
721 double c2Div3 = c2*inv3;
722 double aDiv3 = (c1 - c2*c2Div3)*inv3;
725 double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
726 double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
733 *r0 =
static_cast<T
>(c2Div3 + 2.0*magnitude*cs);
734 *r1 =
static_cast<T
>(c2Div3 - magnitude*(cs + root3*sn));
735 *r2 =
static_cast<T
>(c2Div3 - magnitude*(cs - root3*sn));
747 T ellipticRD(T x, T y, T z)
749 double f = 1.0, s = 0.0, X, Y, Z, m;
752 m = (x + y + 3.0*z) / 5.0;
756 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
768 double d = a - 6.0*b;
769 double e = d + 2.0*c;
770 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0)
771 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5);
775 T ellipticRF(T x, T y, T z)
780 m = (x + y + z) / 3.0;
784 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
791 double d = X*Y -
sq(Z);
818 return s*detail::ellipticRF(c2, 1.0 -
sq(k*s), 1.0);
843 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0));
853 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
854 double ans = t*
VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
855 t*(0.09678418+t*(-0.18628806+t*(0.27886807+
856 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
857 t*0.17087277)))))))));
881 inline double erf(
double x)
883 return detail::erfImpl(x);
897 double a = degreesOfFreedom + noncentrality;
898 double b = (a + noncentrality) /
sq(a);
899 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);
904 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans,
unsigned int & j)
914 dans = dans * arg / j;
921 std::pair<double, double>
noncentralChi2CDF(
unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
923 vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,
924 "noncentralChi2P(): parameters must be positive.");
925 if (arg == 0.0 && degreesOfFreedom > 0)
926 return std::make_pair(0.0, 0.0);
929 double b1 = 0.5 * noncentrality,
932 lnrtpi2 = 0.22579135264473,
933 probability, density, lans, dans, pans,
sum, am, hold;
934 unsigned int maxit = 500,
936 if(degreesOfFreedom % 2)
952 if(degreesOfFreedom == 0)
955 degreesOfFreedom = 2;
957 sum = 1.0 / ao - 1.0 - am;
959 probability = 1.0 + am * pans;
964 degreesOfFreedom = degreesOfFreedom - 1;
966 sum = 1.0 / ao - 1.0;
967 while(i < degreesOfFreedom)
968 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
969 degreesOfFreedom = degreesOfFreedom + 1;
974 for(++m; m<maxit; ++m)
977 detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
979 density = density + am * dans;
981 probability = probability + hold;
982 if((pans * sum < eps2) && (hold < eps2))
986 vigra_fail(
"noncentralChi2P(): no convergence.");
987 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
1001 inline double chi2(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1016 inline double chi2CDF(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1033 double noncentrality,
double arg,
double accuracy = 1e-7)
1051 double noncentrality,
double arg,
double accuracy = 1e-7)
1082 T tmp = NumericTraits<T>::one();
1083 for(T f = l-m+1; f <= l+m; ++f)
1100 template <
class REAL>
1103 vigra_precondition(
abs(x) <= 1.0,
"legendre(): x must be in [-1.0, 1.0].");
1110 return legendre(l,m,x) * s / detail::facLM<REAL>(l,m);
1115 REAL r =
std::sqrt( (1.0-x) * (1.0+x) );
1117 for (
int i=1; i<=m; i++)
1126 REAL result_1 = x * (2.0 * m + 1.0) * result;
1130 for(
unsigned int i = m+2; i <= l; ++i)
1132 other = ( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m);
1147 template <
class REAL>
1162 template <
class REAL>
1170 bool invert =
false;
1184 rem = NumericTraits<REAL>::one();
1200 template <
class REAL>
1208 template <
class REAL>
1209 REAL gammaImpl(REAL x)
1211 int i, k, m, ix = (int)x;
1212 double ga = 0.0, gr = 0.0, r = 0.0, z = 0.0;
1214 static double g[] = {
1217 -0.6558780715202538,
1218 -0.420026350340952e-1,
1220 -0.421977345555443e-1,
1221 -0.9621971527877e-2,
1223 -0.11651675918591e-2,
1224 -0.2152416741149e-3,
1241 vigra_precondition(x <= 171.0,
1242 "gamma(): argument cannot exceed 171.0.");
1249 for (i=2; i<ix; ++i)
1256 vigra_precondition(
false,
1257 "gamma(): gamma function is undefined for 0 and negative integers.");
1267 for (k=1; k<=m; ++k)
1278 for (k=23; k>=0; --k)
1288 ga = -M_PI/(x*ga*
sin_pi(x));
1308 template <
class REAL>
1309 REAL loggammaImpl(REAL x)
1311 vigra_precondition(x > 0.0,
1312 "loggamma(): argument must be positive.");
1314 vigra_precondition(x <= 1.0e307,
1315 "loggamma(): argument must not exceed 1e307.");
1319 if (x < 4.2351647362715017e-22)
1323 else if ((x == 2.0) || (x == 1.0))
1329 static const double a[] = { 7.72156649015328655494e-02,
1330 3.22467033424113591611e-01,
1331 6.73523010531292681824e-02,
1332 2.05808084325167332806e-02,
1333 7.38555086081402883957e-03,
1334 2.89051383673415629091e-03,
1335 1.19270763183362067845e-03,
1336 5.10069792153511336608e-04,
1337 2.20862790713908385557e-04,
1338 1.08011567247583939954e-04,
1339 2.52144565451257326939e-05,
1340 4.48640949618915160150e-05 };
1341 static const double t[] = { 4.83836122723810047042e-01,
1342 -1.47587722994593911752e-01,
1343 6.46249402391333854778e-02,
1344 -3.27885410759859649565e-02,
1345 1.79706750811820387126e-02,
1346 -1.03142241298341437450e-02,
1347 6.10053870246291332635e-03,
1348 -3.68452016781138256760e-03,
1349 2.25964780900612472250e-03,
1350 -1.40346469989232843813e-03,
1351 8.81081882437654011382e-04,
1352 -5.38595305356740546715e-04,
1353 3.15632070903625950361e-04,
1354 -3.12754168375120860518e-04,
1355 3.35529192635519073543e-04 };
1356 static const double u[] = { -7.72156649015328655494e-02,
1357 6.32827064025093366517e-01,
1358 1.45492250137234768737e+00,
1359 9.77717527963372745603e-01,
1360 2.28963728064692451092e-01,
1361 1.33810918536787660377e-02 };
1362 static const double v[] = { 0.0,
1363 2.45597793713041134822e+00,
1364 2.12848976379893395361e+00,
1365 7.69285150456672783825e-01,
1366 1.04222645593369134254e-01,
1367 3.21709242282423911810e-03 };
1368 static const double tc = 1.46163214496836224576e+00;
1369 static const double tf = -1.21486290535849611461e-01;
1370 static const double tt = -3.63867699703950536541e-18;
1378 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1379 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1383 else if (x >= 0.23164)
1385 double y = x-(tc-1.0);
1388 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1389 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1390 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1391 double p = z*p1-(tt-w*(p2+y*p3));
1397 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1398 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1399 res += (-0.5*y + p1/p2);
1409 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1410 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1414 else if(x >= 1.23164)
1419 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1420 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1421 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1422 double p = z*p1-(tt-w*(p2+y*p3));
1428 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1429 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1430 res += (-0.5*y + p1/p2);
1436 static const double s[] = { -7.72156649015328655494e-02,
1437 2.14982415960608852501e-01,
1438 3.25778796408930981787e-01,
1439 1.46350472652464452805e-01,
1440 2.66422703033638609560e-02,
1441 1.84028451407337715652e-03,
1442 3.19475326584100867617e-05 };
1443 static const double r[] = { 0.0,
1444 1.39200533467621045958e+00,
1445 7.21935547567138069525e-01,
1446 1.71933865632803078993e-01,
1447 1.86459191715652901344e-02,
1448 7.77942496381893596434e-04,
1449 7.32668430744625636189e-06 };
1452 double p = y*(s[0]+y*(s[1]+y*(s[2]+y*(s[3]+y*(s[4]+y*(s[5]+y*s[6]))))));
1453 double q = 1.0+y*(r[1]+y*(r[2]+y*(r[3]+y*(r[4]+y*(r[5]+y*r[6])))));
1463 else if (x < 2.8823037615171174e+17)
1465 static const double w[] = { 4.18938533204672725052e-01,
1466 8.33333333333329678849e-02,
1467 -2.77777777728775536470e-03,
1468 7.93650558643019558500e-04,
1469 -5.95187557450339963135e-04,
1470 8.36339918996282139126e-04,
1471 -1.63092934096575273989e-03 };
1475 double yy = w[0]+z*(w[1]+y*(w[2]+y*(w[3]+y*(w[4]+y*(w[5]+y*w[6])))));
1476 res = (x-0.5)*(t-1.0)+yy;
1502 return detail::gammaImpl(x);
1518 return detail::loggammaImpl(x);
1527 FPT safeFloatDivision( FPT f1, FPT f2 )
1529 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
1530 ? NumericTraits<FPT>::max()
1531 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) ||
1532 f1 == NumericTraits<FPT>::zero()
1533 ? NumericTraits<FPT>::zero()
1549 template <
class T1,
class T2>
1553 typedef typename PromoteTraits<T1, T2>::Promote T;
1555 return VIGRA_CSTD::fabs(r) <= epsilon;
1557 return VIGRA_CSTD::fabs(l) <= epsilon;
1558 T diff = VIGRA_CSTD::fabs( l - r );
1559 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
1560 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
1562 return (d1 <= epsilon && d2 <= epsilon);
1565 template <
class T1,
class T2>
1568 typedef typename PromoteTraits<T1, T2>::Promote T;
double ellipticIntegralF(double x, double k)
Definition: mathutil.hxx:814
REAL sin_pi(REAL x)
Definition: mathutil.hxx:1163
FixedPoint16< 2, OverflowHandling > atan2(FixedPoint16< IntBits, OverflowHandling > y, FixedPoint16< IntBits, OverflowHandling > x)
Arctangent. Accuracy better than 1/3 degree (9 significant bits).
Definition: fixedpoint.hxx:1640
R arg(const FFTWComplex< R > &a)
pahse
Definition: fftw3.hxx:1009
Int32 roundi(FixedPoint16< IntBits, OverflowHandling > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:1761
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22, T *r0, T *r1, T *r2)
Definition: mathutil.hxx:713
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
V power(const V &x)
Definition: mathutil.hxx:386
int round(FixedPoint< IntBits, FracBits > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:683
REAL legendre(unsigned int l, int m, REAL x)
Definition: mathutil.hxx:1101
int signi(T t)
Definition: mathutil.hxx:568
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition: fixedpoint.hxx:1626
REAL cos_pi(REAL x)
Definition: mathutil.hxx:1201
double noncentralChi2(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Definition: mathutil.hxx:1032
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg)
Definition: mathutil.hxx:1070
Int32 log2i(UInt32 x)
Definition: mathutil.hxx:320
NumericTraits< T >::Promote sq(T t)
Definition: mathutil.hxx:341
double noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Definition: mathutil.hxx:1050
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:1683
detail::SelectIntegerType< 32, detail::SignedIntTypes >::type Int32
32-bit signed int
Definition: sized_int.hxx:175
void symmetric2x2Eigenvalues(T a00, T a01, T a11, T *r0, T *r1)
Definition: mathutil.hxx:692
double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Definition: mathutil.hxx:1016
double gamma(double x)
Definition: mathutil.hxx:1500
bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Definition: mathutil.hxx:1551
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
double loggamma(double x)
Definition: mathutil.hxx:1516
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition: sized_int.hxx:183
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
Int32 sqrti(Int32 v)
Definition: mathutil.hxx:498
UInt32 floorPower2(UInt32 x)
Definition: mathutil.hxx:275
double chi2(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Definition: mathutil.hxx:1001
UInt32 ceilPower2(UInt32 x)
Definition: mathutil.hxx:253
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
T sign(T t)
Definition: mathutil.hxx:551
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
double ellipticIntegralE(double x, double k)
Definition: mathutil.hxx:838
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616