[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2005 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* The VIGRA Website is */ 00008 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00009 /* Please direct questions, bug reports, and contributions to */ 00010 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00011 /* vigra@informatik.uni-hamburg.de */ 00012 /* */ 00013 /* Permission is hereby granted, free of charge, to any person */ 00014 /* obtaining a copy of this software and associated documentation */ 00015 /* files (the "Software"), to deal in the Software without */ 00016 /* restriction, including without limitation the rights to use, */ 00017 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00018 /* sell copies of the Software, and to permit persons to whom the */ 00019 /* Software is furnished to do so, subject to the following */ 00020 /* conditions: */ 00021 /* */ 00022 /* The above copyright notice and this permission notice shall be */ 00023 /* included in all copies or substantial portions of the */ 00024 /* Software. */ 00025 /* */ 00026 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00027 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00028 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00029 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00030 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00031 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00032 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00033 /* OTHER DEALINGS IN THE SOFTWARE. */ 00034 /* */ 00035 /************************************************************************/ 00036 00037 #ifndef VIGRA_MATHUTIL_HXX 00038 #define VIGRA_MATHUTIL_HXX 00039 00040 #include <cmath> 00041 #include <cstdlib> 00042 #include "config.hxx" 00043 #include "error.hxx" 00044 #include "tuple.hxx" 00045 #include "sized_int.hxx" 00046 #include "numerictraits.hxx" 00047 00048 /*! \page MathConstants Mathematical Constants 00049 00050 <TT>M_PI, M_SQRT2</TT> 00051 00052 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>> 00053 00054 Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized, 00055 we provide definitions here for those compilers that don't support them. 00056 00057 \code 00058 #ifndef M_PI 00059 # define M_PI 3.14159265358979323846 00060 #endif 00061 00062 #ifndef M_SQRT2 00063 # define M_SQRT2 1.41421356237309504880 00064 #endif 00065 \endcode 00066 */ 00067 #ifndef M_PI 00068 # define M_PI 3.14159265358979323846 00069 #endif 00070 00071 #ifndef M_SQRT2 00072 # define M_SQRT2 1.41421356237309504880 00073 #endif 00074 00075 namespace vigra { 00076 00077 /** \addtogroup MathFunctions Mathematical Functions 00078 00079 Useful mathematical functions and functors. 00080 */ 00081 //@{ 00082 // import functions into namespace vigra which VIGRA is going to overload 00083 00084 using VIGRA_CSTD::pow; 00085 using VIGRA_CSTD::floor; 00086 using VIGRA_CSTD::ceil; 00087 00088 // import abs(float), abs(double), abs(long double) from <cmath> 00089 // and abs(int), abs(long), abs(long long) from <cstdlib> 00090 using std::abs; 00091 00092 // define the missing variants of abs() to avoid 'ambigous overload' 00093 // errors in template functions 00094 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \ 00095 inline T abs(T t) { return t; } 00096 00097 VIGRA_DEFINE_UNSIGNED_ABS(bool) 00098 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char) 00099 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short) 00100 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int) 00101 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long) 00102 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long long) 00103 00104 #undef VIGRA_DEFINE_UNSIGNED_ABS 00105 00106 #define VIGRA_DEFINE_MISSING_ABS(T) \ 00107 inline T abs(T t) { return t < 0 ? -t : t; } 00108 00109 VIGRA_DEFINE_MISSING_ABS(signed char) 00110 VIGRA_DEFINE_MISSING_ABS(signed short) 00111 00112 #undef VIGRA_DEFINE_MISSING_ABS 00113 00114 /*! The rounding function. 00115 00116 Defined for all floating point types. Rounds towards the nearest integer 00117 such that <tt>abs(round(t)) == round(abs(t))</tt> for all <tt>t</tt>. 00118 00119 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00120 Namespace: vigra 00121 */ 00122 inline float round(float t) 00123 { 00124 return t >= 0.0 00125 ? floor(t + 0.5) 00126 : ceil(t - 0.5); 00127 } 00128 00129 inline double round(double t) 00130 { 00131 return t >= 0.0 00132 ? floor(t + 0.5) 00133 : ceil(t - 0.5); 00134 } 00135 00136 inline long double round(long double t) 00137 { 00138 return t >= 0.0 00139 ? floor(t + 0.5) 00140 : ceil(t - 0.5); 00141 } 00142 00143 /*! Round up to the nearest power of 2. 00144 00145 Efficient algorithm for finding the smallest power of 2 which is not smaller than \a x 00146 (function clp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003, 00147 see http://www.hackersdelight.org/). 00148 If \a x > 2^31, the function will return 0 because integer arithmetic is defined modulo 2^32. 00149 00150 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00151 Namespace: vigra 00152 */ 00153 inline UInt32 ceilPower2(UInt32 x) 00154 { 00155 if(x == 0) return 0; 00156 00157 x = x - 1; 00158 x = x | (x >> 1); 00159 x = x | (x >> 2); 00160 x = x | (x >> 4); 00161 x = x | (x >> 8); 00162 x = x | (x >>16); 00163 return x + 1; 00164 } 00165 00166 /*! Round down to the nearest power of 2. 00167 00168 Efficient algorithm for finding the largest power of 2 which is not greater than \a x 00169 (function flp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003, 00170 see http://www.hackersdelight.org/). 00171 00172 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00173 Namespace: vigra 00174 */ 00175 inline UInt32 floorPower2(UInt32 x) 00176 { 00177 x = x | (x >> 1); 00178 x = x | (x >> 2); 00179 x = x | (x >> 4); 00180 x = x | (x >> 8); 00181 x = x | (x >>16); 00182 return x - (x >> 1); 00183 } 00184 00185 namespace detail { 00186 00187 template <class T> 00188 class IntLog2 00189 { 00190 public: 00191 static Int32 table[64]; 00192 }; 00193 00194 template <class T> 00195 Int32 IntLog2<T>::table[64] = { 00196 -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21, 00197 29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6, 00198 -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9, 00199 11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25, 00200 -1, 7, 24, -1, 23, -1, 31, -1}; 00201 00202 } // namespace detail 00203 00204 /*! Compute the base-2 logarithm of an integer. 00205 00206 Returns the position of the left-most 1-bit in the given number \a x, or 00207 -1 if \a x == 0. That is, 00208 00209 \code 00210 assert(k >= 0 && k < 32 && log2i(1 << k) == k); 00211 \endcode 00212 00213 The function uses Robert Harley's algorithm to determine the number of leading zeros 00214 in \a x (algorithm nlz10() at http://www.hackersdelight.org/). But note that the functions 00215 \ref floorPower2() or \ref ceilPower2() are more efficient and should be preferred when possible. 00216 00217 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00218 Namespace: vigra 00219 */ 00220 inline Int32 log2i(UInt32 x) 00221 { 00222 // Propagate leftmost 1-bit to the right. 00223 x = x | (x >> 1); 00224 x = x | (x >> 2); 00225 x = x | (x >> 4); 00226 x = x | (x >> 8); 00227 x = x | (x >>16); 00228 x = x*0x06EB14F9; // Multiplier is 7*255**3. 00229 return detail::IntLog2<Int32>::table[x >> 26]; 00230 } 00231 00232 /*! The square function. 00233 00234 <tt>sq(x) = x*x</tt> is needed so often that it makes sense to define it as a function. 00235 00236 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00237 Namespace: vigra 00238 */ 00239 template <class T> 00240 inline 00241 typename NumericTraits<T>::Promote sq(T t) 00242 { 00243 return t*t; 00244 } 00245 00246 namespace detail { 00247 00248 template <class T> 00249 class IntSquareRoot 00250 { 00251 public: 00252 static UInt32 sqq_table[]; 00253 static UInt32 exec(UInt32 v); 00254 }; 00255 00256 template <class T> 00257 UInt32 IntSquareRoot<T>::sqq_table[] = { 00258 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 00259 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 00260 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 00261 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118, 00262 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132, 00263 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145, 00264 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157, 00265 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, 00266 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 00267 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 00268 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 00269 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206, 00270 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215, 00271 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223, 00272 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231, 00273 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, 00274 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 00275 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 00276 253, 254, 254, 255 00277 }; 00278 00279 template <class T> 00280 UInt32 IntSquareRoot<T>::exec(UInt32 x) 00281 { 00282 UInt32 xn; 00283 if (x >= 0x10000) 00284 if (x >= 0x1000000) 00285 if (x >= 0x10000000) 00286 if (x >= 0x40000000) { 00287 if (x >= (UInt32)65535*(UInt32)65535) 00288 return 65535; 00289 xn = sqq_table[x>>24] << 8; 00290 } else 00291 xn = sqq_table[x>>22] << 7; 00292 else 00293 if (x >= 0x4000000) 00294 xn = sqq_table[x>>20] << 6; 00295 else 00296 xn = sqq_table[x>>18] << 5; 00297 else { 00298 if (x >= 0x100000) 00299 if (x >= 0x400000) 00300 xn = sqq_table[x>>16] << 4; 00301 else 00302 xn = sqq_table[x>>14] << 3; 00303 else 00304 if (x >= 0x40000) 00305 xn = sqq_table[x>>12] << 2; 00306 else 00307 xn = sqq_table[x>>10] << 1; 00308 00309 goto nr1; 00310 } 00311 else 00312 if (x >= 0x100) { 00313 if (x >= 0x1000) 00314 if (x >= 0x4000) 00315 xn = (sqq_table[x>>8] >> 0) + 1; 00316 else 00317 xn = (sqq_table[x>>6] >> 1) + 1; 00318 else 00319 if (x >= 0x400) 00320 xn = (sqq_table[x>>4] >> 2) + 1; 00321 else 00322 xn = (sqq_table[x>>2] >> 3) + 1; 00323 00324 goto adj; 00325 } else 00326 return sqq_table[x] >> 4; 00327 00328 /* Run two iterations of the standard convergence formula */ 00329 00330 xn = (xn + 1 + x / xn) / 2; 00331 nr1: 00332 xn = (xn + 1 + x / xn) / 2; 00333 adj: 00334 00335 if (xn * xn > x) /* Correct rounding if necessary */ 00336 xn--; 00337 00338 return xn; 00339 } 00340 00341 } // namespace detail 00342 00343 using VIGRA_CSTD::sqrt; 00344 00345 /*! Signed integer square root. 00346 00347 Useful for fast fixed-point computations. 00348 00349 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00350 Namespace: vigra 00351 */ 00352 inline Int32 sqrti(Int32 v) 00353 { 00354 if(v < 0) 00355 throw std::domain_error("sqrti(Int32): negative argument."); 00356 return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v); 00357 } 00358 00359 /*! Unsigned integer square root. 00360 00361 Useful for fast fixed-point computations. 00362 00363 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00364 Namespace: vigra 00365 */ 00366 inline UInt32 sqrti(UInt32 v) 00367 { 00368 return detail::IntSquareRoot<UInt32>::exec(v); 00369 } 00370 00371 #ifndef VIGRA_HAS_HYPOT 00372 /*! Compute the Euclidean distance (length of the hypothenuse of a right-angled triangle). 00373 00374 The hypot() function returns the sqrt(a*a + b*b). 00375 It is implemented in a way that minimizes round-off error. 00376 00377 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00378 Namespace: vigra 00379 */ 00380 inline double hypot(double a, double b) 00381 { 00382 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b); 00383 if (absa > absb) 00384 return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 00385 else 00386 return absb == 0.0 00387 ? 0.0 00388 : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 00389 } 00390 00391 #else 00392 00393 using ::hypot; 00394 00395 #endif 00396 00397 /*! The sign function. 00398 00399 Returns 1, 0, or -1 depending on the sign of \a t. 00400 00401 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00402 Namespace: vigra 00403 */ 00404 template <class T> 00405 inline T sign(T t) 00406 { 00407 return t > NumericTraits<T>::zero() 00408 ? NumericTraits<T>::one() 00409 : t < NumericTraits<T>::zero() 00410 ? -NumericTraits<T>::one() 00411 : NumericTraits<T>::zero(); 00412 } 00413 00414 /*! The binary sign function. 00415 00416 Transfers the sign of \a t2 to \a t1. 00417 00418 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00419 Namespace: vigra 00420 */ 00421 template <class T1, class T2> 00422 inline T1 sign(T1 t1, T2 t2) 00423 { 00424 return t2 >= NumericTraits<T2>::zero() 00425 ? abs(t1) 00426 : -abs(t1); 00427 } 00428 00429 #define VIGRA_DEFINE_NORM(T) \ 00430 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \ 00431 inline NormTraits<T>::NormType norm(T t) { return abs(t); } 00432 00433 VIGRA_DEFINE_NORM(bool) 00434 VIGRA_DEFINE_NORM(signed char) 00435 VIGRA_DEFINE_NORM(unsigned char) 00436 VIGRA_DEFINE_NORM(short) 00437 VIGRA_DEFINE_NORM(unsigned short) 00438 VIGRA_DEFINE_NORM(int) 00439 VIGRA_DEFINE_NORM(unsigned int) 00440 VIGRA_DEFINE_NORM(long) 00441 VIGRA_DEFINE_NORM(unsigned long) 00442 VIGRA_DEFINE_NORM(float) 00443 VIGRA_DEFINE_NORM(double) 00444 VIGRA_DEFINE_NORM(long double) 00445 00446 #undef VIGRA_DEFINE_NORM 00447 00448 template <class T> 00449 inline typename NormTraits<std::complex<T> >::SquaredNormType 00450 squaredNorm(std::complex<T> const & t) 00451 { 00452 return sq(t.real()) + sq(t.imag()); 00453 } 00454 00455 #ifdef DOXYGEN // only for documentation 00456 /*! The squared norm of a numerical object. 00457 00458 For scalar types: equals <tt>vigra::sq(t)</tt><br>. 00459 For vectorial types: equals <tt>vigra::dot(t, t)</tt><br>. 00460 For complex types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt><br>. 00461 For matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements). 00462 */ 00463 NormTraits<T>::SquaredNormType squaredNorm(T const & t); 00464 00465 #endif 00466 00467 /*! The norm of a numerical object. 00468 00469 For scalar types: implemented as <tt>abs(t)</tt><br> 00470 otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>. 00471 00472 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00473 Namespace: vigra 00474 */ 00475 template <class T> 00476 inline typename NormTraits<T>::NormType 00477 norm(T const & t) 00478 { 00479 typedef typename NormTraits<T>::SquaredNormType SNT; 00480 return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t))); 00481 } 00482 00483 /*! Find the minimum element in a sequence. 00484 00485 The function returns the iterator refering to the minimum element. 00486 00487 <b>Required Interface:</b> 00488 00489 \code 00490 Iterator is a standard forward iterator. 00491 00492 bool f = *first < NumericTraits<typename std::iterator_traits<Iterator>::value_type>::max(); 00493 \endcode 00494 00495 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00496 Namespace: vigra 00497 */ 00498 template <class Iterator> 00499 Iterator argMin(Iterator first, Iterator last) 00500 { 00501 typedef typename std::iterator_traits<Iterator>::value_type Value; 00502 Value vopt = NumericTraits<Value>::max(); 00503 Iterator best = last; 00504 for(; first != last; ++first) 00505 { 00506 if(*first < vopt) 00507 { 00508 vopt = *first; 00509 best = first; 00510 } 00511 } 00512 return best; 00513 } 00514 00515 /*! Find the maximum element in a sequence. 00516 00517 The function returns the iterator refering to the maximum element. 00518 00519 <b>Required Interface:</b> 00520 00521 \code 00522 Iterator is a standard forward iterator. 00523 00524 bool f = NumericTraits<typename std::iterator_traits<Iterator>::value_type>::min() < *first; 00525 \endcode 00526 00527 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00528 Namespace: vigra 00529 */ 00530 template <class Iterator> 00531 Iterator argMax(Iterator first, Iterator last) 00532 { 00533 typedef typename std::iterator_traits<Iterator>::value_type Value; 00534 Value vopt = NumericTraits<Value>::min(); 00535 Iterator best = last; 00536 for(; first != last; ++first) 00537 { 00538 if(vopt < *first) 00539 { 00540 vopt = *first; 00541 best = first; 00542 } 00543 } 00544 return best; 00545 } 00546 00547 /*! Find the minimum element in a sequence conforming to a condition. 00548 00549 The function returns the iterator refering to the minimum element, 00550 where only elements conforming to the condition (i.e. where 00551 <tt>condition(*iterator)</tt> evaluates to <tt>true</tt>) are considered. 00552 If no element conforms to the condition, or the sequence is empty, 00553 the end iterator \a last is returned. 00554 00555 <b>Required Interface:</b> 00556 00557 \code 00558 Iterator is a standard forward iterator. 00559 00560 bool c = condition(*first); 00561 00562 bool f = *first < NumericTraits<typename std::iterator_traits<Iterator>::value_type>::max(); 00563 \endcode 00564 00565 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00566 Namespace: vigra 00567 */ 00568 template <class Iterator, class UnaryFunctor> 00569 Iterator argMinIf(Iterator first, Iterator last, UnaryFunctor condition) 00570 { 00571 typedef typename std::iterator_traits<Iterator>::value_type Value; 00572 Value vopt = NumericTraits<Value>::max(); 00573 Iterator best = last; 00574 for(; first != last; ++first) 00575 { 00576 if(condition(*first) && *first < vopt) 00577 { 00578 vopt = *first; 00579 best = first; 00580 } 00581 } 00582 return best; 00583 } 00584 00585 /*! Find the maximum element in a sequence conforming to a condition. 00586 00587 The function returns the iterator refering to the maximum element, 00588 where only elements conforming to the condition (i.e. where 00589 <tt>condition(*iterator)</tt> evaluates to <tt>true</tt>) are considered. 00590 If no element conforms to the condition, or the sequence is empty, 00591 the end iterator \a last is returned. 00592 00593 <b>Required Interface:</b> 00594 00595 \code 00596 Iterator is a standard forward iterator. 00597 00598 bool c = condition(*first); 00599 00600 bool f = NumericTraits<typename std::iterator_traits<Iterator>::value_type>::min() < *first; 00601 \endcode 00602 00603 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00604 Namespace: vigra 00605 */ 00606 template <class Iterator, class UnaryFunctor> 00607 Iterator argMaxIf(Iterator first, Iterator last, UnaryFunctor condition) 00608 { 00609 typedef typename std::iterator_traits<Iterator>::value_type Value; 00610 Value vopt = NumericTraits<Value>::min(); 00611 Iterator best = last; 00612 for(; first != last; ++first) 00613 { 00614 if(condition(*first) && vopt < *first) 00615 { 00616 vopt = *first; 00617 best = first; 00618 } 00619 } 00620 return best; 00621 } 00622 00623 00624 namespace detail { 00625 00626 template <class T> 00627 T ellipticRD(T x, T y, T z) 00628 { 00629 double f = 1.0, s = 0.0, X, Y, Z, m; 00630 for(;;) 00631 { 00632 m = (x + y + 3.0*z) / 5.0; 00633 X = 1.0 - x/m; 00634 Y = 1.0 - y/m; 00635 Z = 1.0 - z/m; 00636 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01) 00637 break; 00638 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z); 00639 s += f / (VIGRA_CSTD::sqrt(z)*(z + l)); 00640 f /= 4.0; 00641 x = (x + l)/4.0; 00642 y = (y + l)/4.0; 00643 z = (z + l)/4.0; 00644 } 00645 double a = X*Y; 00646 double b = sq(Z); 00647 double c = a - b; 00648 double d = a - 6.0*b; 00649 double e = d + 2.0*c; 00650 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0) 00651 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5); 00652 } 00653 00654 template <class T> 00655 T ellipticRF(T x, T y, T z) 00656 { 00657 double X, Y, Z, m; 00658 for(;;) 00659 { 00660 m = (x + y + z) / 3.0; 00661 X = 1.0 - x/m; 00662 Y = 1.0 - y/m; 00663 Z = 1.0 - z/m; 00664 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01) 00665 break; 00666 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z); 00667 x = (x + l)/4.0; 00668 y = (y + l)/4.0; 00669 z = (z + l)/4.0; 00670 } 00671 double d = X*Y - sq(Z); 00672 double p = X*Y*Z; 00673 return (1.0 - d/10.0 + p/14.0 + sq(d)/24.0 - d*p*3.0/44.0) / VIGRA_CSTD::sqrt(m); 00674 } 00675 00676 } // namespace detail 00677 00678 /*! The incomplete elliptic integral of the first kind. 00679 00680 Computes 00681 00682 \f[ 00683 \mbox{F}(x, k) = \int_0^x \frac{1}{\sqrt{1 - k^2 \sin(t)^2}} dt 00684 \f] 00685 00686 according to the algorithm given in Press et al. "Numerical Recipes". 00687 00688 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral 00689 functions must be k^2 rather than k. Check the documentation when results disagree! 00690 00691 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00692 Namespace: vigra 00693 */ 00694 inline double ellipticIntegralF(double x, double k) 00695 { 00696 double c2 = sq(VIGRA_CSTD::cos(x)); 00697 double s = VIGRA_CSTD::sin(x); 00698 return s*detail::ellipticRF(c2, 1.0 - sq(k*s), 1.0); 00699 } 00700 00701 /*! The incomplete elliptic integral of the second kind. 00702 00703 Computes 00704 00705 \f[ 00706 \mbox{E}(x, k) = \int_0^x \sqrt{1 - k^2 \sin(t)^2} dt 00707 \f] 00708 00709 according to the algorithm given in Press et al. "Numerical Recipes". The 00710 complete elliptic integral of the second kind is simply <tt>ellipticIntegralE(M_PI/2, k)</TT>. 00711 00712 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral 00713 functions must be k^2 rather than k. Check the documentation when results disagree! 00714 00715 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00716 Namespace: vigra 00717 */ 00718 inline double ellipticIntegralE(double x, double k) 00719 { 00720 double c2 = sq(VIGRA_CSTD::cos(x)); 00721 double s = VIGRA_CSTD::sin(x); 00722 k = sq(k*s); 00723 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0)); 00724 } 00725 00726 #ifndef VIGRA_HAS_ERF 00727 00728 namespace detail { 00729 00730 template <class T> 00731 double erfImpl(T x) 00732 { 00733 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x)); 00734 double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+ 00735 t*(0.09678418+t*(-0.18628806+t*(0.27886807+ 00736 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+ 00737 t*0.17087277))))))))); 00738 if (x >= 0.0) 00739 return 1.0 - ans; 00740 else 00741 return ans - 1.0; 00742 } 00743 00744 } // namespace detail 00745 00746 /*! The error function. 00747 00748 If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the 00749 new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error 00750 function 00751 00752 \f[ 00753 \mbox{erf}(x) = \int_0^x e^{-t^2} dt 00754 \f] 00755 00756 according to the formula given in Press et al. "Numerical Recipes". 00757 00758 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00759 Namespace: vigra 00760 */ 00761 inline double erf(double x) 00762 { 00763 return detail::erfImpl(x); 00764 } 00765 00766 #else 00767 00768 using VIGRA_CSTD::erf; 00769 00770 #endif 00771 00772 namespace detail { 00773 00774 template <class T> 00775 double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, T noncentrality, T arg) 00776 { 00777 double a = degreesOfFreedom + noncentrality; 00778 double b = (a + noncentrality) / sq(a); 00779 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); 00780 return 0.5*(1.0 + erf(t/VIGRA_CSTD::sqrt(2.0))); 00781 } 00782 00783 template <class T> 00784 void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j) 00785 { 00786 double tol = -50.0; 00787 if(lans < tol) 00788 { 00789 lans = lans + VIGRA_CSTD::log(arg / j); 00790 dans = VIGRA_CSTD::exp(lans); 00791 } 00792 else 00793 { 00794 dans = dans * arg / j; 00795 } 00796 pans = pans - dans; 00797 j += 2; 00798 } 00799 00800 template <class T> 00801 std::pair<double, double> noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps) 00802 { 00803 vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0, 00804 "noncentralChi2P(): parameters must be positive."); 00805 if (arg == 0.0 && degreesOfFreedom > 0) 00806 return std::make_pair(0.0, 0.0); 00807 00808 // Determine initial values 00809 double b1 = 0.5 * noncentrality, 00810 ao = VIGRA_CSTD::exp(-b1), 00811 eps2 = eps / ao, 00812 lnrtpi2 = 0.22579135264473, 00813 probability, density, lans, dans, pans, sum, am, hold; 00814 unsigned int maxit = 500, 00815 i, m; 00816 if(degreesOfFreedom % 2) 00817 { 00818 i = 1; 00819 lans = -0.5 * (arg + VIGRA_CSTD::log(arg)) - lnrtpi2; 00820 dans = VIGRA_CSTD::exp(lans); 00821 pans = erf(VIGRA_CSTD::sqrt(arg/2.0)); 00822 } 00823 else 00824 { 00825 i = 2; 00826 lans = -0.5 * arg; 00827 dans = VIGRA_CSTD::exp(lans); 00828 pans = 1.0 - dans; 00829 } 00830 00831 // Evaluate first term 00832 if(degreesOfFreedom == 0) 00833 { 00834 m = 1; 00835 degreesOfFreedom = 2; 00836 am = b1; 00837 sum = 1.0 / ao - 1.0 - am; 00838 density = am * dans; 00839 probability = 1.0 + am * pans; 00840 } 00841 else 00842 { 00843 m = 0; 00844 degreesOfFreedom = degreesOfFreedom - 1; 00845 am = 1.0; 00846 sum = 1.0 / ao - 1.0; 00847 while(i < degreesOfFreedom) 00848 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i); 00849 degreesOfFreedom = degreesOfFreedom + 1; 00850 density = dans; 00851 probability = pans; 00852 } 00853 // Evaluate successive terms of the expansion 00854 for(++m; m<maxit; ++m) 00855 { 00856 am = b1 * am / m; 00857 detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom); 00858 sum = sum - am; 00859 density = density + am * dans; 00860 hold = am * pans; 00861 probability = probability + hold; 00862 if((pans * sum < eps2) && (hold < eps2)) 00863 break; // converged 00864 } 00865 if(m == maxit) 00866 vigra_fail("noncentralChi2P(): no convergence."); 00867 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability))); 00868 } 00869 00870 } // namespace detail 00871 00872 /*! Chi square distribution. 00873 00874 Computes the density of a chi square distribution with \a degreesOfFreedom 00875 and tolerance \a accuracy at the given argument \a arg 00876 by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>. 00877 00878 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00879 Namespace: vigra 00880 */ 00881 inline double chi2(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7) 00882 { 00883 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).first; 00884 } 00885 00886 /*! Cumulative chi square distribution. 00887 00888 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom 00889 and tolerance \a accuracy at the given argument \a arg, i.e. the probability that 00890 a random number drawn from the distribution is below \a arg 00891 by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>. 00892 00893 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00894 Namespace: vigra 00895 */ 00896 inline double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7) 00897 { 00898 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).second; 00899 } 00900 00901 /*! Non-central chi square distribution. 00902 00903 Computes the density of a non-central chi square distribution with \a degreesOfFreedom, 00904 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 00905 \a arg. It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 00906 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of 00907 degrees of freedom. 00908 00909 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00910 Namespace: vigra 00911 */ 00912 inline double noncentralChi2(unsigned int degreesOfFreedom, 00913 double noncentrality, double arg, double accuracy = 1e-7) 00914 { 00915 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).first; 00916 } 00917 00918 /*! Cumulative non-central chi square distribution. 00919 00920 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom, 00921 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument 00922 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg 00923 It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from 00924 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of 00925 degrees of freedom (see noncentralChi2CDFApprox() for a constant-time algorithm). 00926 00927 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00928 Namespace: vigra 00929 */ 00930 inline double noncentralChi2CDF(unsigned int degreesOfFreedom, 00931 double noncentrality, double arg, double accuracy = 1e-7) 00932 { 00933 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).second; 00934 } 00935 00936 /*! Cumulative non-central chi square distribution (approximate). 00937 00938 Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom, 00939 and noncentrality parameter \a noncentrality at the given argument 00940 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg 00941 It uses the approximate transform into a normal distribution due to Wilson and Hilferty 00942 (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32). 00943 The algorithm's running time is independent of the inputs, i.e. is should be used 00944 when noncentralChi2CDF() is too slow, and approximate values are sufficient. The accuracy is only 00945 about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5. 00946 00947 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00948 Namespace: vigra 00949 */ 00950 inline double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg) 00951 { 00952 return detail::noncentralChi2CDFApprox(degreesOfFreedom, noncentrality, arg); 00953 } 00954 00955 00956 00957 namespace detail { 00958 00959 // both f1 and f2 are unsigned here 00960 template<class FPT> 00961 inline 00962 FPT safeFloatDivision( FPT f1, FPT f2 ) 00963 { 00964 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max() 00965 ? NumericTraits<FPT>::max() 00966 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) || 00967 f1 == NumericTraits<FPT>::zero() 00968 ? NumericTraits<FPT>::zero() 00969 : f1/f2; 00970 } 00971 00972 } // namespace detail 00973 00974 /*! Tolerance based floating-point comparison. 00975 00976 Check whether two floating point numbers are equal within the given tolerance. 00977 This is useful because floating point numbers that should be equal in theory are 00978 rarely exactly equal in practice. If the tolerance \a epsilon is not given, 00979 twice the machine epsilon is used. 00980 00981 <b>\#include</b> <<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>><br> 00982 Namespace: vigra 00983 */ 00984 template <class T1, class T2> 00985 bool 00986 closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon) 00987 { 00988 typedef typename PromoteTraits<T1, T2>::Promote T; 00989 if(l == 0.0) 00990 return VIGRA_CSTD::fabs(r) <= epsilon; 00991 if(r == 0.0) 00992 return VIGRA_CSTD::fabs(l) <= epsilon; 00993 T diff = VIGRA_CSTD::fabs( l - r ); 00994 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) ); 00995 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) ); 00996 00997 return (d1 <= epsilon && d2 <= epsilon); 00998 } 00999 01000 template <class T1, class T2> 01001 inline bool closeAtTolerance(T1 l, T2 r) 01002 { 01003 typedef typename PromoteTraits<T1, T2>::Promote T; 01004 return closeAtTolerance(l, r, 2.0 * NumericTraits<T>::epsilon()); 01005 } 01006 01007 //@} 01008 01009 } // namespace vigra 01010 01011 #endif /* VIGRA_MATHUTIL_HXX */
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|