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