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
169 template <
class REAL>
170 inline bool isinf(REAL v)
172 return _finite(v) == 0;
175 template <
class REAL>
176 inline bool isnan(REAL v)
178 return _isnan(v) != 0;
186 #define VIGRA_DEFINE_SCALAR_DOT(T) \
187 inline NumericTraits<T>::Promote dot(T l, T r) { return l*r; }
189 VIGRA_DEFINE_SCALAR_DOT(
unsigned char)
190 VIGRA_DEFINE_SCALAR_DOT(
unsigned short)
191 VIGRA_DEFINE_SCALAR_DOT(
unsigned int)
192 VIGRA_DEFINE_SCALAR_DOT(
unsigned long)
193 VIGRA_DEFINE_SCALAR_DOT(
unsigned long long)
194 VIGRA_DEFINE_SCALAR_DOT(
signed char)
195 VIGRA_DEFINE_SCALAR_DOT(
signed short)
196 VIGRA_DEFINE_SCALAR_DOT(
signed int)
197 VIGRA_DEFINE_SCALAR_DOT(
signed long)
198 VIGRA_DEFINE_SCALAR_DOT(
signed long long)
199 VIGRA_DEFINE_SCALAR_DOT(
float)
200 VIGRA_DEFINE_SCALAR_DOT(
double)
201 VIGRA_DEFINE_SCALAR_DOT(
long double)
203 #undef VIGRA_DEFINE_SCALAR_DOT
209 inline float pow(
float v,
double e)
211 return std::pow(v, (
float)e);
214 inline long double pow(
long double v,
double e)
216 return std::pow(v, (
long double)e);
227 #ifdef DOXYGEN // only for documentation
231 inline float round(
float t)
238 inline double round(
double t)
245 inline long double round(
long double t)
315 static Int32 table[64];
319 Int32 IntLog2<T>::table[64] = {
320 -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21,
321 29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6,
322 -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9,
323 11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25,
324 -1, 7, 24, -1, 23, -1, 31, -1};
353 return detail::IntLog2<Int32>::table[x >> 26];
365 typename NumericTraits<T>::Promote
sq(T t)
372 template <
class V,
unsigned>
375 static V call(
const V & x,
const V & y) {
return x * y; }
378 struct cond_mult<V, 0>
380 static V call(
const V &,
const V & y) {
return y; }
383 template <
class V,
unsigned n>
386 static V call(
const V & x)
389 ? cond_mult<V, n & 1>::call(x, power_static<V, n / 2>::call(x * x))
394 struct power_static<V, 0>
396 static V call(
const V & )
409 template <
unsigned n,
class V>
412 return detail::power_static<V, n>::call(x);
421 static UInt32 sqq_table[];
426 UInt32 IntSquareRoot<T>::sqq_table[] = {
427 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57,
428 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83,
429 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102,
430 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
431 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
432 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
433 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
434 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
435 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
436 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
437 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
438 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
439 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
440 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
441 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
442 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
443 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
444 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
455 if (x >= 0x40000000) {
458 xn = sqq_table[x>>24] << 8;
460 xn = sqq_table[x>>22] << 7;
463 xn = sqq_table[x>>20] << 6;
465 xn = sqq_table[x>>18] << 5;
469 xn = sqq_table[x>>16] << 4;
471 xn = sqq_table[x>>14] << 3;
474 xn = sqq_table[x>>12] << 2;
476 xn = sqq_table[x>>10] << 1;
484 xn = (sqq_table[x>>8] >> 0) + 1;
486 xn = (sqq_table[x>>6] >> 1) + 1;
489 xn = (sqq_table[x>>4] >> 2) + 1;
491 xn = (sqq_table[x>>2] >> 3) + 1;
495 return sqq_table[x] >> 4;
499 xn = (xn + 1 + x / xn) / 2;
501 xn = (xn + 1 + x / xn) / 2;
524 throw std::domain_error(
"sqrti(Int32): negative argument.");
525 return (
Int32)detail::IntSquareRoot<UInt32>::exec((
UInt32)v);
537 return detail::IntSquareRoot<UInt32>::exec(v);
540 #ifdef VIGRA_NO_HYPOT
549 inline double hypot(
double a,
double b)
551 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b);
576 return t > NumericTraits<T>::zero()
577 ? NumericTraits<T>::one()
578 : t < NumericTraits<T>::zero()
579 ? -NumericTraits<T>::one()
580 : NumericTraits<T>::zero();
593 return t > NumericTraits<T>::zero()
595 : t < NumericTraits<T>::zero()
607 template <
class T1,
class T2>
610 return t2 >= NumericTraits<T2>::zero()
616 #ifdef DOXYGEN // only for documentation
631 #define VIGRA_DEFINE_ODD_EVEN(T) \
632 inline bool even(T t) { return (t&1) == 0; } \
633 inline bool odd(T t) { return (t&1) == 1; }
635 VIGRA_DEFINE_ODD_EVEN(
char)
636 VIGRA_DEFINE_ODD_EVEN(
short)
637 VIGRA_DEFINE_ODD_EVEN(
int)
638 VIGRA_DEFINE_ODD_EVEN(
long)
639 VIGRA_DEFINE_ODD_EVEN(
long long)
640 VIGRA_DEFINE_ODD_EVEN(
unsigned char)
641 VIGRA_DEFINE_ODD_EVEN(
unsigned short)
642 VIGRA_DEFINE_ODD_EVEN(
unsigned int)
643 VIGRA_DEFINE_ODD_EVEN(
unsigned long)
644 VIGRA_DEFINE_ODD_EVEN(
unsigned long long)
646 #undef VIGRA_DEFINE_ODD_EVEN
648 #define VIGRA_DEFINE_NORM(T) \
649 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
650 inline NormTraits<T>::NormType norm(T t) { return abs(t); }
652 VIGRA_DEFINE_NORM(
bool)
653 VIGRA_DEFINE_NORM(
signed char)
654 VIGRA_DEFINE_NORM(
unsigned char)
655 VIGRA_DEFINE_NORM(
short)
656 VIGRA_DEFINE_NORM(
unsigned short)
657 VIGRA_DEFINE_NORM(
int)
658 VIGRA_DEFINE_NORM(
unsigned int)
659 VIGRA_DEFINE_NORM(
long)
660 VIGRA_DEFINE_NORM(
unsigned long)
661 VIGRA_DEFINE_NORM(
long long)
662 VIGRA_DEFINE_NORM(
unsigned long long)
663 VIGRA_DEFINE_NORM(
float)
664 VIGRA_DEFINE_NORM(
double)
665 VIGRA_DEFINE_NORM(
long double)
667 #undef VIGRA_DEFINE_NORM
670 inline typename NormTraits<std::complex<T> >::SquaredNormType
673 return sq(t.real()) +
sq(t.imag());
676 #ifdef DOXYGEN // only for documentation
686 NormTraits<T>::SquaredNormType
squaredNorm(T
const & t);
699 inline typename NormTraits<T>::NormType
702 typedef typename NormTraits<T>::SquaredNormType SNT;
703 return sqrt(
static_cast<typename SquareRootTraits<SNT>::SquareRootArgument
>(
squaredNorm(t)));
719 double d =
hypot(a00 - a11, 2.0*a01);
720 *r0 =
static_cast<T
>(0.5*(a00 + a11 + d));
721 *r1 =
static_cast<T
>(0.5*(a00 + a11 - d));
738 T * r0, T * r1, T * r2)
740 double inv3 = 1.0 / 3.0, root3 =
std::sqrt(3.0);
742 double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
743 double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
744 double c2 = a00 + a11 + a22;
745 double c2Div3 = c2*inv3;
746 double aDiv3 = (c1 - c2*c2Div3)*inv3;
749 double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
750 double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
757 *r0 =
static_cast<T
>(c2Div3 + 2.0*magnitude*cs);
758 *r1 =
static_cast<T
>(c2Div3 - magnitude*(cs + root3*sn));
759 *r2 =
static_cast<T
>(c2Div3 - magnitude*(cs - root3*sn));
771 T ellipticRD(T x, T y, T z)
773 double f = 1.0, s = 0.0, X, Y, Z, m;
776 m = (x + y + 3.0*z) / 5.0;
780 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
792 double d = a - 6.0*b;
793 double e = d + 2.0*c;
794 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0)
795 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5);
799 T ellipticRF(T x, T y, T z)
804 m = (x + y + z) / 3.0;
808 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
815 double d = X*Y -
sq(Z);
842 return s*detail::ellipticRF(c2, 1.0 -
sq(k*s), 1.0);
867 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0));
870 #if defined(_MSC_VER) && _MSC_VER < 1800
877 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
878 double ans = t*
VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
879 t*(0.09678418+t*(-0.18628806+t*(0.27886807+
880 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
881 t*0.17087277)))))))));
905 inline double erf(
double x)
907 return detail::erfImpl(x);
921 double a = degreesOfFreedom + noncentrality;
922 double b = (a + noncentrality) /
sq(a);
923 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);
928 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans,
unsigned int & j)
938 dans = dans * arg / j;
945 std::pair<double, double>
noncentralChi2CDF(
unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
947 vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,
948 "noncentralChi2P(): parameters must be positive.");
949 if (arg == 0.0 && degreesOfFreedom > 0)
950 return std::make_pair(0.0, 0.0);
953 double b1 = 0.5 * noncentrality,
956 lnrtpi2 = 0.22579135264473,
957 probability, density, lans, dans, pans,
sum, am, hold;
958 unsigned int maxit = 500,
960 if(degreesOfFreedom % 2)
976 if(degreesOfFreedom == 0)
979 degreesOfFreedom = 2;
981 sum = 1.0 / ao - 1.0 - am;
983 probability = 1.0 + am * pans;
988 degreesOfFreedom = degreesOfFreedom - 1;
990 sum = 1.0 / ao - 1.0;
991 while(i < degreesOfFreedom)
992 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
993 degreesOfFreedom = degreesOfFreedom + 1;
998 for(++m; m<maxit; ++m)
1001 detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
1003 density = density + am * dans;
1005 probability = probability + hold;
1006 if((pans * sum < eps2) && (hold < eps2))
1010 vigra_fail(
"noncentralChi2P(): no convergence.");
1011 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
1025 inline double chi2(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1040 inline double chi2CDF(
unsigned int degreesOfFreedom,
double arg,
double accuracy = 1e-7)
1057 double noncentrality,
double arg,
double accuracy = 1e-7)
1075 double noncentrality,
double arg,
double accuracy = 1e-7)
1106 T tmp = NumericTraits<T>::one();
1107 for(T f = l-m+1; f <= l+m; ++f)
1124 template <
class REAL>
1127 vigra_precondition(
abs(x) <= 1.0,
"legendre(): x must be in [-1.0, 1.0].");
1134 return legendre(l,m,x) * s / detail::facLM<REAL>(l,m);
1139 REAL r =
std::sqrt( (1.0-x) * (1.0+x) );
1141 for (
int i=1; i<=m; i++)
1150 REAL result_1 = x * (2.0 * m + 1.0) * result;
1154 for(
unsigned int i = m+2; i <= l; ++i)
1156 other = ( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m);
1171 template <
class REAL>
1186 template <
class REAL>
1194 bool invert =
false;
1208 rem = NumericTraits<REAL>::one();
1224 template <
class REAL>
1232 template <
class REAL>
1235 static REAL
gamma(REAL x);
1248 template <
class REAL>
1249 double GammaImpl<REAL>::g[] = {
1252 -0.6558780715202538,
1253 -0.420026350340952e-1,
1255 -0.421977345555443e-1,
1256 -0.9621971527877e-2,
1258 -0.11651675918591e-2,
1259 -0.2152416741149e-3,
1277 template <
class REAL>
1278 double GammaImpl<REAL>::a[] = {
1279 7.72156649015328655494e-02,
1280 3.22467033424113591611e-01,
1281 6.73523010531292681824e-02,
1282 2.05808084325167332806e-02,
1283 7.38555086081402883957e-03,
1284 2.89051383673415629091e-03,
1285 1.19270763183362067845e-03,
1286 5.10069792153511336608e-04,
1287 2.20862790713908385557e-04,
1288 1.08011567247583939954e-04,
1289 2.52144565451257326939e-05,
1290 4.48640949618915160150e-05
1293 template <
class REAL>
1294 double GammaImpl<REAL>::t[] = {
1295 4.83836122723810047042e-01,
1296 -1.47587722994593911752e-01,
1297 6.46249402391333854778e-02,
1298 -3.27885410759859649565e-02,
1299 1.79706750811820387126e-02,
1300 -1.03142241298341437450e-02,
1301 6.10053870246291332635e-03,
1302 -3.68452016781138256760e-03,
1303 2.25964780900612472250e-03,
1304 -1.40346469989232843813e-03,
1305 8.81081882437654011382e-04,
1306 -5.38595305356740546715e-04,
1307 3.15632070903625950361e-04,
1308 -3.12754168375120860518e-04,
1309 3.35529192635519073543e-04
1312 template <
class REAL>
1313 double GammaImpl<REAL>::u[] = {
1314 -7.72156649015328655494e-02,
1315 6.32827064025093366517e-01,
1316 1.45492250137234768737e+00,
1317 9.77717527963372745603e-01,
1318 2.28963728064692451092e-01,
1319 1.33810918536787660377e-02
1322 template <
class REAL>
1323 double GammaImpl<REAL>::v[] = {
1325 2.45597793713041134822e+00,
1326 2.12848976379893395361e+00,
1327 7.69285150456672783825e-01,
1328 1.04222645593369134254e-01,
1329 3.21709242282423911810e-03
1332 template <
class REAL>
1333 double GammaImpl<REAL>::s[] = {
1334 -7.72156649015328655494e-02,
1335 2.14982415960608852501e-01,
1336 3.25778796408930981787e-01,
1337 1.46350472652464452805e-01,
1338 2.66422703033638609560e-02,
1339 1.84028451407337715652e-03,
1340 3.19475326584100867617e-05
1343 template <
class REAL>
1344 double GammaImpl<REAL>::r[] = {
1346 1.39200533467621045958e+00,
1347 7.21935547567138069525e-01,
1348 1.71933865632803078993e-01,
1349 1.86459191715652901344e-02,
1350 7.77942496381893596434e-04,
1351 7.32668430744625636189e-06
1354 template <
class REAL>
1355 double GammaImpl<REAL>::w[] = {
1356 4.18938533204672725052e-01,
1357 8.33333333333329678849e-02,
1358 -2.77777777728775536470e-03,
1359 7.93650558643019558500e-04,
1360 -5.95187557450339963135e-04,
1361 8.36339918996282139126e-04,
1362 -1.63092934096575273989e-03
1365 template <
class REAL>
1368 int i, k, m, ix = (int)x;
1369 double ga = 0.0, gr = 0.0, r = 0.0, z = 0.0;
1371 vigra_precondition(x <= 171.0,
1372 "gamma(): argument cannot exceed 171.0.");
1379 for (i=2; i<ix; ++i)
1386 vigra_precondition(
false,
1387 "gamma(): gamma function is undefined for 0 and negative integers.");
1397 for (k=1; k<=m; ++k)
1408 for (k=23; k>=0; --k)
1418 ga = -M_PI/(x*ga*
sin_pi(x));
1438 template <
class REAL>
1441 vigra_precondition(x > 0.0,
1442 "loggamma(): argument must be positive.");
1444 vigra_precondition(x <= 1.0e307,
1445 "loggamma(): argument must not exceed 1e307.");
1449 if (x < 4.2351647362715017e-22)
1453 else if ((x == 2.0) || (x == 1.0))
1459 const double tc = 1.46163214496836224576e+00;
1460 const double tf = -1.21486290535849611461e-01;
1461 const double tt = -3.63867699703950536541e-18;
1469 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1470 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1474 else if (x >= 0.23164)
1476 double y = x-(tc-1.0);
1479 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1480 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1481 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1482 double p = z*p1-(tt-w*(p2+y*p3));
1488 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1489 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1490 res += (-0.5*y + p1/p2);
1500 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1501 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1505 else if(x >= 1.23164)
1510 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1511 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1512 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1513 double p = z*p1-(tt-w*(p2+y*p3));
1519 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1520 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1521 res += (-0.5*y + p1/p2);
1529 double p = y*(s[0]+y*(s[1]+y*(s[2]+y*(s[3]+y*(s[4]+y*(s[5]+y*s[6]))))));
1530 double q = 1.0+y*(r[1]+y*(r[2]+y*(r[3]+y*(r[4]+y*(r[5]+y*r[6])))));
1540 else if (x < 2.8823037615171174e+17)
1545 double yy = w[0]+z*(w[1]+y*(w[2]+y*(w[3]+y*(w[4]+y*(w[5]+y*w[6])))));
1546 res = (x-0.5)*(t-1.0)+yy;
1597 FPT safeFloatDivision( FPT f1, FPT f2 )
1599 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
1600 ? NumericTraits<FPT>::max()
1601 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) ||
1602 f1 == NumericTraits<FPT>::zero()
1603 ? NumericTraits<FPT>::zero()
1619 template <
class T1,
class T2>
1623 typedef typename PromoteTraits<T1, T2>::Promote T;
1625 return VIGRA_CSTD::fabs(r) <= epsilon;
1627 return VIGRA_CSTD::fabs(l) <= epsilon;
1628 T diff = VIGRA_CSTD::fabs( l - r );
1629 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
1630 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
1632 return (d1 <= epsilon && d2 <= epsilon);
1635 template <
class T1,
class T2>
1638 typedef typename PromoteTraits<T1, T2>::Promote T;
1653 template <
class T1,
class T2>
1660 template <
class T1,
class T2>
1663 typedef typename PromoteTraits<T1, T2>::Promote T;
1678 template <
class T1,
class T2>
1685 template <
class T1,
class T2>
1688 typedef typename PromoteTraits<T1, T2>::Promote T;
1694 #define VIGRA_MATH_FUNC_HELPER(TYPE) \
1695 inline TYPE clipLower(const TYPE t){ \
1696 return t < static_cast<TYPE>(0.0) ? static_cast<TYPE>(0.0) : t; \
1698 inline TYPE clipLower(const TYPE t,const TYPE valLow){ \
1699 return t < static_cast<TYPE>(valLow) ? static_cast<TYPE>(valLow) : t; \
1701 inline TYPE clipUpper(const TYPE t,const TYPE valHigh){ \
1702 return t > static_cast<TYPE>(valHigh) ? static_cast<TYPE>(valHigh) : t; \
1704 inline TYPE clip(const TYPE t,const TYPE valLow, const TYPE valHigh){ \
1707 else if(t>valHigh) \
1712 inline TYPE sum(const TYPE t){ \
1715 inline NumericTraits<TYPE>::RealPromote mean(const TYPE t){ \
1718 inline TYPE isZero(const TYPE t){ \
1719 return t==static_cast<TYPE>(0); \
1721 inline NumericTraits<TYPE>::RealPromote sizeDividedSquaredNorm(const TYPE t){ \
1722 return squaredNorm(t); \
1724 inline NumericTraits<TYPE>::RealPromote sizeDividedNorm(const TYPE t){ \
1729 VIGRA_MATH_FUNC_HELPER(
unsigned char)
1730 VIGRA_MATH_FUNC_HELPER(
unsigned short)
1731 VIGRA_MATH_FUNC_HELPER(
unsigned int)
1732 VIGRA_MATH_FUNC_HELPER(
unsigned long)
1733 VIGRA_MATH_FUNC_HELPER(
unsigned long long)
1734 VIGRA_MATH_FUNC_HELPER(
signed char)
1735 VIGRA_MATH_FUNC_HELPER(
signed short)
1736 VIGRA_MATH_FUNC_HELPER(
signed int)
1737 VIGRA_MATH_FUNC_HELPER(
signed long)
1738 VIGRA_MATH_FUNC_HELPER(
signed long long)
1739 VIGRA_MATH_FUNC_HELPER(
float)
1740 VIGRA_MATH_FUNC_HELPER(
double)
1741 VIGRA_MATH_FUNC_HELPER(
long double)
1745 #undef VIGRA_MATH_FUNC_HELPER