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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)