[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/mathutil.hxx VIGRA

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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (20 Sep 2011)