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 110 # define M_E 2.71828182845904523536 113 #ifndef M_EULER_GAMMA 114 # define M_EULER_GAMMA 0.5772156649015329 126 using VIGRA_CSTD::pow;
138 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \ 139 inline T abs(T t) { return t; } 141 VIGRA_DEFINE_UNSIGNED_ABS(
bool)
142 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned char)
143 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned short)
144 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned int)
145 VIGRA_DEFINE_UNSIGNED_ABS(
unsigned long)
146 VIGRA_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; } 153 VIGRA_DEFINE_MISSING_ABS(
signed char)
154 VIGRA_DEFINE_MISSING_ABS(
signed short)
156 #if defined(_MSC_VER) && _MSC_VER < 1600 157 VIGRA_DEFINE_MISSING_ABS(
signed long long)
160 #undef VIGRA_DEFINE_MISSING_ABS 165 #define VIGRA_DEFINE_SCALAR_DOT(T) \ 166 inline NumericTraits<T>::Promote dot(T l, T r) { return l*r; } 168 VIGRA_DEFINE_SCALAR_DOT(
unsigned char)
169 VIGRA_DEFINE_SCALAR_DOT(
unsigned short)
170 VIGRA_DEFINE_SCALAR_DOT(
unsigned int)
171 VIGRA_DEFINE_SCALAR_DOT(
unsigned long)
172 VIGRA_DEFINE_SCALAR_DOT(
unsigned long long)
173 VIGRA_DEFINE_SCALAR_DOT(
signed char)
174 VIGRA_DEFINE_SCALAR_DOT(
signed short)
175 VIGRA_DEFINE_SCALAR_DOT(
signed int)
176 VIGRA_DEFINE_SCALAR_DOT(
signed long)
177 VIGRA_DEFINE_SCALAR_DOT(
signed long long)
178 VIGRA_DEFINE_SCALAR_DOT(
float)
179 VIGRA_DEFINE_SCALAR_DOT(
double)
180 VIGRA_DEFINE_SCALAR_DOT(
long double)
182 #undef VIGRA_DEFINE_SCALAR_DOT 188 inline float pow(
float v,
double e)
190 return std::pow(v, (
float)e);
193 inline long double pow(
long double v,
double e)
195 return std::pow(v, (
long double)e);
206 #ifdef DOXYGEN // only for documentation 210 inline float round(
float t)
217 inline double round(
double t)
224 inline long double round(
long double t)
294 static Int32 table[64];
298 Int32 IntLog2<T>::table[64] = {
299 -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21,
300 29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6,
301 -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9,
302 11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25,
303 -1, 7, 24, -1, 23, -1, 31, -1};
332 return detail::IntLog2<Int32>::table[x >> 26];
344 typename NumericTraits<T>::Promote
sq(T t)
351 template <
class V,
unsigned>
354 static V call(
const V & x,
const V & y) {
return x * y; }
357 struct cond_mult<V, 0>
359 static V call(
const V &,
const V & y) {
return y; }
362 template <
class V,
unsigned n>
365 static V call(
const V & x)
368 ? cond_mult<V, n & 1>::call(x, power_static<V, n / 2>::call(x * x))
373 struct power_static<V, 0>
375 static V call(
const V & x)
388 template <
unsigned n,
class V>
391 return detail::power_static<V, n>::call(x);
400 static UInt32 sqq_table[];
405 UInt32 IntSquareRoot<T>::sqq_table[] = {
406 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57,
407 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83,
408 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102,
409 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
410 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
411 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
412 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
413 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
414 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
415 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
416 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
417 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
418 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
419 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
420 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
421 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
422 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
423 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
434 if (x >= 0x40000000) {
437 xn = sqq_table[x>>24] << 8;
439 xn = sqq_table[x>>22] << 7;
442 xn = sqq_table[x>>20] << 6;
444 xn = sqq_table[x>>18] << 5;
448 xn = sqq_table[x>>16] << 4;
450 xn = sqq_table[x>>14] << 3;
453 xn = sqq_table[x>>12] << 2;
455 xn = sqq_table[x>>10] << 1;
463 xn = (sqq_table[x>>8] >> 0) + 1;
465 xn = (sqq_table[x>>6] >> 1) + 1;
468 xn = (sqq_table[x>>4] >> 2) + 1;
470 xn = (sqq_table[x>>2] >> 3) + 1;
474 return sqq_table[x] >> 4;
478 xn = (xn + 1 + x / xn) / 2;
480 xn = (xn + 1 + x / xn) / 2;
503 throw std::domain_error(
"sqrti(Int32): negative argument.");
504 return (
Int32)detail::IntSquareRoot<UInt32>::exec((
UInt32)v);
516 return detail::IntSquareRoot<UInt32>::exec(v);
519 #ifdef VIGRA_NO_HYPOT 528 inline double hypot(
double a,
double b)
530 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b);
555 return t > NumericTraits<T>::zero()
556 ? NumericTraits<T>::one()
557 : t < NumericTraits<T>::zero()
558 ? -NumericTraits<T>::one()
559 : NumericTraits<T>::zero();
572 return t > NumericTraits<T>::zero()
574 : t < NumericTraits<T>::zero()
586 template <
class T1,
class T2>
589 return t2 >= NumericTraits<T2>::zero()
595 #ifdef DOXYGEN // only for documentation 610 #define VIGRA_DEFINE_ODD_EVEN(T) \ 611 inline bool even(T t) { return (t&1) == 0; } \ 612 inline bool odd(T t) { return (t&1) == 1; } 614 VIGRA_DEFINE_ODD_EVEN(
char)
615 VIGRA_DEFINE_ODD_EVEN(
short)
616 VIGRA_DEFINE_ODD_EVEN(
int)
617 VIGRA_DEFINE_ODD_EVEN(
long)
618 VIGRA_DEFINE_ODD_EVEN(
long long)
619 VIGRA_DEFINE_ODD_EVEN(
unsigned char)
620 VIGRA_DEFINE_ODD_EVEN(
unsigned short)
621 VIGRA_DEFINE_ODD_EVEN(
unsigned int)
622 VIGRA_DEFINE_ODD_EVEN(
unsigned long)
623 VIGRA_DEFINE_ODD_EVEN(
unsigned long long)
625 #undef VIGRA_DEFINE_ODD_EVEN 627 #define VIGRA_DEFINE_NORM(T) \ 628 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \ 629 inline NormTraits<T>::NormType norm(T t) { return abs(t); } 631 VIGRA_DEFINE_NORM(
bool)
632 VIGRA_DEFINE_NORM(
signed char)
633 VIGRA_DEFINE_NORM(
unsigned char)
634 VIGRA_DEFINE_NORM(
short)
635 VIGRA_DEFINE_NORM(
unsigned short)
636 VIGRA_DEFINE_NORM(
int)
637 VIGRA_DEFINE_NORM(
unsigned int)
638 VIGRA_DEFINE_NORM(
long)
639 VIGRA_DEFINE_NORM(
unsigned long)
640 VIGRA_DEFINE_NORM(
long long)
641 VIGRA_DEFINE_NORM(
unsigned long long)
642 VIGRA_DEFINE_NORM(
float)
643 VIGRA_DEFINE_NORM(
double)
644 VIGRA_DEFINE_NORM(
long double)
646 #undef VIGRA_DEFINE_NORM 649 inline typename NormTraits<std::complex<T> >::SquaredNormType
652 return sq(t.real()) +
sq(t.imag());
655 #ifdef DOXYGEN // only for documentation 665 NormTraits<T>::SquaredNormType
squaredNorm(T
const & t);
678 inline typename NormTraits<T>::NormType
681 typedef typename NormTraits<T>::SquaredNormType SNT;
682 return sqrt(
static_cast<typename SquareRootTraits<SNT>::SquareRootArgument
>(
squaredNorm(t)));
698 double d =
hypot(a00 - a11, 2.0*a01);
699 *r0 =
static_cast<T
>(0.5*(a00 + a11 + d));
700 *r1 =
static_cast<T
>(0.5*(a00 + a11 - d));
717 T * r0, T * r1, T * r2)
719 double inv3 = 1.0 / 3.0, root3 =
std::sqrt(3.0);
721 double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
722 double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
723 double c2 = a00 + a11 + a22;
724 double c2Div3 = c2*inv3;
725 double aDiv3 = (c1 - c2*c2Div3)*inv3;
728 double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
729 double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
736 *r0 =
static_cast<T
>(c2Div3 + 2.0*magnitude*cs);
737 *r1 =
static_cast<T
>(c2Div3 - magnitude*(cs + root3*sn));
738 *r2 =
static_cast<T
>(c2Div3 - magnitude*(cs - root3*sn));
750 T ellipticRD(T x, T y, T z)
752 double f = 1.0, s = 0.0, X, Y, Z, m;
755 m = (x + y + 3.0*z) / 5.0;
759 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
771 double d = a - 6.0*b;
772 double e = d + 2.0*c;
773 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0)
774 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5);
778 T ellipticRF(T x, T y, T z)
783 m = (x + y + z) / 3.0;
787 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
794 double d = X*Y -
sq(Z);
821 return s*detail::ellipticRF(c2, 1.0 -
sq(k*s), 1.0);
846 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0));
856 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
857 double ans = t*
VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
858 t*(0.09678418+t*(-0.18628806+t*(0.27886807+
859 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
860 t*0.17087277)))))))));
884 inline double erf(
double x)
886 return detail::erfImpl(x);
900 double a = degreesOfFreedom + noncentrality;
901 double b = (a + noncentrality) /
sq(a);
902 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);
907 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans,
unsigned int & j)
917 dans = dans * arg / j;
924 std::pair<double, double>
noncentralChi2CDF(
unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
926 vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,
927 "noncentralChi2P(): parameters must be positive.");
928 if (arg == 0.0 && degreesOfFreedom > 0)
929 return std::make_pair(0.0, 0.0);
932 double b1 = 0.5 * noncentrality,
935 lnrtpi2 = 0.22579135264473,
936 probability, density, lans, dans, pans,
sum, am, hold;
937 unsigned int maxit = 500,
939 if(degreesOfFreedom % 2)
955 if(degreesOfFreedom == 0)
958 degreesOfFreedom = 2;
960 sum = 1.0 / ao - 1.0 - am;
962 probability = 1.0 + am * pans;
967 degreesOfFreedom = degreesOfFreedom - 1;
969 sum = 1.0 / ao - 1.0;
970 while(i < degreesOfFreedom)
971 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
972 degreesOfFreedom = degreesOfFreedom + 1;
977 for(++m; m<maxit; ++m)
980 detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
982 density = density + am * dans;
984 probability = probability + hold;
985 if((pans * sum < eps2) && (hold < eps2))
989 vigra_fail(
"noncentralChi2P(): no convergence.");
990 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
1004 inline double chi2(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1019 inline double chi2CDF(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1036 double noncentrality,
double arg,
double accuracy = 1e-7)
1054 double noncentrality,
double arg,
double accuracy = 1e-7)
1085 T tmp = NumericTraits<T>::one();
1086 for(T f = l-m+1; f <= l+m; ++f)
1103 template <
class REAL>
1106 vigra_precondition(
abs(x) <= 1.0,
"legendre(): x must be in [-1.0, 1.0].");
1113 return legendre(l,m,x) * s / detail::facLM<REAL>(l,m);
1118 REAL r =
std::sqrt( (1.0-x) * (1.0+x) );
1120 for (
int i=1; i<=m; i++)
1129 REAL result_1 = x * (2.0 * m + 1.0) * result;
1133 for(
unsigned int i = m+2; i <= l; ++i)
1135 other = ( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m);
1150 template <
class REAL>
1165 template <
class REAL>
1173 bool invert =
false;
1187 rem = NumericTraits<REAL>::one();
1203 template <
class REAL>
1211 template <
class REAL>
1214 static REAL
gamma(REAL x);
1227 template <
class REAL>
1228 double GammaImpl<REAL>::g[] = {
1231 -0.6558780715202538,
1232 -0.420026350340952e-1,
1234 -0.421977345555443e-1,
1235 -0.9621971527877e-2,
1237 -0.11651675918591e-2,
1238 -0.2152416741149e-3,
1256 template <
class REAL>
1257 double GammaImpl<REAL>::a[] = {
1258 7.72156649015328655494e-02,
1259 3.22467033424113591611e-01,
1260 6.73523010531292681824e-02,
1261 2.05808084325167332806e-02,
1262 7.38555086081402883957e-03,
1263 2.89051383673415629091e-03,
1264 1.19270763183362067845e-03,
1265 5.10069792153511336608e-04,
1266 2.20862790713908385557e-04,
1267 1.08011567247583939954e-04,
1268 2.52144565451257326939e-05,
1269 4.48640949618915160150e-05
1272 template <
class REAL>
1273 double GammaImpl<REAL>::t[] = {
1274 4.83836122723810047042e-01,
1275 -1.47587722994593911752e-01,
1276 6.46249402391333854778e-02,
1277 -3.27885410759859649565e-02,
1278 1.79706750811820387126e-02,
1279 -1.03142241298341437450e-02,
1280 6.10053870246291332635e-03,
1281 -3.68452016781138256760e-03,
1282 2.25964780900612472250e-03,
1283 -1.40346469989232843813e-03,
1284 8.81081882437654011382e-04,
1285 -5.38595305356740546715e-04,
1286 3.15632070903625950361e-04,
1287 -3.12754168375120860518e-04,
1288 3.35529192635519073543e-04
1291 template <
class REAL>
1292 double GammaImpl<REAL>::u[] = {
1293 -7.72156649015328655494e-02,
1294 6.32827064025093366517e-01,
1295 1.45492250137234768737e+00,
1296 9.77717527963372745603e-01,
1297 2.28963728064692451092e-01,
1298 1.33810918536787660377e-02
1301 template <
class REAL>
1302 double GammaImpl<REAL>::v[] = {
1304 2.45597793713041134822e+00,
1305 2.12848976379893395361e+00,
1306 7.69285150456672783825e-01,
1307 1.04222645593369134254e-01,
1308 3.21709242282423911810e-03
1311 template <
class REAL>
1312 double GammaImpl<REAL>::s[] = {
1313 -7.72156649015328655494e-02,
1314 2.14982415960608852501e-01,
1315 3.25778796408930981787e-01,
1316 1.46350472652464452805e-01,
1317 2.66422703033638609560e-02,
1318 1.84028451407337715652e-03,
1319 3.19475326584100867617e-05
1322 template <
class REAL>
1323 double GammaImpl<REAL>::r[] = {
1325 1.39200533467621045958e+00,
1326 7.21935547567138069525e-01,
1327 1.71933865632803078993e-01,
1328 1.86459191715652901344e-02,
1329 7.77942496381893596434e-04,
1330 7.32668430744625636189e-06
1333 template <
class REAL>
1334 double GammaImpl<REAL>::w[] = {
1335 4.18938533204672725052e-01,
1336 8.33333333333329678849e-02,
1337 -2.77777777728775536470e-03,
1338 7.93650558643019558500e-04,
1339 -5.95187557450339963135e-04,
1340 8.36339918996282139126e-04,
1341 -1.63092934096575273989e-03
1344 template <
class REAL>
1347 int i, k, m, ix = (int)x;
1348 double ga = 0.0, gr = 0.0, r = 0.0, z = 0.0;
1350 vigra_precondition(x <= 171.0,
1351 "gamma(): argument cannot exceed 171.0.");
1358 for (i=2; i<ix; ++i)
1365 vigra_precondition(
false,
1366 "gamma(): gamma function is undefined for 0 and negative integers.");
1376 for (k=1; k<=m; ++k)
1387 for (k=23; k>=0; --k)
1397 ga = -M_PI/(x*ga*
sin_pi(x));
1417 template <
class REAL>
1420 vigra_precondition(x > 0.0,
1421 "loggamma(): argument must be positive.");
1423 vigra_precondition(x <= 1.0e307,
1424 "loggamma(): argument must not exceed 1e307.");
1428 if (x < 4.2351647362715017e-22)
1432 else if ((x == 2.0) || (x == 1.0))
1438 const double tc = 1.46163214496836224576e+00;
1439 const double tf = -1.21486290535849611461e-01;
1440 const double tt = -3.63867699703950536541e-18;
1448 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1449 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1453 else if (x >= 0.23164)
1455 double y = x-(tc-1.0);
1458 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1459 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1460 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1461 double p = z*p1-(tt-w*(p2+y*p3));
1467 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1468 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1469 res += (-0.5*y + p1/p2);
1479 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1480 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1484 else if(x >= 1.23164)
1489 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1490 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1491 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1492 double p = z*p1-(tt-w*(p2+y*p3));
1498 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1499 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1500 res += (-0.5*y + p1/p2);
1508 double p = y*(s[0]+y*(s[1]+y*(s[2]+y*(s[3]+y*(s[4]+y*(s[5]+y*s[6]))))));
1509 double q = 1.0+y*(r[1]+y*(r[2]+y*(r[3]+y*(r[4]+y*(r[5]+y*r[6])))));
1519 else if (x < 2.8823037615171174e+17)
1524 double yy = w[0]+z*(w[1]+y*(w[2]+y*(w[3]+y*(w[4]+y*(w[5]+y*w[6])))));
1525 res = (x-0.5)*(t-1.0)+yy;
1576 FPT safeFloatDivision( FPT f1, FPT f2 )
1578 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
1579 ? NumericTraits<FPT>::max()
1580 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) ||
1581 f1 == NumericTraits<FPT>::zero()
1582 ? NumericTraits<FPT>::zero()
1598 template <
class T1,
class T2>
1600 closeAtTolerance(T1 l, T2 r,
typename PromoteTraits<T1, T2>::Promote epsilon)
1602 typedef typename PromoteTraits<T1, T2>::Promote T;
1604 return VIGRA_CSTD::fabs(r) <= epsilon;
1606 return VIGRA_CSTD::fabs(l) <= epsilon;
1607 T diff = VIGRA_CSTD::fabs( l - r );
1608 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
1609 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
1611 return (d1 <= epsilon && d2 <= epsilon);
1614 template <
class T1,
class T2>
1615 inline bool closeAtTolerance(T1 l, T2 r)
1617 typedef typename PromoteTraits<T1, T2>::Promote T;
1618 return closeAtTolerance(l, r, T(2.0) * NumericTraits<T>::epsilon());
double ellipticIntegralF(double x, double k)
The incomplete elliptic integral of the first kind.
Definition: mathutil.hxx:817
REAL sin_pi(REAL x)
sin(pi*x).
Definition: mathutil.hxx:1166
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:1654
NormTraits< T >::NormType norm(T const &t)
The norm of a numerical object.
Definition: mathutil.hxx:679
R arg(const FFTWComplex< R > &a)
pahse
Definition: fftw3.hxx:1009
T1 sign(T1 t1, T2 t2)
The binary sign function.
Definition: mathutil.hxx:587
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)
Compute the eigenvalues of a 3x3 real symmetric matrix.
Definition: mathutil.hxx:716
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
V power(const V &x)
Exponentiation to a positive integer power by squaring.
Definition: mathutil.hxx:389
int round(FixedPoint< IntBits, FracBits > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:683
int signi(T t)
The integer sign function.
Definition: mathutil.hxx:570
Definition: accessor.hxx:43
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition: fixedpoint.hxx:1640
bool odd(int t)
Check if an integer is odd.
REAL cos_pi(REAL x)
cos(pi*x).
Definition: mathutil.hxx:1204
double noncentralChi2(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Non-central chi square distribution.
Definition: mathutil.hxx:1035
double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg)
Cumulative non-central chi square distribution (approximate).
Definition: mathutil.hxx:1073
Int32 log2i(UInt32 x)
Compute the base-2 logarithm of an integer.
Definition: mathutil.hxx:323
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:344
double noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Cumulative non-central chi square distribution.
Definition: mathutil.hxx:1053
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:1871
detail::SelectIntegerType< 32, detail::SignedIntTypes >::type Int32
32-bit signed int
Definition: sized_int.hxx:175
REAL legendre(unsigned int l, REAL x)
Legendre polynomial.
Definition: mathutil.hxx:1151
NormTraits< T >::SquaredNormType squaredNorm(T const &t)
The squared norm of a numerical object.
void symmetric2x2Eigenvalues(T a00, T a01, T a11, T *r0, T *r1)
Compute the eigenvalues of a 2x2 real symmetric matrix.
Definition: mathutil.hxx:696
bool even(int t)
Check if an integer is even.
double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Cumulative chi square distribution.
Definition: mathutil.hxx:1019
double gamma(double x)
The gamma function.
Definition: mathutil.hxx:1549
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
double loggamma(double x)
The natural logarithm of the gamma function.
Definition: mathutil.hxx:1565
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)
Signed integer square root.
Definition: mathutil.hxx:500
UInt32 floorPower2(UInt32 x)
Round down to the nearest power of 2.
Definition: mathutil.hxx:279
double chi2(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Chi square distribution.
Definition: mathutil.hxx:1004
UInt32 ceilPower2(UInt32 x)
Round up to the nearest power of 2.
Definition: mathutil.hxx:257
int roundi(double t)
Round and cast to integer.
Definition: mathutil.hxx:240
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
double ellipticIntegralE(double x, double k)
The incomplete elliptic integral of the second kind.
Definition: mathutil.hxx:841
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616