[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/mathutil.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2011 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_MATHUTIL_HXX 00037 #define VIGRA_MATHUTIL_HXX 00038 00039 #ifdef _MSC_VER 00040 # pragma warning (disable: 4996) // hypot/_hypot confusion 00041 #endif 00042 00043 #include <cmath> 00044 #include <cstdlib> 00045 #include <complex> 00046 #include "config.hxx" 00047 #include "error.hxx" 00048 #include "tuple.hxx" 00049 #include "sized_int.hxx" 00050 #include "numerictraits.hxx" 00051 #include "algorithm.hxx" 00052 00053 /*! \page MathConstants Mathematical Constants 00054 00055 <TT>M_PI, M_SQRT2 etc.</TT> 00056 00057 <b>\#include</b> <vigra/mathutil.hxx> 00058 00059 Since mathematical constants such as <TT>M_PI</TT> and <TT>M_SQRT2</TT> 00060 are not officially standardized, we provide definitions here for those 00061 compilers that don't support them. 00062 00063 \code 00064 #ifndef M_PI 00065 # define M_PI 3.14159265358979323846 00066 #endif 00067 00068 #ifndef M_SQRT2 00069 # define M_2_PI 0.63661977236758134308 00070 #endif 00071 00072 #ifndef M_PI_2 00073 # define M_PI_2 1.57079632679489661923 00074 #endif 00075 00076 #ifndef M_PI_4 00077 # define M_PI_4 0.78539816339744830962 00078 #endif 00079 00080 #ifndef M_SQRT2 00081 # define M_SQRT2 1.41421356237309504880 00082 #endif 00083 00084 #ifndef M_EULER_GAMMA 00085 # define M_EULER_GAMMA 0.5772156649015329 00086 #endif 00087 \endcode 00088 */ 00089 #ifndef M_PI 00090 # define M_PI 3.14159265358979323846 00091 #endif 00092 00093 #ifndef M_2_PI 00094 # define M_2_PI 0.63661977236758134308 00095 #endif 00096 00097 #ifndef M_PI_2 00098 # define M_PI_2 1.57079632679489661923 00099 #endif 00100 00101 #ifndef M_PI_4 00102 # define M_PI_4 0.78539816339744830962 00103 #endif 00104 00105 #ifndef M_SQRT2 00106 # define M_SQRT2 1.41421356237309504880 00107 #endif 00108 00109 #ifndef M_EULER_GAMMA 00110 # define M_EULER_GAMMA 0.5772156649015329 00111 #endif 00112 00113 namespace vigra { 00114 00115 /** \addtogroup MathFunctions Mathematical Functions 00116 00117 Useful mathematical functions and functors. 00118 */ 00119 //@{ 00120 // import functions into namespace vigra which VIGRA is going to overload 00121 00122 using VIGRA_CSTD::pow; 00123 using VIGRA_CSTD::floor; 00124 using VIGRA_CSTD::ceil; 00125 using VIGRA_CSTD::exp; 00126 00127 // import abs(float), abs(double), abs(long double) from <cmath> 00128 // abs(int), abs(long), abs(long long) from <cstdlib> 00129 // abs(std::complex<T>) from <complex> 00130 using std::abs; 00131 00132 // define the missing variants of abs() to avoid 'ambiguous overload' 00133 // errors in template functions 00134 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \ 00135 inline T abs(T t) { return t; } 00136 00137 VIGRA_DEFINE_UNSIGNED_ABS(bool) 00138 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char) 00139 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short) 00140 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int) 00141 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long) 00142 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long long) 00143 00144 #undef VIGRA_DEFINE_UNSIGNED_ABS 00145 00146 #define VIGRA_DEFINE_MISSING_ABS(T) \ 00147 inline T abs(T t) { return t < 0 ? static_cast<T>(-t) : t; } 00148 00149 VIGRA_DEFINE_MISSING_ABS(signed char) 00150 VIGRA_DEFINE_MISSING_ABS(signed short) 00151 00152 #if defined(_MSC_VER) && _MSC_VER < 1600 00153 VIGRA_DEFINE_MISSING_ABS(signed long long) 00154 #endif 00155 00156 #undef VIGRA_DEFINE_MISSING_ABS 00157 00158 // scalar dot is needed for generic functions that should work with 00159 // scalars and vectors alike 00160 00161 #define VIGRA_DEFINE_SCALAR_DOT(T) \ 00162 inline NumericTraits<T>::Promote dot(T l, T r) { return l*r; } 00163 00164 VIGRA_DEFINE_SCALAR_DOT(unsigned char) 00165 VIGRA_DEFINE_SCALAR_DOT(unsigned short) 00166 VIGRA_DEFINE_SCALAR_DOT(unsigned int) 00167 VIGRA_DEFINE_SCALAR_DOT(unsigned long) 00168 VIGRA_DEFINE_SCALAR_DOT(unsigned long long) 00169 VIGRA_DEFINE_SCALAR_DOT(signed char) 00170 VIGRA_DEFINE_SCALAR_DOT(signed short) 00171 VIGRA_DEFINE_SCALAR_DOT(signed int) 00172 VIGRA_DEFINE_SCALAR_DOT(signed long) 00173 VIGRA_DEFINE_SCALAR_DOT(signed long long) 00174 VIGRA_DEFINE_SCALAR_DOT(float) 00175 VIGRA_DEFINE_SCALAR_DOT(double) 00176 VIGRA_DEFINE_SCALAR_DOT(long double) 00177 00178 #undef VIGRA_DEFINE_SCALAR_DOT 00179 00180 using std::pow; 00181 00182 // support 'double' exponents for all floating point versions of pow() 00183 00184 inline float pow(float v, double e) 00185 { 00186 return std::pow(v, (float)e); 00187 } 00188 00189 inline long double pow(long double v, double e) 00190 { 00191 return std::pow(v, (long double)e); 00192 } 00193 00194 /*! The rounding function. 00195 00196 Defined for all floating point types. Rounds towards the nearest integer 00197 such that <tt>abs(round(t)) == round(abs(t))</tt> for all <tt>t</tt>. 00198 00199 <b>\#include</b> <vigra/mathutil.hxx><br> 00200 Namespace: vigra 00201 */ 00202 #ifdef DOXYGEN // only for documentation 00203 REAL round(REAL v); 00204 #endif 00205 00206 inline float round(float t) 00207 { 00208 return t >= 0.0 00209 ? floor(t + 0.5f) 00210 : ceil(t - 0.5f); 00211 } 00212 00213 inline double round(double t) 00214 { 00215 return t >= 0.0 00216 ? floor(t + 0.5) 00217 : ceil(t - 0.5); 00218 } 00219 00220 inline long double round(long double t) 00221 { 00222 return t >= 0.0 00223 ? floor(t + 0.5) 00224 : ceil(t - 0.5); 00225 } 00226 00227 00228 /*! Round and cast to integer. 00229 00230 Rounds to the nearest integer like round(), but casts the result to 00231 <tt>int</tt> (this will be faster and is usually needed anyway). 00232 00233 <b>\#include</b> <vigra/mathutil.hxx><br> 00234 Namespace: vigra 00235 */ 00236 inline int roundi(double t) 00237 { 00238 return t >= 0.0 00239 ? int(t + 0.5) 00240 : int(t - 0.5); 00241 } 00242 00243 /*! Round up to the nearest power of 2. 00244 00245 Efficient algorithm for finding the smallest power of 2 which is not smaller than \a x 00246 (function clp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003, 00247 see http://www.hackersdelight.org/). 00248 If \a x > 2^31, the function will return 0 because integer arithmetic is defined modulo 2^32. 00249 00250 <b>\#include</b> <vigra/mathutil.hxx><br> 00251 Namespace: vigra 00252 */ 00253 inline UInt32 ceilPower2(UInt32 x) 00254 { 00255 if(x == 0) return 0; 00256 00257 x = x - 1; 00258 x = x | (x >> 1); 00259 x = x | (x >> 2); 00260 x = x | (x >> 4); 00261 x = x | (x >> 8); 00262 x = x | (x >>16); 00263 return x + 1; 00264 } 00265 00266 /*! Round down to the nearest power of 2. 00267 00268 Efficient algorithm for finding the largest power of 2 which is not greater than \a x 00269 (function flp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003, 00270 see http://www.hackersdelight.org/). 00271 00272 <b>\#include</b> <vigra/mathutil.hxx><br> 00273 Namespace: vigra 00274 */ 00275 inline UInt32 floorPower2(UInt32 x) 00276 { 00277 x = x | (x >> 1); 00278 x = x | (x >> 2); 00279 x = x | (x >> 4); 00280 x = x | (x >> 8); 00281 x = x | (x >>16); 00282 return x - (x >> 1); 00283 } 00284 00285 namespace detail { 00286 00287 template <class T> 00288 class IntLog2 00289 { 00290 public: 00291 static Int32 table[64]; 00292 }; 00293 00294 template <class T> 00295 Int32 IntLog2<T>::table[64] = { 00296 -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21, 00297 29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6, 00298 -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9, 00299 11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25, 00300 -1, 7, 24, -1, 23, -1, 31, -1}; 00301 00302 } // namespace detail 00303 00304 /*! Compute the base-2 logarithm of an integer. 00305 00306 Returns the position of the left-most 1-bit in the given number \a x, or 00307 -1 if \a x == 0. That is, 00308 00309 \code 00310 assert(k >= 0 && k < 32 && log2i(1 << k) == k); 00311 \endcode 00312 00313 The function uses Robert Harley's algorithm to determine the number of leading zeros 00314 in \a x (algorithm nlz10() at http://www.hackersdelight.org/). But note that the functions 00315 \ref floorPower2() or \ref ceilPower2() are more efficient and should be preferred when possible. 00316 00317 <b>\#include</b> <vigra/mathutil.hxx><br> 00318 Namespace: vigra 00319 */ 00320 inline Int32 log2i(UInt32 x) 00321 { 00322 // Propagate leftmost 1-bit to the right. 00323 x = x | (x >> 1); 00324 x = x | (x >> 2); 00325 x = x | (x >> 4); 00326 x = x | (x >> 8); 00327 x = x | (x >>16); 00328 x = x*0x06EB14F9; // Multiplier is 7*255**3. 00329 return detail::IntLog2<Int32>::table[x >> 26]; 00330 } 00331 00332 /*! The square function. 00333 00334 <tt>sq(x) = x*x</tt> is needed so often that it makes sense to define it as a function. 00335 00336 <b>\#include</b> <vigra/mathutil.hxx><br> 00337 Namespace: vigra 00338 */ 00339 template <class T> 00340 inline 00341 typename NumericTraits<T>::Promote sq(T t) 00342 { 00343 return t*t; 00344 } 00345 00346 namespace detail { 00347 00348 template <class V, unsigned> 00349 struct cond_mult 00350 { 00351 static V call(const V & x, const V & y) { return x * y; } 00352 }; 00353 template <class V> 00354 struct cond_mult<V, 0> 00355 { 00356 static V call(const V &, const V & y) { return y; } 00357 }; 00358 00359 template <class V, unsigned n> 00360 struct power_static 00361 { 00362 static V call(const V & x) 00363 { 00364 return n / 2 00365 ? cond_mult<V, n & 1>::call(x, power_static<V, n / 2>::call(x * x)) 00366 : n & 1 ? x : V(); 00367 } 00368 }; 00369 template <class V> 00370 struct power_static<V, 0> 00371 { 00372 static V call(const V & x) 00373 { 00374 return V(1); 00375 } 00376 }; 00377 00378 } // namespace detail 00379 00380 /*! Exponentiation to a positive integer power by squaring. 00381 00382 <b>\#include</b> <vigra/mathutil.hxx><br> 00383 Namespace: vigra 00384 */ 00385 template <unsigned n, class V> 00386 inline V power(const V & x) 00387 { 00388 return detail::power_static<V, n>::call(x); 00389 } 00390 //doxygen_overloaded_function(template <unsigned n, class V> power(const V & x)) 00391 00392 namespace detail { 00393 00394 template <class T> 00395 class IntSquareRoot 00396 { 00397 public: 00398 static UInt32 sqq_table[]; 00399 static UInt32 exec(UInt32 v); 00400 }; 00401 00402 template <class T> 00403 UInt32 IntSquareRoot<T>::sqq_table[] = { 00404 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 00405 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 00406 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 00407 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118, 00408 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132, 00409 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145, 00410 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157, 00411 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, 00412 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 00413 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 00414 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 00415 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206, 00416 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215, 00417 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223, 00418 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 00419 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, 00420 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 00421 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 00422 253, 254, 254, 255 00423 }; 00424 00425 template <class T> 00426 UInt32 IntSquareRoot<T>::exec(UInt32 x) 00427 { 00428 UInt32 xn; 00429 if (x >= 0x10000) 00430 if (x >= 0x1000000) 00431 if (x >= 0x10000000) 00432 if (x >= 0x40000000) { 00433 if (x >= (UInt32)65535*(UInt32)65535) 00434 return 65535; 00435 xn = sqq_table[x>>24] << 8; 00436 } else 00437 xn = sqq_table[x>>22] << 7; 00438 else 00439 if (x >= 0x4000000) 00440 xn = sqq_table[x>>20] << 6; 00441 else 00442 xn = sqq_table[x>>18] << 5; 00443 else { 00444 if (x >= 0x100000) 00445 if (x >= 0x400000) 00446 xn = sqq_table[x>>16] << 4; 00447 else 00448 xn = sqq_table[x>>14] << 3; 00449 else 00450 if (x >= 0x40000) 00451 xn = sqq_table[x>>12] << 2; 00452 else 00453 xn = sqq_table[x>>10] << 1; 00454 00455 goto nr1; 00456 } 00457 else 00458 if (x >= 0x100) { 00459 if (x >= 0x1000) 00460 if (x >= 0x4000) 00461 xn = (sqq_table[x>>8] >> 0) + 1; 00462 else 00463 xn = (sqq_table[x>>6] >> 1) + 1; 00464 else 00465 if (x >= 0x400) 00466 xn = (sqq_table[x>>4] >> 2) + 1; 00467 else 00468 xn = (sqq_table[x>>2] >> 3) + 1; 00469 00470 goto adj; 00471 } else 00472 return sqq_table[x] >> 4; 00473 00474 /* Run two iterations of the standard convergence formula */ 00475 00476 xn = (xn + 1 + x / xn) / 2; 00477 nr1: 00478 xn = (xn + 1 + x / xn) / 2; 00479 adj: 00480 00481 if (xn * xn > x) /* Correct rounding if necessary */ 00482 xn--; 00483 00484 return xn; 00485 } 00486 00487 } // namespace detail 00488 00489 using VIGRA_CSTD::sqrt; 00490 00491 /*! Signed integer square root. 00492 00493 Useful for fast fixed-point computations. 00494 00495 <b>\#include</b> <vigra/mathutil.hxx><br> 00496 Namespace: vigra 00497 */ 00498 inline Int32 sqrti(Int32 v) 00499 { 00500 if(v < 0) 00501 throw std::domain_error("sqrti(Int32): negative argument."); 00502 return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v); 00503 } 00504 00505 /*! Unsigned integer square root. 00506 00507 Useful for fast fixed-point computations. 00508 00509 <b>\#include</b> <vigra/mathutil.hxx><br> 00510 Namespace: vigra 00511 */ 00512 inline UInt32 sqrti(UInt32 v) 00513 { 00514 return detail::IntSquareRoot<UInt32>::exec(v); 00515 } 00516 00517 #ifdef VIGRA_NO_HYPOT 00518 /*! Compute the Euclidean distance (length of the hypotenuse of a right-angled triangle). 00519 00520 The hypot() function returns the sqrt(a*a + b*b). 00521 It is implemented in a way that minimizes round-off error. 00522 00523 <b>\#include</b> <vigra/mathutil.hxx><br> 00524 Namespace: vigra 00525 */ 00526 inline double hypot(double a, double b) 00527 { 00528 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b); 00529 if (absa > absb) 00530 return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 00531 else 00532 return absb == 0.0 00533 ? 0.0 00534 : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 00535 } 00536 00537 #else 00538 00539 using ::hypot; 00540 00541 #endif 00542 00543 /*! The sign function. 00544 00545 Returns 1, 0, or -1 depending on the sign of \a t, but with the same type as \a t. 00546 00547 <b>\#include</b> <vigra/mathutil.hxx><br> 00548 Namespace: vigra 00549 */ 00550 template <class T> 00551 inline T sign(T t) 00552 { 00553 return t > NumericTraits<T>::zero() 00554 ? NumericTraits<T>::one() 00555 : t < NumericTraits<T>::zero() 00556 ? -NumericTraits<T>::one() 00557 : NumericTraits<T>::zero(); 00558 } 00559 00560 /*! The integer sign function. 00561 00562 Returns 1, 0, or -1 depending on the sign of \a t, converted to int. 00563 00564 <b>\#include</b> <vigra/mathutil.hxx><br> 00565 Namespace: vigra 00566 */ 00567 template <class T> 00568 inline int signi(T t) 00569 { 00570 return t > NumericTraits<T>::zero() 00571 ? 1 00572 : t < NumericTraits<T>::zero() 00573 ? -1 00574 : 0; 00575 } 00576 00577 /*! The binary sign function. 00578 00579 Transfers the sign of \a t2 to \a t1. 00580 00581 <b>\#include</b> <vigra/mathutil.hxx><br> 00582 Namespace: vigra 00583 */ 00584 template <class T1, class T2> 00585 inline T1 sign(T1 t1, T2 t2) 00586 { 00587 return t2 >= NumericTraits<T2>::zero() 00588 ? abs(t1) 00589 : -abs(t1); 00590 } 00591 00592 00593 #ifdef DOXYGEN // only for documentation 00594 /*! Check if an integer is even. 00595 00596 Defined for all integral types. 00597 */ 00598 bool even(int t); 00599 00600 /*! Check if an integer is odd. 00601 00602 Defined for all integral types. 00603 */ 00604 bool odd(int t); 00605 00606 #endif 00607 00608 #define VIGRA_DEFINE_ODD_EVEN(T) \ 00609 inline bool even(T t) { return (t&1) == 0; } \ 00610 inline bool odd(T t) { return (t&1) == 1; } 00611 00612 VIGRA_DEFINE_ODD_EVEN(char) 00613 VIGRA_DEFINE_ODD_EVEN(short) 00614 VIGRA_DEFINE_ODD_EVEN(int) 00615 VIGRA_DEFINE_ODD_EVEN(long) 00616 VIGRA_DEFINE_ODD_EVEN(long long) 00617 VIGRA_DEFINE_ODD_EVEN(unsigned char) 00618 VIGRA_DEFINE_ODD_EVEN(unsigned short) 00619 VIGRA_DEFINE_ODD_EVEN(unsigned int) 00620 VIGRA_DEFINE_ODD_EVEN(unsigned long) 00621 VIGRA_DEFINE_ODD_EVEN(unsigned long long) 00622 00623 #undef VIGRA_DEFINE_ODD_EVEN 00624 00625 #define VIGRA_DEFINE_NORM(T) \ 00626 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \ 00627 inline NormTraits<T>::NormType norm(T t) { return abs(t); } 00628 00629 VIGRA_DEFINE_NORM(bool) 00630 VIGRA_DEFINE_NORM(signed char) 00631 VIGRA_DEFINE_NORM(unsigned char) 00632 VIGRA_DEFINE_NORM(short) 00633 VIGRA_DEFINE_NORM(unsigned short) 00634 VIGRA_DEFINE_NORM(int) 00635 VIGRA_DEFINE_NORM(unsigned int) 00636 VIGRA_DEFINE_NORM(long) 00637 VIGRA_DEFINE_NORM(unsigned long) 00638 VIGRA_DEFINE_NORM(long long) 00639 VIGRA_DEFINE_NORM(unsigned long long) 00640 VIGRA_DEFINE_NORM(float) 00641 VIGRA_DEFINE_NORM(double) 00642 VIGRA_DEFINE_NORM(long double) 00643 00644 #undef VIGRA_DEFINE_NORM 00645 00646 template <class T> 00647 inline typename NormTraits<std::complex<T> >::SquaredNormType 00648 squaredNorm(std::complex<T> const & t) 00649 { 00650 return sq(t.real()) + sq(t.imag()); 00651 } 00652 00653 #ifdef DOXYGEN // only for documentation 00654 /*! The squared norm of a numerical object. 00655 00656 For scalar types: equals <tt>vigra::sq(t)</tt><br>. 00657 For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>. 00658 For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>. 00659 For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements). 00660 */ 00661 NormTraits<T>::SquaredNormType squaredNorm(T const & t); 00662 00663 #endif 00664 00665 /*! The norm of a numerical object. 00666 00667 For scalar types: implemented as <tt>abs(t)</tt><br> 00668 otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>. 00669 00670 <b>\#include</b> <vigra/mathutil.hxx><br> 00671 Namespace: vigra 00672 */ 00673 template <class T> 00674 inline typename NormTraits<T>::NormType 00675 norm(T const & t) 00676 { 00677 typedef typename NormTraits<T>::SquaredNormType SNT; 00678 return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t))); 00679 } 00680 00681 /*! Compute the eigenvalues of a 2x2 real symmetric matrix. 00682 00683 This uses the analytical eigenvalue formula 00684 \f[ 00685 \lambda_{1,2} = \frac{1}{2}\left(a_{00} + a_{11} \pm \sqrt{(a_{00} - a_{11})^2 + 4 a_{01}^2}\right) 00686 \f] 00687 00688 <b>\#include</b> <vigra/mathutil.hxx><br> 00689 Namespace: vigra 00690 */ 00691 template <class T> 00692 void symmetric2x2Eigenvalues(T a00, T a01, T a11, T * r0, T * r1) 00693 { 00694 double d = hypot(a00 - a11, 2.0*a01); 00695 *r0 = static_cast<T>(0.5*(a00 + a11 + d)); 00696 *r1 = static_cast<T>(0.5*(a00 + a11 - d)); 00697 if(*r0 < *r1) 00698 std::swap(*r0, *r1); 00699 } 00700 00701 /*! Compute the eigenvalues of a 3x3 real symmetric matrix. 00702 00703 This uses a numerically stable version of the analytical eigenvalue formula according to 00704 <p> 00705 David Eberly: <a href="http://www.geometrictools.com/Documentation/EigenSymmetric3x3.pdf"> 00706 <em>"Eigensystems for 3 × 3 Symmetric Matrices (Revisited)"</em></a>, Geometric Tools Documentation, 2006 00707 00708 00709 <b>\#include</b> <vigra/mathutil.hxx><br> 00710 Namespace: vigra 00711 */ 00712 template <class T> 00713 void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22, 00714 T * r0, T * r1, T * r2) 00715 { 00716 static double inv3 = 1.0 / 3.0, root3 = std::sqrt(3.0); 00717 00718 double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01; 00719 double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12; 00720 double c2 = a00 + a11 + a22; 00721 double c2Div3 = c2*inv3; 00722 double aDiv3 = (c1 - c2*c2Div3)*inv3; 00723 if (aDiv3 > 0.0) 00724 aDiv3 = 0.0; 00725 double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1)); 00726 double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3; 00727 if (q > 0.0) 00728 q = 0.0; 00729 double magnitude = std::sqrt(-aDiv3); 00730 double angle = std::atan2(std::sqrt(-q),mbDiv2)*inv3; 00731 double cs = std::cos(angle); 00732 double sn = std::sin(angle); 00733 *r0 = static_cast<T>(c2Div3 + 2.0*magnitude*cs); 00734 *r1 = static_cast<T>(c2Div3 - magnitude*(cs + root3*sn)); 00735 *r2 = static_cast<T>(c2Div3 - magnitude*(cs - root3*sn)); 00736 if(*r0 < *r1) 00737 std::swap(*r0, *r1); 00738 if(*r0 < *r2) 00739 std::swap(*r0, *r2); 00740 if(*r1 < *r2) 00741 std::swap(*r1, *r2); 00742 } 00743 00744 namespace detail { 00745 00746 template <class T> 00747 T ellipticRD(T x, T y, T z) 00748 { 00749 double f = 1.0, s = 0.0, X, Y, Z, m; 00750 for(;;) 00751 { 00752 m = (x + y + 3.0*z) / 5.0; 00753 X = 1.0 - x/m; 00754 Y = 1.0 - y/m; 00755 Z = 1.0 - z/m; 00756 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01) 00757 break; 00758 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z); 00759 s += f / (VIGRA_CSTD::sqrt(z)*(z + l)); 00760 f /= 4.0; 00761 x = (x + l)/4.0; 00762 y = (y + l)/4.0; 00763 z = (z + l)/4.0; 00764 } 00765 double a = X*Y; 00766 double b = sq(Z); 00767 double c = a - b; 00768 double d = a - 6.0*b; 00769 double e = d + 2.0*c; 00770 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0) 00771 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5); 00772 } 00773 00774 template <class T> 00775 T ellipticRF(T x, T y, T z) 00776 { 00777 double X, Y, Z, m; 00778 for(;;) 00779 { 00780 m = (x + y + z) / 3.0; 00781 X = 1.0 - x/m; 00782 Y = 1.0 - y/m; 00783 Z = 1.0 - z/m; 00784 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01) 00785 break; 00786 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z); 00787 x = (x + l)/4.0; 00788 y = (y + l)/4.0; 00789 z = (z + l)/4.0; 00790 } 00791 double d = X*Y - sq(Z); 00792 double p = X*Y*Z; 00793 return (1.0 - d/10.0 + p/14.0 + sq(d)/24.0 - d*p*3.0/44.0) / VIGRA_CSTD::sqrt(m); 00794 } 00795 00796 } // namespace detail 00797 00798 /*! The incomplete elliptic integral of the first kind. 00799 00800 Computes 00801 00802 \f[ 00803 \mbox{F}(x, k) = \int_0^x \frac{1}{\sqrt{1 - k^2 \sin(t)^2}} dt 00804 \f] 00805 00806 according to the algorithm given in Press et al. "Numerical Recipes". 00807 00808 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral 00809 functions must be k^2 rather than k. Check the documentation when results disagree! 00810 00811 <b>\#include</b> <vigra/mathutil.hxx><br> 00812 Namespace: vigra 00813 */ 00814 inline double ellipticIntegralF(double x, double k) 00815 { 00816 double c2 = sq(VIGRA_CSTD::cos(x)); 00817 double s = VIGRA_CSTD::sin(x); 00818 return s*detail::ellipticRF(c2, 1.0 - sq(k*s), 1.0); 00819 } 00820 00821 /*! The incomplete elliptic integral of the second kind. 00822 00823 Computes 00824 00825 \f[ 00826 \mbox{E}(x, k) = \int_0^x \sqrt{1 - k^2 \sin(t)^2} dt 00827 \f] 00828 00829 according to the algorithm given in Press et al. "Numerical Recipes". The 00830 complete elliptic integral of the second kind is simply <tt>ellipticIntegralE(M_PI/2, k)</TT>. 00831 00832 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral 00833 functions must be k^2 rather than k. Check the documentation when results disagree! 00834 00835 <b>\#include</b> <vigra/mathutil.hxx><br> 00836 Namespace: vigra 00837 */ 00838 inline double ellipticIntegralE(double x, double k) 00839 { 00840 double c2 = sq(VIGRA_CSTD::cos(x)); 00841 double s = VIGRA_CSTD::sin(x); 00842 k = sq(k*s); 00843 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0)); 00844 } 00845 00846 #ifdef _MSC_VER 00847 00848 namespace detail { 00849 00850 template <class T> 00851 double erfImpl(T x) 00852 { 00853 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x)); 00854 double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+ 00855 t*(0.09678418+t*(-0.18628806+t*(0.27886807+ 00856 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+ 00857 t*0.17087277))))))))); 00858 if (x >= 0.0) 00859 return 1.0 - ans; 00860 else 00861 return ans - 1.0; 00862 } 00863 00864 } // namespace detail 00865 00866 /*! The error function. 00867 00868 If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the 00869 new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error 00870 function 00871 00872 \f[ 00873 \mbox{erf}(x) = \int_0^x e^{-t^2} dt 00874 \f] 00875 00876 according to the formula given in Press et al. "Numerical Recipes". 00877 00878 <b>\#include</b> <vigra/mathutil.hxx><br> 00879 Namespace: vigra 00880 */ 00881 inline double erf(double x) 00882 { 00883 return detail::erfImpl(x); 00884 } 00885 00886 #else 00887 00888 using ::erf; 00889 00890 #endif 00891 00892 namespace detail { 00893 00894 template <class T> 00895 double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, T noncentrality, T arg) 00896 { 00897 double a = degreesOfFreedom + noncentrality; 00898 double b = (a + noncentrality) / sq(a); 00899 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); 00900 return 0.5*(1.0 + erf(t/VIGRA_CSTD::sqrt(2.0))); 00901 } 00902 00903 template <class T> 00904 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j) 00905 { 00906 double tol = -50.0; 00907 if(lans < tol) 00908 { 00909 lans = lans + VIGRA_CSTD::log(arg / j); 00910 dans = VIGRA_CSTD::exp(lans); 00911 } 00912 else 00913 { 00914 dans = dans * arg / j; 00915 } 00916 pans = pans - dans; 00917 j += 2; 00918 } 00919 00920 template <class T> 00921 std::pair<double, double> noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps) 00922 { 00923 vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0, 00924 "noncentralChi2P(): parameters must be positive."); 00925 if (arg == 0.0 && degreesOfFreedom > 0) 00926 return std::make_pair(0.0, 0.0); 00927 00928 // Determine initial values 00929 double b1 = 0.5 * noncentrality, 00930 ao = VIGRA_CSTD::exp(-b1), 00931 eps2 = eps / ao, 00932 lnrtpi2 = 0.22579135264473, 00933 probability, density, lans, dans, pans, sum, am, hold; 00934 unsigned int maxit = 500, 00935 i, m; 00936 if(degreesOfFreedom % 2) 00937 { 00938 i = 1; 00939 lans = -0.5 * (arg + VIGRA_CSTD::log(arg)) - lnrtpi2; 00940 dans = VIGRA_CSTD::exp(lans); 00941 pans = erf(VIGRA_CSTD::sqrt(arg/2.0)); 00942 } 00943 else 00944 { 00945 i = 2; 00946 lans = -0.5 * arg; 00947 dans = VIGRA_CSTD::exp(lans); 00948 pans = 1.0 - dans; 00949 } 00950 00951 // Evaluate first term 00952 if(degreesOfFreedom == 0) 00953 { 00954 m = 1; 00955 degreesOfFreedom = 2; 00956 am = b1; 00957 sum = 1.0 / ao - 1.0 - am; 00958 density = am * dans; 00959 probability = 1.0 + am * pans; 00960 } 00961 else 00962 { 00963 m = 0; 00964 degreesOfFreedom = degreesOfFreedom - 1; 00965 am = 1.0; 00966 sum = 1.0 / ao - 1.0; 00967 while(i < degreesOfFreedom) 00968 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i); 00969 degreesOfFreedom = degreesOfFreedom + 1; 00970 density = dans; 00971 probability = pans; 00972 } 00973 // Evaluate successive terms of the expansion 00974 for(++m; m<maxit; ++m) 00975 { 00976 am = b1 * am / m; 00977 detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom); 00978 sum = sum - am; 00979 density = density + am * dans; 00980 hold = am * pans; 00981 probability = probability + hold; 00982 if((pans * sum < eps2) && (hold < eps2)) 00983 break; // converged 00984 } 00985 if(m == maxit) 00986 vigra_fail("noncentralChi2P(): no convergence."); 00987 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability))); 00988 } 00989 00990 } // namespace detail 00991 00992 /*! Chi square distribution. 00993 00994 Computes the density of a chi square distribution with \a degreesOfFreedom 00995 and tolerance \a accuracy at the given argument \a arg 00996 by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>. 00997 00998 <b>\#include</b> <vigra/mathutil.hxx><br> 00999 Namespace: vigra 01000 */ 01001 inline double chi2(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7) 01002 { 01003 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).first; 01004 } 01005 01006 /*! Cumulative chi square distribution. 01007 01008 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom 01009 and tolerance \a accuracy at the given argument \a arg, i.e. the probability that 01010 a random number drawn from the distribution is below \a arg 01011 by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>. 01012 01013 <b>\#include</b> <vigra/mathutil.hxx><br> 01014 Namespace: vigra 01015 */ 01016 inline double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7) 01017 { 01018 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).second; 01019 } 01020 01021 /*! Non-central chi square distribution. 01022 01023 Computes the density of a non-central chi square distribution with \a degreesOfFreedom, 01024 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 01025 \a arg. It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 01026 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of 01027 degrees of freedom. 01028 01029 <b>\#include</b> <vigra/mathutil.hxx><br> 01030 Namespace: vigra 01031 */ 01032 inline double noncentralChi2(unsigned int degreesOfFreedom, 01033 double noncentrality, double arg, double accuracy = 1e-7) 01034 { 01035 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).first; 01036 } 01037 01038 /*! Cumulative non-central chi square distribution. 01039 01040 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom, 01041 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 01042 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg 01043 It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 01044 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of 01045 degrees of freedom (see noncentralChi2CDFApprox() for a constant-time algorithm). 01046 01047 <b>\#include</b> <vigra/mathutil.hxx><br> 01048 Namespace: vigra 01049 */ 01050 inline double noncentralChi2CDF(unsigned int degreesOfFreedom, 01051 double noncentrality, double arg, double accuracy = 1e-7) 01052 { 01053 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).second; 01054 } 01055 01056 /*! Cumulative non-central chi square distribution (approximate). 01057 01058 Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom, 01059 and noncentrality parameter \a noncentrality at the given argument 01060 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg 01061 It uses the approximate transform into a normal distribution due to Wilson and Hilferty 01062 (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32). 01063 The algorithm's running time is independent of the inputs, i.e. is should be used 01064 when noncentralChi2CDF() is too slow, and approximate values are sufficient. The accuracy is only 01065 about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5. 01066 01067 <b>\#include</b> <vigra/mathutil.hxx><br> 01068 Namespace: vigra 01069 */ 01070 inline double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg) 01071 { 01072 return detail::noncentralChi2CDFApprox(degreesOfFreedom, noncentrality, arg); 01073 } 01074 01075 namespace detail { 01076 01077 // computes (l+m)! / (l-m)! 01078 // l and m must be positive 01079 template <class T> 01080 T facLM(T l, T m) 01081 { 01082 T tmp = NumericTraits<T>::one(); 01083 for(T f = l-m+1; f <= l+m; ++f) 01084 tmp *= f; 01085 return tmp; 01086 } 01087 01088 } // namespace detail 01089 01090 /*! Associated Legendre polynomial. 01091 01092 Computes the value of the associated Legendre polynomial of order <tt>l, m</tt> 01093 for argument <tt>x</tt>. <tt>x</tt> must be in the range <tt>[-1.0, 1.0]</tt>, 01094 otherwise an exception is thrown. The standard Legendre polynomials are the 01095 special case <tt>m == 0</tt>. 01096 01097 <b>\#include</b> <vigra/mathutil.hxx><br> 01098 Namespace: vigra 01099 */ 01100 template <class REAL> 01101 REAL legendre(unsigned int l, int m, REAL x) 01102 { 01103 vigra_precondition(abs(x) <= 1.0, "legendre(): x must be in [-1.0, 1.0]."); 01104 if (m < 0) 01105 { 01106 m = -m; 01107 REAL s = odd(m) 01108 ? -1.0 01109 : 1.0; 01110 return legendre(l,m,x) * s / detail::facLM<REAL>(l,m); 01111 } 01112 REAL result = 1.0; 01113 if (m > 0) 01114 { 01115 REAL r = std::sqrt( (1.0-x) * (1.0+x) ); 01116 REAL f = 1.0; 01117 for (int i=1; i<=m; i++) 01118 { 01119 result *= (-f) * r; 01120 f += 2.0; 01121 } 01122 } 01123 if((int)l == m) 01124 return result; 01125 01126 REAL result_1 = x * (2.0 * m + 1.0) * result; 01127 if((int)l == m+1) 01128 return result_1; 01129 REAL other = 0.0; 01130 for(unsigned int i = m+2; i <= l; ++i) 01131 { 01132 other = ( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m); 01133 result = result_1; 01134 result_1 = other; 01135 } 01136 return other; 01137 } 01138 01139 /*! Legendre polynomial. 01140 01141 Computes the value of the Legendre polynomial of order <tt>l</tt> for argument <tt>x</tt>. 01142 <tt>x</tt> must be in the range <tt>[-1.0, 1.0]</tt>, otherwise an exception is thrown. 01143 01144 <b>\#include</b> <vigra/mathutil.hxx><br> 01145 Namespace: vigra 01146 */ 01147 template <class REAL> 01148 REAL legendre(unsigned int l, REAL x) 01149 { 01150 return legendre(l, 0, x); 01151 } 01152 01153 /*! sin(pi*x). 01154 01155 Essentially calls <tt>std::sin(M_PI*x)</tt> but uses a more accurate implementation 01156 to make sure that <tt>sin_pi(1.0) == 0.0</tt> (which does not hold for 01157 <tt>std::sin(M_PI)</tt> due to round-off error), and <tt>sin_pi(0.5) == 1.0</tt>. 01158 01159 <b>\#include</b> <vigra/mathutil.hxx><br> 01160 Namespace: vigra 01161 */ 01162 template <class REAL> 01163 REAL sin_pi(REAL x) 01164 { 01165 if(x < 0.0) 01166 return -sin_pi(-x); 01167 if(x < 0.5) 01168 return std::sin(M_PI * x); 01169 01170 bool invert = false; 01171 if(x < 1.0) 01172 { 01173 invert = true; 01174 x = -x; 01175 } 01176 01177 REAL rem = std::floor(x); 01178 if(odd((int)rem)) 01179 invert = !invert; 01180 rem = x - rem; 01181 if(rem > 0.5) 01182 rem = 1.0 - rem; 01183 if(rem == 0.5) 01184 rem = NumericTraits<REAL>::one(); 01185 else 01186 rem = std::sin(M_PI * rem); 01187 return invert 01188 ? -rem 01189 : rem; 01190 } 01191 01192 /*! cos(pi*x). 01193 01194 Essentially calls <tt>std::cos(M_PI*x)</tt> but uses a more accurate implementation 01195 to make sure that <tt>cos_pi(1.0) == -1.0</tt> and <tt>cos_pi(0.5) == 0.0</tt>. 01196 01197 <b>\#include</b> <vigra/mathutil.hxx><br> 01198 Namespace: vigra 01199 */ 01200 template <class REAL> 01201 REAL cos_pi(REAL x) 01202 { 01203 return sin_pi(x+0.5); 01204 } 01205 01206 namespace detail { 01207 01208 template <class REAL> 01209 REAL gammaImpl(REAL x) 01210 { 01211 int i, k, m, ix = (int)x; 01212 double ga = 0.0, gr = 0.0, r = 0.0, z = 0.0; 01213 01214 static double g[] = { 01215 1.0, 01216 0.5772156649015329, 01217 -0.6558780715202538, 01218 -0.420026350340952e-1, 01219 0.1665386113822915, 01220 -0.421977345555443e-1, 01221 -0.9621971527877e-2, 01222 0.7218943246663e-2, 01223 -0.11651675918591e-2, 01224 -0.2152416741149e-3, 01225 0.1280502823882e-3, 01226 -0.201348547807e-4, 01227 -0.12504934821e-5, 01228 0.1133027232e-5, 01229 -0.2056338417e-6, 01230 0.6116095e-8, 01231 0.50020075e-8, 01232 -0.11812746e-8, 01233 0.1043427e-9, 01234 0.77823e-11, 01235 -0.36968e-11, 01236 0.51e-12, 01237 -0.206e-13, 01238 -0.54e-14, 01239 0.14e-14}; 01240 01241 vigra_precondition(x <= 171.0, 01242 "gamma(): argument cannot exceed 171.0."); 01243 01244 if (x == ix) 01245 { 01246 if (ix > 0) 01247 { 01248 ga = 1.0; // use factorial 01249 for (i=2; i<ix; ++i) 01250 { 01251 ga *= i; 01252 } 01253 } 01254 else 01255 { 01256 vigra_precondition(false, 01257 "gamma(): gamma function is undefined for 0 and negative integers."); 01258 } 01259 } 01260 else 01261 { 01262 if (abs(x) > 1.0) 01263 { 01264 z = abs(x); 01265 m = (int)z; 01266 r = 1.0; 01267 for (k=1; k<=m; ++k) 01268 { 01269 r *= (z-k); 01270 } 01271 z -= m; 01272 } 01273 else 01274 { 01275 z = x; 01276 } 01277 gr = g[24]; 01278 for (k=23; k>=0; --k) 01279 { 01280 gr = gr*z+g[k]; 01281 } 01282 ga = 1.0/(gr*z); 01283 if (abs(x) > 1.0) 01284 { 01285 ga *= r; 01286 if (x < 0.0) 01287 { 01288 ga = -M_PI/(x*ga*sin_pi(x)); 01289 } 01290 } 01291 } 01292 return ga; 01293 } 01294 01295 /* 01296 * the following code is derived from lgamma_r() by Sun 01297 * 01298 * ==================================================== 01299 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 01300 * 01301 * Developed at SunPro, a Sun Microsystems, Inc. business. 01302 * Permission to use, copy, modify, and distribute this 01303 * software is freely granted, provided that this notice 01304 * is preserved. 01305 * ==================================================== 01306 * 01307 */ 01308 template <class REAL> 01309 REAL loggammaImpl(REAL x) 01310 { 01311 vigra_precondition(x > 0.0, 01312 "loggamma(): argument must be positive."); 01313 01314 vigra_precondition(x <= 1.0e307, 01315 "loggamma(): argument must not exceed 1e307."); 01316 01317 double res; 01318 01319 if (x < 4.2351647362715017e-22) 01320 { 01321 res = -std::log(x); 01322 } 01323 else if ((x == 2.0) || (x == 1.0)) 01324 { 01325 res = 0.0; 01326 } 01327 else if (x < 2.0) 01328 { 01329 static const double a[] = { 7.72156649015328655494e-02, 01330 3.22467033424113591611e-01, 01331 6.73523010531292681824e-02, 01332 2.05808084325167332806e-02, 01333 7.38555086081402883957e-03, 01334 2.89051383673415629091e-03, 01335 1.19270763183362067845e-03, 01336 5.10069792153511336608e-04, 01337 2.20862790713908385557e-04, 01338 1.08011567247583939954e-04, 01339 2.52144565451257326939e-05, 01340 4.48640949618915160150e-05 }; 01341 static const double t[] = { 4.83836122723810047042e-01, 01342 -1.47587722994593911752e-01, 01343 6.46249402391333854778e-02, 01344 -3.27885410759859649565e-02, 01345 1.79706750811820387126e-02, 01346 -1.03142241298341437450e-02, 01347 6.10053870246291332635e-03, 01348 -3.68452016781138256760e-03, 01349 2.25964780900612472250e-03, 01350 -1.40346469989232843813e-03, 01351 8.81081882437654011382e-04, 01352 -5.38595305356740546715e-04, 01353 3.15632070903625950361e-04, 01354 -3.12754168375120860518e-04, 01355 3.35529192635519073543e-04 }; 01356 static const double u[] = { -7.72156649015328655494e-02, 01357 6.32827064025093366517e-01, 01358 1.45492250137234768737e+00, 01359 9.77717527963372745603e-01, 01360 2.28963728064692451092e-01, 01361 1.33810918536787660377e-02 }; 01362 static const double v[] = { 0.0, 01363 2.45597793713041134822e+00, 01364 2.12848976379893395361e+00, 01365 7.69285150456672783825e-01, 01366 1.04222645593369134254e-01, 01367 3.21709242282423911810e-03 }; 01368 static const double tc = 1.46163214496836224576e+00; 01369 static const double tf = -1.21486290535849611461e-01; 01370 static const double tt = -3.63867699703950536541e-18; 01371 if (x <= 0.9) 01372 { 01373 res = -std::log(x); 01374 if (x >= 0.7316) 01375 { 01376 double y = 1.0-x; 01377 double z = y*y; 01378 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10])))); 01379 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11]))))); 01380 double p = y*p1+p2; 01381 res += (p-0.5*y); 01382 } 01383 else if (x >= 0.23164) 01384 { 01385 double y = x-(tc-1.0); 01386 double z = y*y; 01387 double w = z*y; 01388 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12]))); 01389 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13]))); 01390 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14]))); 01391 double p = z*p1-(tt-w*(p2+y*p3)); 01392 res += (tf + p); 01393 } 01394 else 01395 { 01396 double y = x; 01397 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5]))))); 01398 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5])))); 01399 res += (-0.5*y + p1/p2); 01400 } 01401 } 01402 else 01403 { 01404 res = 0.0; 01405 if (x >= 1.7316) 01406 { 01407 double y = 2.0-x; 01408 double z = y*y; 01409 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10])))); 01410 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11]))))); 01411 double p = y*p1+p2; 01412 res += (p-0.5*y); 01413 } 01414 else if(x >= 1.23164) 01415 { 01416 double y = x-tc; 01417 double z = y*y; 01418 double w = z*y; 01419 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12]))); 01420 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13]))); 01421 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14]))); 01422 double p = z*p1-(tt-w*(p2+y*p3)); 01423 res += (tf + p); 01424 } 01425 else 01426 { 01427 double y = x-1.0; 01428 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5]))))); 01429 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5])))); 01430 res += (-0.5*y + p1/p2); 01431 } 01432 } 01433 } 01434 else if(x < 8.0) 01435 { 01436 static const double s[] = { -7.72156649015328655494e-02, 01437 2.14982415960608852501e-01, 01438 3.25778796408930981787e-01, 01439 1.46350472652464452805e-01, 01440 2.66422703033638609560e-02, 01441 1.84028451407337715652e-03, 01442 3.19475326584100867617e-05 }; 01443 static const double r[] = { 0.0, 01444 1.39200533467621045958e+00, 01445 7.21935547567138069525e-01, 01446 1.71933865632803078993e-01, 01447 1.86459191715652901344e-02, 01448 7.77942496381893596434e-04, 01449 7.32668430744625636189e-06 }; 01450 double i = std::floor(x); 01451 double y = x-i; 01452 double p = y*(s[0]+y*(s[1]+y*(s[2]+y*(s[3]+y*(s[4]+y*(s[5]+y*s[6])))))); 01453 double q = 1.0+y*(r[1]+y*(r[2]+y*(r[3]+y*(r[4]+y*(r[5]+y*r[6]))))); 01454 res = 0.5*y+p/q; 01455 double z = 1.0; 01456 while (i > 2.0) 01457 { 01458 --i; 01459 z *= (y+i); 01460 } 01461 res += std::log(z); 01462 } 01463 else if (x < 2.8823037615171174e+17) 01464 { 01465 static const double w[] = { 4.18938533204672725052e-01, 01466 8.33333333333329678849e-02, 01467 -2.77777777728775536470e-03, 01468 7.93650558643019558500e-04, 01469 -5.95187557450339963135e-04, 01470 8.36339918996282139126e-04, 01471 -1.63092934096575273989e-03 }; 01472 double t = std::log(x); 01473 double z = 1.0/x; 01474 double y = z*z; 01475 double yy = w[0]+z*(w[1]+y*(w[2]+y*(w[3]+y*(w[4]+y*(w[5]+y*w[6]))))); 01476 res = (x-0.5)*(t-1.0)+yy; 01477 } 01478 else 01479 { 01480 res = x*(std::log(x) - 1.0); 01481 } 01482 01483 return res; 01484 } 01485 01486 01487 } // namespace detail 01488 01489 /*! The gamma function. 01490 01491 This function implements the algorithm from<br> 01492 Zhang and Jin: "Computation of Special Functions", John Wiley and Sons, 1996. 01493 01494 The argument must be <= 171.0 and cannot be zero or a negative integer. An 01495 exception is thrown when these conditions are violated. 01496 01497 <b>\#include</b> <vigra/mathutil.hxx><br> 01498 Namespace: vigra 01499 */ 01500 inline double gamma(double x) 01501 { 01502 return detail::gammaImpl(x); 01503 } 01504 01505 /*! The natural logarithm of the gamma function. 01506 01507 This function is based on a free implementation by Sun Microsystems, Inc., see 01508 <a href="http://www.sourceware.org/cgi-bin/cvsweb.cgi/~checkout~/src/newlib/libm/mathfp/er_lgamma.c?rev=1.6&content-type=text/plain&cvsroot=src">sourceware.org</a> archive. It can be removed once all compilers support the new C99 01509 math functions. 01510 01511 The argument must be positive and < 1e30. An exception is thrown when these conditions are violated. 01512 01513 <b>\#include</b> <vigra/mathutil.hxx><br> 01514 Namespace: vigra 01515 */ 01516 inline double loggamma(double x) 01517 { 01518 return detail::loggammaImpl(x); 01519 } 01520 01521 01522 namespace detail { 01523 01524 // both f1 and f2 are unsigned here 01525 template<class FPT> 01526 inline 01527 FPT safeFloatDivision( FPT f1, FPT f2 ) 01528 { 01529 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max() 01530 ? NumericTraits<FPT>::max() 01531 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) || 01532 f1 == NumericTraits<FPT>::zero() 01533 ? NumericTraits<FPT>::zero() 01534 : f1/f2; 01535 } 01536 01537 } // namespace detail 01538 01539 /*! Tolerance based floating-point comparison. 01540 01541 Check whether two floating point numbers are equal within the given tolerance. 01542 This is useful because floating point numbers that should be equal in theory are 01543 rarely exactly equal in practice. If the tolerance \a epsilon is not given, 01544 twice the machine epsilon is used. 01545 01546 <b>\#include</b> <vigra/mathutil.hxx><br> 01547 Namespace: vigra 01548 */ 01549 template <class T1, class T2> 01550 bool 01551 closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon) 01552 { 01553 typedef typename PromoteTraits<T1, T2>::Promote T; 01554 if(l == 0.0) 01555 return VIGRA_CSTD::fabs(r) <= epsilon; 01556 if(r == 0.0) 01557 return VIGRA_CSTD::fabs(l) <= epsilon; 01558 T diff = VIGRA_CSTD::fabs( l - r ); 01559 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) ); 01560 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) ); 01561 01562 return (d1 <= epsilon && d2 <= epsilon); 01563 } 01564 01565 template <class T1, class T2> 01566 inline bool closeAtTolerance(T1 l, T2 r) 01567 { 01568 typedef typename PromoteTraits<T1, T2>::Promote T; 01569 return closeAtTolerance(l, r, T(2.0) * NumericTraits<T>::epsilon()); 01570 } 01571 01572 //@} 01573 01574 } // namespace vigra 01575 01576 #endif /* VIGRA_MATHUTIL_HXX */
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|