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

details vigra/mathutil.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)