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

vigra/mathutil.hxx |

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2011 by Ullrich Koethe                  */
    00004 /*                                                                      */
    00005 /*    This file is part of the VIGRA computer vision library.           */
    00006 /*    The VIGRA Website is                                              */
    00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
    00008 /*    Please direct questions, bug reports, and contributions to        */
    00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00010 /*        vigra@informatik.uni-hamburg.de                               */
    00011 /*                                                                      */
    00012 /*    Permission is hereby granted, free of charge, to any person       */
    00013 /*    obtaining a copy of this software and associated documentation    */
    00014 /*    files (the "Software"), to deal in the Software without           */
    00015 /*    restriction, including without limitation the rights to use,      */
    00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00017 /*    sell copies of the Software, and to permit persons to whom the    */
    00018 /*    Software is furnished to do so, subject to the following          */
    00019 /*    conditions:                                                       */
    00020 /*                                                                      */
    00021 /*    The above copyright notice and this permission notice shall be    */
    00022 /*    included in all copies or substantial portions of the             */
    00023 /*    Software.                                                         */
    00024 /*                                                                      */
    00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00033 /*                                                                      */
    00034 /************************************************************************/
    00035
    00036 #ifndef VIGRA_MATHUTIL_HXX
    00037 #define VIGRA_MATHUTIL_HXX
    00038
    00039 #ifdef _MSC_VER
    00040 # pragma warning (disable: 4996) // hypot/_hypot confusion
    00041 #endif
    00042
    00043 #include <cmath>
    00044 #include <cstdlib>
    00045 #include <complex>
    00046 #include "config.hxx"
    00047 #include "error.hxx"
    00048 #include "tuple.hxx"
    00049 #include "sized_int.hxx"
    00050 #include "numerictraits.hxx"
    00051 #include "algorithm.hxx"
    00052
    00053 /*! \page MathConstants Mathematical Constants
    00054
    00055     <TT>M_PI, M_SQRT2 etc.</TT>
    00056
    00057     <b>\#include</b> <vigra/mathutil.hxx>
    00058
    00059     Since mathematical constants such as <TT>M_PI</TT> and <TT>M_SQRT2</TT>
    00060     are not officially standardized, we provide definitions here for those
    00061     compilers that don't support them.
    00062
    00063     \code
    00064     #ifndef M_PI
    00065     #    define M_PI     3.14159265358979323846
    00066     #endif
    00067
    00068     #ifndef M_SQRT2
    00069     #    define M_2_PI   0.63661977236758134308
    00070     #endif
    00071
    00072     #ifndef M_PI_2
    00073     #    define M_PI_2   1.57079632679489661923
    00074     #endif
    00075
    00076     #ifndef M_PI_4
    00077     #    define M_PI_4   0.78539816339744830962
    00078     #endif
    00079
    00080     #ifndef M_SQRT2
    00081     #    define M_SQRT2  1.41421356237309504880
    00082     #endif
    00083
    00084     #ifndef M_EULER_GAMMA
    00085     #    define M_EULER_GAMMA  0.5772156649015329
    00086     #endif
    00087     \endcode
    00088 */
    00089 #ifndef M_PI
    00090 #    define M_PI     3.14159265358979323846
    00091 #endif
    00092
    00093 #ifndef M_2_PI
    00094 #    define M_2_PI   0.63661977236758134308
    00095 #endif
    00096
    00097 #ifndef M_PI_2
    00098 #    define M_PI_2   1.57079632679489661923
    00099 #endif
    00100
    00101 #ifndef M_PI_4
    00102 #    define M_PI_4   0.78539816339744830962
    00103 #endif
    00104
    00105 #ifndef M_SQRT2
    00106 #    define M_SQRT2  1.41421356237309504880
    00107 #endif
    00108
    00109 #ifndef M_EULER_GAMMA
    00110 #    define M_EULER_GAMMA  0.5772156649015329
    00111 #endif
    00112
    00113 namespace vigra {
    00114
    00115 /** \addtogroup MathFunctions Mathematical Functions
    00116
    00117     Useful mathematical functions and functors.
    00118 */
    00119 //@{
    00120 // import functions into namespace vigra which VIGRA is going to overload
    00121
    00122 using VIGRA_CSTD::pow;
    00123 using VIGRA_CSTD::floor;
    00124 using VIGRA_CSTD::ceil;
    00125 using VIGRA_CSTD::exp;
    00126
    00127 // import abs(float), abs(double), abs(long double) from <cmath>
    00128 //        abs(int), abs(long), abs(long long) from <cstdlib>
    00129 //        abs(std::complex<T>) from <complex>
    00130 using std::abs;
    00131
    00132 // define the missing variants of abs() to avoid 'ambiguous overload'
    00133 // errors in template functions
    00134 #define VIGRA_DEFINE_UNSIGNED_ABS(T) \
    00135     inline T abs(T t) { return t; }
    00136
    00137 VIGRA_DEFINE_UNSIGNED_ABS(bool)
    00138 VIGRA_DEFINE_UNSIGNED_ABS(unsigned char)
    00139 VIGRA_DEFINE_UNSIGNED_ABS(unsigned short)
    00140 VIGRA_DEFINE_UNSIGNED_ABS(unsigned int)
    00141 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long)
    00142 VIGRA_DEFINE_UNSIGNED_ABS(unsigned long long)
    00143
    00144 #undef VIGRA_DEFINE_UNSIGNED_ABS
    00145
    00146 #define VIGRA_DEFINE_MISSING_ABS(T) \
    00147     inline T abs(T t) { return t < 0 ? 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)