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