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

mathutil.hxx
1/************************************************************************/
2/* */
3/* Copyright 1998-2011 by Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36#ifndef VIGRA_MATHUTIL_HXX
37#define VIGRA_MATHUTIL_HXX
38
39#ifdef _MSC_VER
40# pragma warning (disable: 4996) // hypot/_hypot confusion
41#endif
42
43#include <cmath>
44#include <cstdlib>
45#include <complex>
46#include "config.hxx"
47#include "error.hxx"
48#include "tuple.hxx"
49#include "sized_int.hxx"
50#include "numerictraits.hxx"
51#include "algorithm.hxx"
52
53/** \page MathConstants Mathematical Constants
54
55 <TT>M_PI, M_SQRT2 etc.</TT>
56
57 <b>\#include</b> <vigra/mathutil.hxx>
58
59 Since mathematical constants such as <TT>M_PI</TT> and <TT>M_SQRT2</TT>
60 are not officially standardized, we provide definitions here for those
61 compilers that don't support them.
62
63 \code
64 #ifndef M_PI
65 # define M_PI 3.14159265358979323846
66 #endif
67
68 #ifndef M_SQRT2
69 # define M_2_PI 0.63661977236758134308
70 #endif
71
72 #ifndef M_PI_2
73 # define M_PI_2 1.57079632679489661923
74 #endif
75
76 #ifndef M_PI_4
77 # define M_PI_4 0.78539816339744830962
78 #endif
79
80 #ifndef M_SQRT2
81 # define M_SQRT2 1.41421356237309504880
82 #endif
83
84 #ifndef M_EULER_GAMMA
85 # define M_EULER_GAMMA 0.5772156649015329
86 #endif
87 \endcode
88*/
89#ifndef M_PI
90# define M_PI 3.14159265358979323846
91#endif
92
93#ifndef M_2_PI
94# define M_2_PI 0.63661977236758134308
95#endif
96
97#ifndef M_PI_2
98# define M_PI_2 1.57079632679489661923
99#endif
100
101#ifndef M_PI_4
102# define M_PI_4 0.78539816339744830962
103#endif
104
105#ifndef M_SQRT2
106# define M_SQRT2 1.41421356237309504880
107#endif
108
109#ifndef M_E
110# define M_E 2.71828182845904523536
111#endif
112
113#ifndef M_EULER_GAMMA
114# define M_EULER_GAMMA 0.5772156649015329
115#endif
116
117namespace vigra {
118
119/** \addtogroup MathFunctions Mathematical Functions
120
121 Useful mathematical functions and functors.
122*/
123//@{
124// import functions into namespace vigra which VIGRA is going to overload
125
126using VIGRA_CSTD::pow;
127using VIGRA_CSTD::floor;
128using VIGRA_CSTD::ceil;
129using VIGRA_CSTD::exp;
130
131// import abs(float), abs(double), abs(long double) from <cmath>
132// abs(int), abs(long), abs(long long) from <cstdlib>
133// abs(std::complex<T>) from <complex>
134using std::abs;
135
136// define the missing variants of abs() to avoid 'ambiguous overload'
137// errors in template functions
138#define VIGRA_DEFINE_UNSIGNED_ABS(T) \
139 inline T abs(T t) { return t; }
140
141VIGRA_DEFINE_UNSIGNED_ABS(bool)
142VIGRA_DEFINE_UNSIGNED_ABS(unsigned char)
143VIGRA_DEFINE_UNSIGNED_ABS(unsigned short)
144VIGRA_DEFINE_UNSIGNED_ABS(unsigned int)
145VIGRA_DEFINE_UNSIGNED_ABS(unsigned long)
146VIGRA_DEFINE_UNSIGNED_ABS(unsigned long long)
147
148#undef VIGRA_DEFINE_UNSIGNED_ABS
149
150#define VIGRA_DEFINE_MISSING_ABS(T) \
151 inline T abs(T t) { return t < 0 ? static_cast<T>(-t) : t; }
152
153VIGRA_DEFINE_MISSING_ABS(signed char)
154VIGRA_DEFINE_MISSING_ABS(signed short)
155
156#if defined(_MSC_VER) && _MSC_VER < 1600
157VIGRA_DEFINE_MISSING_ABS(signed long long)
158#endif
159
160#undef VIGRA_DEFINE_MISSING_ABS
161
162#ifndef _MSC_VER
163
164using std::isinf;
165using std::isnan;
166using std::isfinite;
167
168#else
169
170template <class REAL>
171inline bool isinf(REAL v)
172{
173 return _finite(v) == 0 && _isnan(v) == 0;
174}
175
176template <class REAL>
177inline bool isnan(REAL v)
178{
179 return _isnan(v) != 0;
180}
181
182template <class REAL>
183inline bool isfinite(REAL v)
184{
185 return _finite(v) != 0;
186}
187
188#endif
189
190// scalar dot is needed for generic functions that should work with
191// scalars and vectors alike
192
193#define VIGRA_DEFINE_SCALAR_DOT(T) \
194 inline NumericTraits<T>::Promote dot(T l, T r) { return l*r; }
195
196VIGRA_DEFINE_SCALAR_DOT(unsigned char)
197VIGRA_DEFINE_SCALAR_DOT(unsigned short)
198VIGRA_DEFINE_SCALAR_DOT(unsigned int)
199VIGRA_DEFINE_SCALAR_DOT(unsigned long)
200VIGRA_DEFINE_SCALAR_DOT(unsigned long long)
201VIGRA_DEFINE_SCALAR_DOT(signed char)
202VIGRA_DEFINE_SCALAR_DOT(signed short)
203VIGRA_DEFINE_SCALAR_DOT(signed int)
204VIGRA_DEFINE_SCALAR_DOT(signed long)
205VIGRA_DEFINE_SCALAR_DOT(signed long long)
206VIGRA_DEFINE_SCALAR_DOT(float)
207VIGRA_DEFINE_SCALAR_DOT(double)
208VIGRA_DEFINE_SCALAR_DOT(long double)
209
210#undef VIGRA_DEFINE_SCALAR_DOT
211
212using std::pow;
213
214// support 'double' exponents for all floating point versions of pow()
215
216inline float pow(float v, double e)
217{
218 return std::pow(v, (float)e);
219}
220
221inline long double pow(long double v, double e)
222{
223 return std::pow(v, (long double)e);
224}
225
226 /** \brief The rounding function.
227
228 Defined for all floating point types. Rounds towards the nearest integer
229 such that <tt>abs(round(t)) == round(abs(t))</tt> for all <tt>t</tt>.
230
231 <b>\#include</b> <vigra/mathutil.hxx><br>
232 Namespace: vigra
233 */
234#ifdef DOXYGEN // only for documentation
236#endif
237
238inline float round(float t)
239{
240 return t >= 0.0
241 ? floor(t + 0.5f)
242 : ceil(t - 0.5f);
243}
244
245inline double round(double t)
246{
247 return t >= 0.0
248 ? floor(t + 0.5)
249 : ceil(t - 0.5);
250}
251
252inline long double round(long double t)
253{
254 return t >= 0.0
255 ? floor(t + 0.5)
256 : ceil(t - 0.5);
257}
258
259
260 /** \brief Round and cast to integer.
261
262 Rounds to the nearest integer like round(), but casts the result to
263 <tt>long long</tt> (this will be faster and is usually needed anyway).
264
265 <b>\#include</b> <vigra/mathutil.hxx><br>
266 Namespace: vigra
267 */
268inline long long roundi(double t)
269{
270 return t >= 0.0
271 ? (long long)(t + 0.5)
272 : (long long)(t - 0.5);
273}
274
275 /** \brief Determine whether x is a power of 2
276 Bit twiddle from https://graphics.stanford.edu/~seander/bithacks.html#DetermineIfPowerOf2
277 */
278inline bool isPower2(UInt32 x)
279{
280 return x && !(x & (x - 1));
281}
282
283
284 /** \brief Round up to the nearest power of 2.
285
286 Efficient algorithm for finding the smallest power of 2 which is not smaller than \a x
287 (function clp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003,
288 see http://www.hackersdelight.org/).
289 If \a x > 2^31, the function will return 0 because integer arithmetic is defined modulo 2^32.
290
291 <b>\#include</b> <vigra/mathutil.hxx><br>
292 Namespace: vigra
293 */
295{
296 if(x == 0) return 0;
297
298 x = x - 1;
299 x = x | (x >> 1);
300 x = x | (x >> 2);
301 x = x | (x >> 4);
302 x = x | (x >> 8);
303 x = x | (x >>16);
304 return x + 1;
305}
306
307
308 /** \brief Round down to the nearest power of 2.
309
310 Efficient algorithm for finding the largest power of 2 which is not greater than \a x
311 (function flp2() from Henry Warren: "Hacker's Delight", Addison-Wesley, 2003,
312 see http://www.hackersdelight.org/).
313
314 <b>\#include</b> <vigra/mathutil.hxx><br>
315 Namespace: vigra
316 */
318{
319 x = x | (x >> 1);
320 x = x | (x >> 2);
321 x = x | (x >> 4);
322 x = x | (x >> 8);
323 x = x | (x >>16);
324 return x - (x >> 1);
325}
326
327namespace detail {
328
329template <class T>
330struct IntLog2
331{
332 static Int32 table[64];
333};
334
335template <class T>
336Int32 IntLog2<T>::table[64] = {
337 -1, 0, -1, 15, -1, 1, 28, -1, 16, -1, -1, -1, 2, 21,
338 29, -1, -1, -1, 19, 17, 10, -1, 12, -1, -1, 3, -1, 6,
339 -1, 22, 30, -1, 14, -1, 27, -1, -1, -1, 20, -1, 18, 9,
340 11, -1, 5, -1, -1, 13, 26, -1, -1, 8, -1, 4, -1, 25,
341 -1, 7, 24, -1, 23, -1, 31, -1};
342
343} // namespace detail
344
345 /** \brief Compute the base-2 logarithm of an integer.
346
347 Returns the position of the left-most 1-bit in the given number \a x, or
348 -1 if \a x == 0. That is,
349
350 \code
351 assert(k >= 0 && k < 32 && log2i(1 << k) == k);
352 \endcode
353
354 The function uses Robert Harley's algorithm to determine the number of leading zeros
355 in \a x (algorithm nlz10() at http://www.hackersdelight.org/). But note that the functions
356 \ref floorPower2() or \ref ceilPower2() are more efficient and should be preferred when possible.
357
358 <b>\#include</b> <vigra/mathutil.hxx><br>
359 Namespace: vigra
360 */
362{
363 // Propagate leftmost 1-bit to the right.
364 x = x | (x >> 1);
365 x = x | (x >> 2);
366 x = x | (x >> 4);
367 x = x | (x >> 8);
368 x = x | (x >>16);
369 x = x*0x06EB14F9; // Multiplier is 7*255**3.
370 return detail::IntLog2<Int32>::table[x >> 26];
371}
372
373 /** \brief The square function.
374
375 <tt>sq(x) = x*x</tt> is needed so often that it makes sense to define it as a function.
376
377 <b>\#include</b> <vigra/mathutil.hxx><br>
378 Namespace: vigra
379 */
380template <class T>
381inline
382typename NumericTraits<T>::Promote sq(T t)
383{
384 return t*t;
385}
386
387namespace detail {
388
389template <class V, unsigned>
390struct cond_mult
391{
392 static V call(const V & x, const V & y) { return x * y; }
393};
394template <class V>
395struct cond_mult<V, 0>
396{
397 static V call(const V &, const V & y) { return y; }
398};
399
400template <class V, unsigned n>
401struct power_static
402{
403 static V call(const V & x)
404 {
405 return n / 2
406 ? cond_mult<V, n & 1>::call(x, power_static<V, n / 2>::call(x * x))
407 : n & 1 ? x : V();
408 }
409};
410template <class V>
411struct power_static<V, 0>
412{
413 static V call(const V & /* x */)
414 {
415 return V(1);
416 }
417};
418
419} // namespace detail
420
421 /** \brief Exponentiation to a positive integer power by squaring.
422
423 <b>\#include</b> <vigra/mathutil.hxx><br>
424 Namespace: vigra
425 */
426template <unsigned n, class V>
427inline V power(const V & x)
428{
429 return detail::power_static<V, n>::call(x);
430}
431//doxygen_overloaded_function(template <unsigned n, class V> power(const V & x))
432
433namespace detail {
434
435template <class T>
436struct IntSquareRoot
437{
438 static UInt32 sqq_table[];
439 static UInt32 exec(UInt32 v);
440};
441
442template <class T>
443UInt32 IntSquareRoot<T>::sqq_table[] = {
444 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57,
445 59, 61, 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83,
446 84, 86, 87, 89, 90, 91, 93, 94, 96, 97, 98, 99, 101, 102,
447 103, 104, 106, 107, 108, 109, 110, 112, 113, 114, 115, 116, 117, 118,
448 119, 120, 121, 122, 123, 124, 125, 126, 128, 128, 129, 130, 131, 132,
449 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 144, 145,
450 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, 156, 157,
451 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168,
452 169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178,
453 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188,
454 189, 189, 190, 191, 192, 192, 193, 193, 194, 195, 195, 196, 197, 197,
455 198, 199, 199, 200, 201, 201, 202, 203, 203, 204, 204, 205, 206, 206,
456 207, 208, 208, 209, 209, 210, 211, 211, 212, 212, 213, 214, 214, 215,
457 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, 221, 222, 222, 223,
458 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, 230, 231,
459 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238,
460 239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246,
461 246, 247, 247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253,
462 253, 254, 254, 255
463};
464
465template <class T>
466UInt32 IntSquareRoot<T>::exec(UInt32 x)
467{
468 UInt32 xn;
469 if (x >= 0x10000)
470 if (x >= 0x1000000)
471 if (x >= 0x10000000)
472 if (x >= 0x40000000) {
473 if (x >= (UInt32)65535*(UInt32)65535)
474 return 65535;
475 xn = sqq_table[x>>24] << 8;
476 } else
477 xn = sqq_table[x>>22] << 7;
478 else
479 if (x >= 0x4000000)
480 xn = sqq_table[x>>20] << 6;
481 else
482 xn = sqq_table[x>>18] << 5;
483 else {
484 if (x >= 0x100000)
485 if (x >= 0x400000)
486 xn = sqq_table[x>>16] << 4;
487 else
488 xn = sqq_table[x>>14] << 3;
489 else
490 if (x >= 0x40000)
491 xn = sqq_table[x>>12] << 2;
492 else
493 xn = sqq_table[x>>10] << 1;
494
495 goto nr1;
496 }
497 else
498 if (x >= 0x100) {
499 if (x >= 0x1000)
500 if (x >= 0x4000)
501 xn = (sqq_table[x>>8] >> 0) + 1;
502 else
503 xn = (sqq_table[x>>6] >> 1) + 1;
504 else
505 if (x >= 0x400)
506 xn = (sqq_table[x>>4] >> 2) + 1;
507 else
508 xn = (sqq_table[x>>2] >> 3) + 1;
509
510 goto adj;
511 } else
512 return sqq_table[x] >> 4;
513
514 /* Run two iterations of the standard convergence formula */
515
516 xn = (xn + 1 + x / xn) / 2;
517 nr1:
518 xn = (xn + 1 + x / xn) / 2;
519 adj:
520
521 if (xn * xn > x) /* Correct rounding if necessary */
522 xn--;
523
524 return xn;
525}
526
527} // namespace detail
528
529using VIGRA_CSTD::sqrt;
530
531 /** \brief Signed integer square root.
532
533 Useful for fast fixed-point computations.
534
535 <b>\#include</b> <vigra/mathutil.hxx><br>
536 Namespace: vigra
537 */
539{
540 if(v < 0)
541 throw std::domain_error("sqrti(Int32): negative argument.");
542 return (Int32)detail::IntSquareRoot<UInt32>::exec((UInt32)v);
543}
544
545 /** \brief Unsigned integer square root.
546
547 Useful for fast fixed-point computations.
548
549 <b>\#include</b> <vigra/mathutil.hxx><br>
550 Namespace: vigra
551 */
552inline UInt32 sqrti(UInt32 v)
553{
554 return detail::IntSquareRoot<UInt32>::exec(v);
555}
556
557#ifdef VIGRA_NO_HYPOT
558 /** \brief Compute the Euclidean distance (length of the hypotenuse of a right-angled triangle).
559
560 The hypot() function returns the sqrt(a*a + b*b).
561 It is implemented in a way that minimizes round-off error.
562
563 <b>\#include</b> <vigra/mathutil.hxx><br>
564 Namespace: vigra
565 */
566inline double hypot(double a, double b)
567{
568 double absa = VIGRA_CSTD::fabs(a), absb = VIGRA_CSTD::fabs(b);
569 if (absa > absb)
570 return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa));
571 else
572 return absb == 0.0
573 ? 0.0
574 : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb));
575}
576
577#else
578
579using ::hypot;
580
581#endif
582
583 /** \brief The sign function.
584
585 Returns 1, 0, or -1 depending on the sign of \a t, but with the same type as \a t.
586
587 <b>\#include</b> <vigra/mathutil.hxx><br>
588 Namespace: vigra
589 */
590template <class T>
591inline T sign(T t)
592{
593 return t > NumericTraits<T>::zero()
594 ? NumericTraits<T>::one()
596 ? -NumericTraits<T>::one()
597 : NumericTraits<T>::zero();
598}
599
600 /** \brief The integer sign function.
601
602 Returns 1, 0, or -1 depending on the sign of \a t, converted to int.
603
604 <b>\#include</b> <vigra/mathutil.hxx><br>
605 Namespace: vigra
606 */
607template <class T>
608inline int signi(T t)
609{
610 return t > NumericTraits<T>::zero()
611 ? 1
613 ? -1
614 : 0;
615}
616
617 /** \brief The binary sign function.
618
619 Transfers the sign of \a t2 to \a t1.
620
621 <b>\#include</b> <vigra/mathutil.hxx><br>
622 Namespace: vigra
623 */
624template <class T1, class T2>
625inline T1 sign(T1 t1, T2 t2)
626{
627 return t2 >= NumericTraits<T2>::zero()
628 ? abs(t1)
629 : -abs(t1);
630}
631
632
633#ifdef DOXYGEN // only for documentation
634 /** \brief Check if an integer is even.
635
636 Defined for all integral types.
637 */
638bool even(int t);
639
640 /** \brief Check if an integer is odd.
641
642 Defined for all integral types.
643 */
644bool odd(int t);
645
646#endif
647
648#define VIGRA_DEFINE_ODD_EVEN(T) \
649 inline bool even(T t) { return (t&1) == 0; } \
650 inline bool odd(T t) { return (t&1) == 1; }
651
652VIGRA_DEFINE_ODD_EVEN(char)
653VIGRA_DEFINE_ODD_EVEN(short)
654VIGRA_DEFINE_ODD_EVEN(int)
655VIGRA_DEFINE_ODD_EVEN(long)
656VIGRA_DEFINE_ODD_EVEN(long long)
657VIGRA_DEFINE_ODD_EVEN(unsigned char)
658VIGRA_DEFINE_ODD_EVEN(unsigned short)
659VIGRA_DEFINE_ODD_EVEN(unsigned int)
660VIGRA_DEFINE_ODD_EVEN(unsigned long)
661VIGRA_DEFINE_ODD_EVEN(unsigned long long)
662
663#undef VIGRA_DEFINE_ODD_EVEN
664
665#define VIGRA_DEFINE_NORM(T) \
666 inline NormTraits<T>::SquaredNormType squaredNorm(T t) { return sq(t); } \
667 inline NormTraits<T>::NormType norm(T t) { return abs(t); }
668
669VIGRA_DEFINE_NORM(bool)
670VIGRA_DEFINE_NORM(signed char)
671VIGRA_DEFINE_NORM(unsigned char)
672VIGRA_DEFINE_NORM(short)
673VIGRA_DEFINE_NORM(unsigned short)
674VIGRA_DEFINE_NORM(int)
675VIGRA_DEFINE_NORM(unsigned int)
676VIGRA_DEFINE_NORM(long)
677VIGRA_DEFINE_NORM(unsigned long)
678VIGRA_DEFINE_NORM(long long)
679VIGRA_DEFINE_NORM(unsigned long long)
680VIGRA_DEFINE_NORM(float)
681VIGRA_DEFINE_NORM(double)
682VIGRA_DEFINE_NORM(long double)
683
684#undef VIGRA_DEFINE_NORM
685
686template <class T>
687inline typename NormTraits<std::complex<T> >::SquaredNormType
688squaredNorm(std::complex<T> const & t)
689{
690 return sq(t.real()) + sq(t.imag());
691}
692
693#ifdef DOXYGEN // only for documentation
694 /** \brief The squared norm of a numerical object.
695
696 <ul>
697 <li>For scalar types: equals <tt>vigra::sq(t)</tt>.
698 <li>For vectorial types (including TinyVector): equals <tt>vigra::dot(t, t)</tt>.
699 <li>For complex number types: equals <tt>vigra::sq(t.real()) + vigra::sq(t.imag())</tt>.
700 <li>For array and matrix types: results in the squared Frobenius norm (sum of squares of the matrix elements).
701 </ul>
702 */
703NormTraits<T>::SquaredNormType squaredNorm(T const & t);
704
705#endif
706
707 /** \brief The norm of a numerical object.
708
709 For scalar types: implemented as <tt>abs(t)</tt><br>
710 otherwise: implemented as <tt>sqrt(squaredNorm(t))</tt>.
711
712 <b>\#include</b> <vigra/mathutil.hxx><br>
713 Namespace: vigra
714 */
715template <class T>
716inline typename NormTraits<T>::NormType
717norm(T const & t)
718{
719 typedef typename NormTraits<T>::SquaredNormType SNT;
720 return sqrt(static_cast<typename SquareRootTraits<SNT>::SquareRootArgument>(squaredNorm(t)));
721}
722
723 /** \brief Compute the eigenvalues of a 2x2 real symmetric matrix.
724
725 This uses the analytical eigenvalue formula
726 \f[
727 \lambda_{1,2} = \frac{1}{2}\left(a_{00} + a_{11} \pm \sqrt{(a_{00} - a_{11})^2 + 4 a_{01}^2}\right)
728 \f]
729
730 <b>\#include</b> <vigra/mathutil.hxx><br>
731 Namespace: vigra
732 */
733template <class T>
734void symmetric2x2Eigenvalues(T a00, T a01, T a11, T * r0, T * r1)
735{
736 double d = hypot(a00 - a11, 2.0*a01);
737 *r0 = static_cast<T>(0.5*(a00 + a11 + d));
738 *r1 = static_cast<T>(0.5*(a00 + a11 - d));
739 if(*r0 < *r1)
740 std::swap(*r0, *r1);
741}
742
743 /** \brief Compute the eigenvalues of a 3x3 real symmetric matrix.
744
745 This uses a numerically stable version of the analytical eigenvalue formula according to
746 <p>
747 David Eberly: <a href="http://www.geometrictools.com/Documentation/EigenSymmetric3x3.pdf">
748 <em>"Eigensystems for 3 x 3 Symmetric Matrices (Revisited)"</em></a>, Geometric Tools Documentation, 2006
749
750 <b>\#include</b> <vigra/mathutil.hxx><br>
751 Namespace: vigra
752 */
753template <class T>
755 T * r0, T * r1, T * r2)
756{
757 double inv3 = 1.0 / 3.0, root3 = std::sqrt(3.0);
758
759 double c0 = a00*a11*a22 + 2.0*a01*a02*a12 - a00*a12*a12 - a11*a02*a02 - a22*a01*a01;
760 double c1 = a00*a11 - a01*a01 + a00*a22 - a02*a02 + a11*a22 - a12*a12;
761 double c2 = a00 + a11 + a22;
762 double c2Div3 = c2*inv3;
763 double aDiv3 = (c1 - c2*c2Div3)*inv3;
764 if (aDiv3 > 0.0)
765 aDiv3 = 0.0;
766 double mbDiv2 = 0.5*(c0 + c2Div3*(2.0*c2Div3*c2Div3 - c1));
767 double q = mbDiv2*mbDiv2 + aDiv3*aDiv3*aDiv3;
768 if (q > 0.0)
769 q = 0.0;
770 double magnitude = std::sqrt(-aDiv3);
771 double angle = std::atan2(std::sqrt(-q),mbDiv2)*inv3;
772 double cs = std::cos(angle);
773 double sn = std::sin(angle);
774 *r0 = static_cast<T>(c2Div3 + 2.0*magnitude*cs);
775 *r1 = static_cast<T>(c2Div3 - magnitude*(cs + root3*sn));
776 *r2 = static_cast<T>(c2Div3 - magnitude*(cs - root3*sn));
777 if(*r0 < *r1)
778 std::swap(*r0, *r1);
779 if(*r0 < *r2)
780 std::swap(*r0, *r2);
781 if(*r1 < *r2)
782 std::swap(*r1, *r2);
783}
784
785namespace detail {
786
787template <class T>
788T ellipticRD(T x, T y, T z)
789{
790 double f = 1.0, s = 0.0, X, Y, Z, m;
791 for(;;)
792 {
793 m = (x + y + 3.0*z) / 5.0;
794 X = 1.0 - x/m;
795 Y = 1.0 - y/m;
796 Z = 1.0 - z/m;
797 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
798 break;
799 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
800 s += f / (VIGRA_CSTD::sqrt(z)*(z + l));
801 f /= 4.0;
802 x = (x + l)/4.0;
803 y = (y + l)/4.0;
804 z = (z + l)/4.0;
805 }
806 double a = X*Y;
807 double b = sq(Z);
808 double c = a - b;
809 double d = a - 6.0*b;
810 double e = d + 2.0*c;
811 return 3.0*s + f*(1.0+d*(-3.0/14.0+d*9.0/88.0-Z*e*4.5/26.0)
812 +Z*(e/6.0+Z*(-c*9.0/22.0+a*Z*3.0/26.0))) / VIGRA_CSTD::pow(m,1.5);
813}
814
815template <class T>
816T ellipticRF(T x, T y, T z)
817{
818 double X, Y, Z, m;
819 for(;;)
820 {
821 m = (x + y + z) / 3.0;
822 X = 1.0 - x/m;
823 Y = 1.0 - y/m;
824 Z = 1.0 - z/m;
825 if(std::max(std::max(VIGRA_CSTD::fabs(X), VIGRA_CSTD::fabs(Y)), VIGRA_CSTD::fabs(Z)) < 0.01)
826 break;
827 double l = VIGRA_CSTD::sqrt(x*y) + VIGRA_CSTD::sqrt(x*z) + VIGRA_CSTD::sqrt(y*z);
828 x = (x + l)/4.0;
829 y = (y + l)/4.0;
830 z = (z + l)/4.0;
831 }
832 double d = X*Y - sq(Z);
833 double p = X*Y*Z;
834 return (1.0 - d/10.0 + p/14.0 + sq(d)/24.0 - d*p*3.0/44.0) / VIGRA_CSTD::sqrt(m);
835}
836
837} // namespace detail
838
839 /** \brief The incomplete elliptic integral of the first kind.
840
841 This function computes
842
843 \f[
844 \mbox{F}(x, k) = \int_0^x \frac{1}{\sqrt{1 - k^2 \sin(t)^2}} dt
845 \f]
846
847 according to the algorithm given in Press et al. "Numerical Recipes".
848
849 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral
850 functions must be k^2 rather than k. Check the documentation when results disagree!
851
852 <b>\#include</b> <vigra/mathutil.hxx><br>
853 Namespace: vigra
854 */
855inline double ellipticIntegralF(double x, double k)
856{
857 double c2 = sq(VIGRA_CSTD::cos(x));
858 double s = VIGRA_CSTD::sin(x);
859 return s*detail::ellipticRF(c2, 1.0 - sq(k*s), 1.0);
860}
861
862 /** \brief The incomplete elliptic integral of the second kind.
863
864 This function computes
865
866 \f[
867 \mbox{E}(x, k) = \int_0^x \sqrt{1 - k^2 \sin(t)^2} dt
868 \f]
869
870 according to the algorithm given in Press et al. "Numerical Recipes". The
871 complete elliptic integral of the second kind is simply <tt>ellipticIntegralE(M_PI/2, k)</TT>.
872
873 Note: In some libraries (e.g. Mathematica), the second parameter of the elliptic integral
874 functions must be k^2 rather than k. Check the documentation when results disagree!
875
876 <b>\#include</b> <vigra/mathutil.hxx><br>
877 Namespace: vigra
878 */
879inline double ellipticIntegralE(double x, double k)
880{
881 double c2 = sq(VIGRA_CSTD::cos(x));
882 double s = VIGRA_CSTD::sin(x);
883 k = sq(k*s);
884 return s*(detail::ellipticRF(c2, 1.0-k, 1.0) - k/3.0*detail::ellipticRD(c2, 1.0-k, 1.0));
885}
886
887#if defined(_MSC_VER) && _MSC_VER < 1800
888
889namespace detail {
890
891template <class T>
892double erfImpl(T x)
893{
894 double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
895 double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
896 t*(0.09678418+t*(-0.18628806+t*(0.27886807+
897 t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
898 t*0.17087277)))))))));
899 if (x >= 0.0)
900 return 1.0 - ans;
901 else
902 return ans - 1.0;
903}
904
905} // namespace detail
906
907 /** \brief The error function.
908
909 If <tt>erf()</tt> is not provided in the C standard math library (as it should according to the
910 new C99 standard ?), VIGRA implements <tt>erf()</tt> as an approximation of the error
911 function
912
913 \f[
914 \mbox{erf}(x) = \int_0^x e^{-t^2} dt
915 \f]
916
917 according to the formula given in Press et al. "Numerical Recipes".
918
919 <b>\#include</b> <vigra/mathutil.hxx><br>
920 Namespace: vigra
921 */
922inline double erf(double x)
923{
924 return detail::erfImpl(x);
925}
926
927#else
928
929using ::erf;
930
931#endif
932
933namespace detail {
934
935template <class T>
936double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, T noncentrality, T arg)
937{
938 double a = degreesOfFreedom + noncentrality;
939 double b = (a + noncentrality) / sq(a);
940 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);
941 return 0.5*(1.0 + erf(t/VIGRA_CSTD::sqrt(2.0)));
942}
943
944template <class T>
945void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j)
946{
947 double tol = -50.0;
948 if(lans < tol)
949 {
950 lans = lans + VIGRA_CSTD::log(arg / j);
951 dans = VIGRA_CSTD::exp(lans);
952 }
953 else
954 {
955 dans = dans * arg / j;
956 }
957 pans = pans - dans;
958 j += 2;
959}
960
961template <class T>
962std::pair<double, double> noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
963{
964 vigra_precondition(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,
965 "noncentralChi2P(): parameters must be positive.");
966 if (arg == 0.0 && degreesOfFreedom > 0)
967 return std::make_pair(0.0, 0.0);
968
969 // Determine initial values
970 double b1 = 0.5 * noncentrality,
971 ao = VIGRA_CSTD::exp(-b1),
972 eps2 = eps / ao,
973 lnrtpi2 = 0.22579135264473,
974 probability, density, lans, dans, pans, sum, am, hold;
975 unsigned int maxit = 500,
976 i, m;
977 if(degreesOfFreedom % 2)
978 {
979 i = 1;
980 lans = -0.5 * (arg + VIGRA_CSTD::log(arg)) - lnrtpi2;
981 dans = VIGRA_CSTD::exp(lans);
982 pans = erf(VIGRA_CSTD::sqrt(arg/2.0));
983 }
984 else
985 {
986 i = 2;
987 lans = -0.5 * arg;
988 dans = VIGRA_CSTD::exp(lans);
989 pans = 1.0 - dans;
990 }
991
992 // Evaluate first term
993 if(degreesOfFreedom == 0)
994 {
995 m = 1;
996 degreesOfFreedom = 2;
997 am = b1;
998 sum = 1.0 / ao - 1.0 - am;
999 density = am * dans;
1000 probability = 1.0 + am * pans;
1001 }
1002 else
1003 {
1004 m = 0;
1005 degreesOfFreedom = degreesOfFreedom - 1;
1006 am = 1.0;
1007 sum = 1.0 / ao - 1.0;
1008 while(i < degreesOfFreedom)
1009 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
1010 degreesOfFreedom = degreesOfFreedom + 1;
1011 density = dans;
1012 probability = pans;
1013 }
1014 // Evaluate successive terms of the expansion
1015 for(++m; m<maxit; ++m)
1016 {
1017 am = b1 * am / m;
1018 detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
1019 sum = sum - am;
1020 density = density + am * dans;
1021 hold = am * pans;
1022 probability = probability + hold;
1023 if((pans * sum < eps2) && (hold < eps2))
1024 break; // converged
1025 }
1026 if(m == maxit)
1027 vigra_fail("noncentralChi2P(): no convergence.");
1028 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
1029}
1030
1031} // namespace detail
1032
1033 /** \brief Chi square distribution.
1034
1035 Computes the density of a chi square distribution with \a degreesOfFreedom
1036 and tolerance \a accuracy at the given argument \a arg
1037 by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
1038
1039 <b>\#include</b> <vigra/mathutil.hxx><br>
1040 Namespace: vigra
1041 */
1042inline double chi2(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
1043{
1044 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).first;
1045}
1046
1047 /** \brief Cumulative chi square distribution.
1048
1049 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom
1050 and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
1051 a random number drawn from the distribution is below \a arg
1052 by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
1053
1054 <b>\#include</b> <vigra/mathutil.hxx><br>
1055 Namespace: vigra
1056 */
1057inline double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
1058{
1059 return detail::noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy).second;
1060}
1061
1062 /** \brief Non-central chi square distribution.
1063
1064 Computes the density of a non-central chi square distribution with \a degreesOfFreedom,
1065 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument
1066 \a arg. It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from
1067 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of
1068 degrees of freedom.
1069
1070 <b>\#include</b> <vigra/mathutil.hxx><br>
1071 Namespace: vigra
1072 */
1073inline double noncentralChi2(unsigned int degreesOfFreedom,
1074 double noncentrality, double arg, double accuracy = 1e-7)
1075{
1076 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).first;
1077}
1078
1079 /** \brief Cumulative non-central chi square distribution.
1080
1081 Computes the cumulative density of a chi square distribution with \a degreesOfFreedom,
1082 noncentrality parameter \a noncentrality and tolerance \a accuracy at the given argument
1083 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
1084 It uses Algorithm AS 231 from Appl. Statist. (1987) Vol.36, No.3 (code ported from
1085 http://lib.stat.cmu.edu/apstat/231). The algorithm has linear complexity in the number of
1086 degrees of freedom (see noncentralChi2CDFApprox() for a constant-time algorithm).
1087
1088 <b>\#include</b> <vigra/mathutil.hxx><br>
1089 Namespace: vigra
1090 */
1091inline double noncentralChi2CDF(unsigned int degreesOfFreedom,
1092 double noncentrality, double arg, double accuracy = 1e-7)
1093{
1094 return detail::noncentralChi2CDF(degreesOfFreedom, noncentrality, arg, accuracy).second;
1095}
1096
1097 /** \brief Cumulative non-central chi square distribution (approximate).
1098
1099 Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom,
1100 and noncentrality parameter \a noncentrality at the given argument
1101 \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
1102 It uses the approximate transform into a normal distribution due to Wilson and Hilferty
1103 (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32).
1104 The algorithm's running time is independent of the inputs, i.e. is should be used
1105 when noncentralChi2CDF() is too slow, and approximate values are sufficient. The accuracy is only
1106 about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
1107
1108 <b>\#include</b> <vigra/mathutil.hxx><br>
1109 Namespace: vigra
1110 */
1111inline double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg)
1112{
1113 return detail::noncentralChi2CDFApprox(degreesOfFreedom, noncentrality, arg);
1114}
1115
1116namespace detail {
1117
1118// computes (l+m)! / (l-m)!
1119// l and m must be positive
1120template <class T>
1121T facLM(T l, T m)
1122{
1123 T tmp = NumericTraits<T>::one();
1124 for(T f = l-m+1; f <= l+m; ++f)
1125 tmp *= f;
1126 return tmp;
1127}
1128
1129} // namespace detail
1130
1131 /** \brief Associated Legendre polynomial.
1132
1133 Computes the value of the associated Legendre polynomial of order <tt>l, m</tt>
1134 for argument <tt>x</tt>. <tt>x</tt> must be in the range <tt>[-1.0, 1.0]</tt>,
1135 otherwise an exception is thrown. The standard Legendre polynomials are the
1136 special case <tt>m == 0</tt>.
1137
1138 <b>\#include</b> <vigra/mathutil.hxx><br>
1139 Namespace: vigra
1140 */
1141template <class REAL>
1142REAL legendre(unsigned int l, int m, REAL x)
1143{
1144 vigra_precondition(abs(x) <= 1.0, "legendre(): x must be in [-1.0, 1.0].");
1145 if (m < 0)
1146 {
1147 m = -m;
1148 REAL s = odd(m)
1149 ? -1.0
1150 : 1.0;
1151 return legendre(l,m,x) * s / detail::facLM<REAL>(l,m);
1152 }
1153 REAL result = 1.0;
1154 if (m > 0)
1155 {
1156 REAL r = std::sqrt( (1.0-x) * (1.0+x) );
1157 REAL f = 1.0;
1158 for (int i=1; i<=m; i++)
1159 {
1160 result *= (-f) * r;
1161 f += 2.0;
1162 }
1163 }
1164 if((int)l == m)
1165 return result;
1166
1167 REAL result_1 = x * (2.0 * m + 1.0) * result;
1168 if((int)l == m+1)
1169 return result_1;
1170 REAL other = 0.0;
1171 for(unsigned int i = m+2; i <= l; ++i)
1172 {
1173 other = ( (2.0*i-1.0) * x * result_1 - (i+m-1.0)*result) / (i-m);
1174 result = result_1;
1175 result_1 = other;
1176 }
1177 return other;
1178}
1179
1180 /** \brief \brief Legendre polynomial.
1181
1182 Computes the value of the Legendre polynomial of order <tt>l</tt> for argument <tt>x</tt>.
1183 <tt>x</tt> must be in the range <tt>[-1.0, 1.0]</tt>, otherwise an exception is thrown.
1184
1185 <b>\#include</b> <vigra/mathutil.hxx><br>
1186 Namespace: vigra
1187 */
1188template <class REAL>
1189REAL legendre(unsigned int l, REAL x)
1190{
1191 return legendre(l, 0, x);
1192}
1193
1194 /** \brief sin(pi*x).
1195
1196 Essentially calls <tt>std::sin(M_PI*x)</tt> but uses a more accurate implementation
1197 to make sure that <tt>sin_pi(1.0) == 0.0</tt> (which does not hold for
1198 <tt>std::sin(M_PI)</tt> due to round-off error), and <tt>sin_pi(0.5) == 1.0</tt>.
1199
1200 <b>\#include</b> <vigra/mathutil.hxx><br>
1201 Namespace: vigra
1202 */
1203template <class REAL>
1204REAL sin_pi(REAL x)
1205{
1206 if(x < 0.0)
1207 return -sin_pi(-x);
1208 if(x < 0.5)
1209 return std::sin(M_PI * x);
1210
1211 bool invert = false;
1212 if(x < 1.0)
1213 {
1214 invert = true;
1215 x = -x;
1216 }
1217
1218 REAL rem = std::floor(x);
1219 if(odd((int)rem))
1220 invert = !invert;
1221 rem = x - rem;
1222 if(rem > 0.5)
1223 rem = 1.0 - rem;
1224 if(rem == 0.5)
1225 rem = NumericTraits<REAL>::one();
1226 else
1227 rem = std::sin(M_PI * rem);
1228 return invert
1229 ? -rem
1230 : rem;
1231}
1232
1233 /** \brief cos(pi*x).
1234
1235 Essentially calls <tt>std::cos(M_PI*x)</tt> but uses a more accurate implementation
1236 to make sure that <tt>cos_pi(1.0) == -1.0</tt> and <tt>cos_pi(0.5) == 0.0</tt>.
1237
1238 <b>\#include</b> <vigra/mathutil.hxx><br>
1239 Namespace: vigra
1240 */
1241template <class REAL>
1242REAL cos_pi(REAL x)
1243{
1244 return sin_pi(x+0.5);
1245}
1246
1247namespace detail {
1248
1249template <class REAL>
1250struct GammaImpl
1251{
1252 static REAL gamma(REAL x);
1253 static REAL loggamma(REAL x);
1254
1255 static double g[];
1256 static double a[];
1257 static double t[];
1258 static double u[];
1259 static double v[];
1260 static double s[];
1261 static double r[];
1262 static double w[];
1263};
1264
1265template <class REAL>
1266double GammaImpl<REAL>::g[] = {
1267 1.0,
1268 0.5772156649015329,
1269 -0.6558780715202538,
1270 -0.420026350340952e-1,
1271 0.1665386113822915,
1272 -0.421977345555443e-1,
1273 -0.9621971527877e-2,
1274 0.7218943246663e-2,
1275 -0.11651675918591e-2,
1276 -0.2152416741149e-3,
1277 0.1280502823882e-3,
1278 -0.201348547807e-4,
1279 -0.12504934821e-5,
1280 0.1133027232e-5,
1281 -0.2056338417e-6,
1282 0.6116095e-8,
1283 0.50020075e-8,
1284 -0.11812746e-8,
1285 0.1043427e-9,
1286 0.77823e-11,
1287 -0.36968e-11,
1288 0.51e-12,
1289 -0.206e-13,
1290 -0.54e-14,
1291 0.14e-14
1292};
1293
1294template <class REAL>
1295double GammaImpl<REAL>::a[] = {
1296 7.72156649015328655494e-02,
1297 3.22467033424113591611e-01,
1298 6.73523010531292681824e-02,
1299 2.05808084325167332806e-02,
1300 7.38555086081402883957e-03,
1301 2.89051383673415629091e-03,
1302 1.19270763183362067845e-03,
1303 5.10069792153511336608e-04,
1304 2.20862790713908385557e-04,
1305 1.08011567247583939954e-04,
1306 2.52144565451257326939e-05,
1307 4.48640949618915160150e-05
1308};
1309
1310template <class REAL>
1311double GammaImpl<REAL>::t[] = {
1312 4.83836122723810047042e-01,
1313 -1.47587722994593911752e-01,
1314 6.46249402391333854778e-02,
1315 -3.27885410759859649565e-02,
1316 1.79706750811820387126e-02,
1317 -1.03142241298341437450e-02,
1318 6.10053870246291332635e-03,
1319 -3.68452016781138256760e-03,
1320 2.25964780900612472250e-03,
1321 -1.40346469989232843813e-03,
1322 8.81081882437654011382e-04,
1323 -5.38595305356740546715e-04,
1324 3.15632070903625950361e-04,
1325 -3.12754168375120860518e-04,
1326 3.35529192635519073543e-04
1327};
1328
1329template <class REAL>
1330double GammaImpl<REAL>::u[] = {
1331 -7.72156649015328655494e-02,
1332 6.32827064025093366517e-01,
1333 1.45492250137234768737e+00,
1334 9.77717527963372745603e-01,
1335 2.28963728064692451092e-01,
1336 1.33810918536787660377e-02
1337};
1338
1339template <class REAL>
1340double GammaImpl<REAL>::v[] = {
1341 0.0,
1342 2.45597793713041134822e+00,
1343 2.12848976379893395361e+00,
1344 7.69285150456672783825e-01,
1345 1.04222645593369134254e-01,
1346 3.21709242282423911810e-03
1347};
1348
1349template <class REAL>
1350double GammaImpl<REAL>::s[] = {
1351 -7.72156649015328655494e-02,
1352 2.14982415960608852501e-01,
1353 3.25778796408930981787e-01,
1354 1.46350472652464452805e-01,
1355 2.66422703033638609560e-02,
1356 1.84028451407337715652e-03,
1357 3.19475326584100867617e-05
1358};
1359
1360template <class REAL>
1361double GammaImpl<REAL>::r[] = {
1362 0.0,
1363 1.39200533467621045958e+00,
1364 7.21935547567138069525e-01,
1365 1.71933865632803078993e-01,
1366 1.86459191715652901344e-02,
1367 7.77942496381893596434e-04,
1368 7.32668430744625636189e-06
1369};
1370
1371template <class REAL>
1372double GammaImpl<REAL>::w[] = {
1373 4.18938533204672725052e-01,
1374 8.33333333333329678849e-02,
1375 -2.77777777728775536470e-03,
1376 7.93650558643019558500e-04,
1377 -5.95187557450339963135e-04,
1378 8.36339918996282139126e-04,
1379 -1.63092934096575273989e-03
1380};
1381
1382template <class REAL>
1383REAL GammaImpl<REAL>::gamma(REAL x)
1384{
1385 int i, k, m, ix = (int)x;
1386 double ga = 0.0, gr = 0.0, r = 0.0, z = 0.0;
1387
1388 vigra_precondition(x <= 171.0,
1389 "gamma(): argument cannot exceed 171.0.");
1390
1391 if (x == ix)
1392 {
1393 if (ix > 0)
1394 {
1395 ga = 1.0; // use factorial
1396 for (i=2; i<ix; ++i)
1397 {
1398 ga *= i;
1399 }
1400 }
1401 else
1402 {
1403 vigra_precondition(false,
1404 "gamma(): gamma function is undefined for 0 and negative integers.");
1405 }
1406 }
1407 else
1408 {
1409 if (abs(x) > 1.0)
1410 {
1411 z = abs(x);
1412 m = (int)z;
1413 r = 1.0;
1414 for (k=1; k<=m; ++k)
1415 {
1416 r *= (z-k);
1417 }
1418 z -= m;
1419 }
1420 else
1421 {
1422 z = x;
1423 }
1424 gr = g[24];
1425 for (k=23; k>=0; --k)
1426 {
1427 gr = gr*z+g[k];
1428 }
1429 ga = 1.0/(gr*z);
1430 if (abs(x) > 1.0)
1431 {
1432 ga *= r;
1433 if (x < 0.0)
1434 {
1435 ga = -M_PI/(x*ga*sin_pi(x));
1436 }
1437 }
1438 }
1439 return ga;
1440}
1441
1442/*
1443 * the following code is derived from lgamma_r() by Sun
1444 *
1445 * ====================================================
1446 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
1447 *
1448 * Developed at SunPro, a Sun Microsystems, Inc. business.
1449 * Permission to use, copy, modify, and distribute this
1450 * software is freely granted, provided that this notice
1451 * is preserved.
1452 * ====================================================
1453 *
1454 */
1455template <class REAL>
1456REAL GammaImpl<REAL>::loggamma(REAL x)
1457{
1458 vigra_precondition(x > 0.0,
1459 "loggamma(): argument must be positive.");
1460
1461 vigra_precondition(x <= 1.0e307,
1462 "loggamma(): argument must not exceed 1e307.");
1463
1464 double res;
1465
1466 if (x < 4.2351647362715017e-22)
1467 {
1468 res = -std::log(x);
1469 }
1470 else if ((x == 2.0) || (x == 1.0))
1471 {
1472 res = 0.0;
1473 }
1474 else if (x < 2.0)
1475 {
1476 const double tc = 1.46163214496836224576e+00;
1477 const double tf = -1.21486290535849611461e-01;
1478 const double tt = -3.63867699703950536541e-18;
1479 if (x <= 0.9)
1480 {
1481 res = -std::log(x);
1482 if (x >= 0.7316)
1483 {
1484 double y = 1.0-x;
1485 double z = y*y;
1486 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1487 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1488 double p = y*p1+p2;
1489 res += (p-0.5*y);
1490 }
1491 else if (x >= 0.23164)
1492 {
1493 double y = x-(tc-1.0);
1494 double z = y*y;
1495 double w = z*y;
1496 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1497 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1498 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1499 double p = z*p1-(tt-w*(p2+y*p3));
1500 res += (tf + p);
1501 }
1502 else
1503 {
1504 double y = x;
1505 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1506 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1507 res += (-0.5*y + p1/p2);
1508 }
1509 }
1510 else
1511 {
1512 res = 0.0;
1513 if (x >= 1.7316)
1514 {
1515 double y = 2.0-x;
1516 double z = y*y;
1517 double p1 = a[0]+z*(a[2]+z*(a[4]+z*(a[6]+z*(a[8]+z*a[10]))));
1518 double p2 = z*(a[1]+z*(a[3]+z*(a[5]+z*(a[7]+z*(a[9]+z*a[11])))));
1519 double p = y*p1+p2;
1520 res += (p-0.5*y);
1521 }
1522 else if(x >= 1.23164)
1523 {
1524 double y = x-tc;
1525 double z = y*y;
1526 double w = z*y;
1527 double p1 = t[0]+w*(t[3]+w*(t[6]+w*(t[9] +w*t[12])));
1528 double p2 = t[1]+w*(t[4]+w*(t[7]+w*(t[10]+w*t[13])));
1529 double p3 = t[2]+w*(t[5]+w*(t[8]+w*(t[11]+w*t[14])));
1530 double p = z*p1-(tt-w*(p2+y*p3));
1531 res += (tf + p);
1532 }
1533 else
1534 {
1535 double y = x-1.0;
1536 double p1 = y*(u[0]+y*(u[1]+y*(u[2]+y*(u[3]+y*(u[4]+y*u[5])))));
1537 double p2 = 1.0+y*(v[1]+y*(v[2]+y*(v[3]+y*(v[4]+y*v[5]))));
1538 res += (-0.5*y + p1/p2);
1539 }
1540 }
1541 }
1542 else if(x < 8.0)
1543 {
1544 double i = std::floor(x);
1545 double y = x-i;
1546 double p = y*(s[0]+y*(s[1]+y*(s[2]+y*(s[3]+y*(s[4]+y*(s[5]+y*s[6]))))));
1547 double q = 1.0+y*(r[1]+y*(r[2]+y*(r[3]+y*(r[4]+y*(r[5]+y*r[6])))));
1548 res = 0.5*y+p/q;
1549 double z = 1.0;
1550 while (i > 2.0)
1551 {
1552 --i;
1553 z *= (y+i);
1554 }
1555 res += std::log(z);
1556 }
1557 else if (x < 2.8823037615171174e+17)
1558 {
1559 double t = std::log(x);
1560 double z = 1.0/x;
1561 double y = z*z;
1562 double yy = w[0]+z*(w[1]+y*(w[2]+y*(w[3]+y*(w[4]+y*(w[5]+y*w[6])))));
1563 res = (x-0.5)*(t-1.0)+yy;
1564 }
1565 else
1566 {
1567 res = x*(std::log(x) - 1.0);
1568 }
1569
1570 return res;
1571}
1572
1573
1574} // namespace detail
1575
1576 /** \brief The gamma function.
1577
1578 This function implements the algorithm from<br>
1579 Zhang and Jin: "Computation of Special Functions", John Wiley and Sons, 1996.
1580
1581 The argument must be <= 171.0 and cannot be zero or a negative integer. An
1582 exception is thrown when these conditions are violated.
1583
1584 <b>\#include</b> <vigra/mathutil.hxx><br>
1585 Namespace: vigra
1586 */
1587inline double gamma(double x)
1588{
1589 return detail::GammaImpl<double>::gamma(x);
1590}
1591
1592 /** \brief The natural logarithm of the gamma function.
1593
1594 This function is based on a free implementation by Sun Microsystems, Inc., see
1595 <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
1596 math functions.
1597
1598 The argument must be positive and < 1e30. An exception is thrown when these conditions are violated.
1599
1600 <b>\#include</b> <vigra/mathutil.hxx><br>
1601 Namespace: vigra
1602 */
1603inline double loggamma(double x)
1604{
1605 return detail::GammaImpl<double>::loggamma(x);
1606}
1607
1608
1609namespace detail {
1610
1611// both f1 and f2 are unsigned here
1612template<class FPT>
1613inline
1614FPT safeFloatDivision( FPT f1, FPT f2 )
1615{
1616 return f2 < NumericTraits<FPT>::one() && f1 > f2 * NumericTraits<FPT>::max()
1617 ? NumericTraits<FPT>::max()
1618 : (f2 > NumericTraits<FPT>::one() && f1 < f2 * NumericTraits<FPT>::smallestPositive()) ||
1619 f1 == NumericTraits<FPT>::zero()
1620 ? NumericTraits<FPT>::zero()
1621 : f1/f2;
1622}
1623
1624} // namespace detail
1625
1626 /** \brief Tolerance based floating-point equality.
1627
1628 Check whether two floating point numbers are equal within the given tolerance.
1629 This is useful because floating point numbers that should be equal in theory are
1630 rarely exactly equal in practice. If the tolerance \a epsilon is not given,
1631 twice the machine epsilon is used.
1632
1633 <b>\#include</b> <vigra/mathutil.hxx><br>
1634 Namespace: vigra
1635 */
1636template <class T1, class T2>
1637bool
1638closeAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon)
1639{
1640 typedef typename PromoteTraits<T1, T2>::Promote T;
1641 if(l == 0.0)
1642 return VIGRA_CSTD::fabs(r) <= epsilon;
1643 if(r == 0.0)
1644 return VIGRA_CSTD::fabs(l) <= epsilon;
1645 T diff = VIGRA_CSTD::fabs( l - r );
1646 T d1 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( r ) );
1647 T d2 = detail::safeFloatDivision<T>( diff, VIGRA_CSTD::fabs( l ) );
1648
1649 return (d1 <= epsilon && d2 <= epsilon);
1650}
1651
1652template <class T1, class T2>
1653inline bool closeAtTolerance(T1 l, T2 r)
1654{
1655 typedef typename PromoteTraits<T1, T2>::Promote T;
1656 return closeAtTolerance(l, r, T(2.0) * NumericTraits<T>::epsilon());
1657}
1658
1659 /** \brief Tolerance based floating-point less-or-equal.
1660
1661 Check whether two floating point numbers are less or equal within the given tolerance.
1662 That is, \a l can actually be greater than \a r within the given \a epsilon.
1663 This is useful because floating point numbers that should be equal in theory are
1664 rarely exactly equal in practice. If the tolerance \a epsilon is not given,
1665 twice the machine epsilon is used.
1666
1667 <b>\#include</b> <vigra/mathutil.hxx><br>
1668 Namespace: vigra
1669 */
1670template <class T1, class T2>
1671inline bool
1672lessEqualAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon)
1673{
1674 return l < r || closeAtTolerance(l, r, epsilon);
1675}
1676
1677template <class T1, class T2>
1678inline bool lessEqualAtTolerance(T1 l, T2 r)
1679{
1680 typedef typename PromoteTraits<T1, T2>::Promote T;
1681 return lessEqualAtTolerance(l, r, T(2.0) * NumericTraits<T>::epsilon());
1682}
1683
1684 /** \brief Tolerance based floating-point greater-or-equal.
1685
1686 Check whether two floating point numbers are greater or equal within the given tolerance.
1687 That is, \a l can actually be less than \a r within the given \a epsilon.
1688 This is useful because floating point numbers that should be equal in theory are
1689 rarely exactly equal in practice. If the tolerance \a epsilon is not given,
1690 twice the machine epsilon is used.
1691
1692 <b>\#include</b> <vigra/mathutil.hxx><br>
1693 Namespace: vigra
1694 */
1695template <class T1, class T2>
1696inline bool
1697greaterEqualAtTolerance(T1 l, T2 r, typename PromoteTraits<T1, T2>::Promote epsilon)
1698{
1699 return r < l || closeAtTolerance(l, r, epsilon);
1700}
1701
1702template <class T1, class T2>
1703inline bool greaterEqualAtTolerance(T1 l, T2 r)
1704{
1705 typedef typename PromoteTraits<T1, T2>::Promote T;
1706 return greaterEqualAtTolerance(l, r, T(2.0) * NumericTraits<T>::epsilon());
1707}
1708
1709//@}
1710
1711#define VIGRA_MATH_FUNC_HELPER(TYPE) \
1712 inline TYPE clipLower(const TYPE t){ \
1713 return t < static_cast<TYPE>(0.0) ? static_cast<TYPE>(0.0) : t; \
1714 } \
1715 inline TYPE clipLower(const TYPE t,const TYPE valLow){ \
1716 return t < static_cast<TYPE>(valLow) ? static_cast<TYPE>(valLow) : t; \
1717 } \
1718 inline TYPE clipUpper(const TYPE t,const TYPE valHigh){ \
1719 return t > static_cast<TYPE>(valHigh) ? static_cast<TYPE>(valHigh) : t; \
1720 } \
1721 inline TYPE clip(const TYPE t,const TYPE valLow, const TYPE valHigh){ \
1722 if(t<valLow) \
1723 return valLow; \
1724 else if(t>valHigh) \
1725 return valHigh; \
1726 else \
1727 return t; \
1728 } \
1729 inline TYPE sum(const TYPE t){ \
1730 return t; \
1731 }\
1732 inline NumericTraits<TYPE>::RealPromote mean(const TYPE t){ \
1733 return t; \
1734 }\
1735 inline TYPE isZero(const TYPE t){ \
1736 return t==static_cast<TYPE>(0); \
1737 } \
1738 inline NumericTraits<TYPE>::RealPromote sizeDividedSquaredNorm(const TYPE t){ \
1739 return squaredNorm(t); \
1740 } \
1741 inline NumericTraits<TYPE>::RealPromote sizeDividedNorm(const TYPE t){ \
1742 return norm(t); \
1743 }
1744
1745
1746#ifdef __GNUC__
1747#pragma GCC diagnostic push
1748#pragma GCC diagnostic ignored "-Wtype-limits"
1749#endif
1750
1751VIGRA_MATH_FUNC_HELPER(unsigned char)
1752VIGRA_MATH_FUNC_HELPER(unsigned short)
1753VIGRA_MATH_FUNC_HELPER(unsigned int)
1754VIGRA_MATH_FUNC_HELPER(unsigned long)
1755VIGRA_MATH_FUNC_HELPER(unsigned long long)
1756VIGRA_MATH_FUNC_HELPER(signed char)
1757VIGRA_MATH_FUNC_HELPER(signed short)
1758VIGRA_MATH_FUNC_HELPER(signed int)
1759VIGRA_MATH_FUNC_HELPER(signed long)
1760VIGRA_MATH_FUNC_HELPER(signed long long)
1761VIGRA_MATH_FUNC_HELPER(float)
1762VIGRA_MATH_FUNC_HELPER(double)
1763VIGRA_MATH_FUNC_HELPER(long double)
1764
1765#ifdef __GNUC__
1766#pragma GCC diagnostic pop
1767#endif
1768
1769#undef VIGRA_MATH_FUNC_HELPER
1770
1771
1772} // namespace vigra
1773
1774#endif /* VIGRA_MATHUTIL_HXX */
Class for a single RGB value.
Definition rgbvalue.hxx:128
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition fixedpoint.hxx:667
double chi2(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Chi square distribution.
Definition mathutil.hxx:1042
detail::SelectIntegerType< 32, detail::UnsignedIntTypes >::type UInt32
32-bit unsigned int
Definition sized_int.hxx:183
int round(FixedPoint< IntBits, FracBits > v)
rounding to the nearest integer.
Definition fixedpoint.hxx:683
V power(const V &x)
Exponentiation to a positive integer power by squaring.
Definition mathutil.hxx:427
double gamma(double x)
The gamma function.
Definition mathutil.hxx:1587
R arg(const FFTWComplex< R > &a)
phase
Definition fftw3.hxx:1009
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
double ellipticIntegralF(double x, double k)
The incomplete elliptic integral of the first kind.
Definition mathutil.hxx:855
Int32 sqrti(Int32 v)
Signed integer square root.
Definition mathutil.hxx:538
double chi2CDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Cumulative chi square distribution.
Definition mathutil.hxx:1057
Int32 roundi(FixedPoint16< IntBits, OverflowHandling > v)
rounding to the nearest integer.
Definition fixedpoint.hxx:1775
UInt32 ceilPower2(UInt32 x)
Round up to the nearest power of 2.
Definition mathutil.hxx:294
REAL legendre(unsigned int l, int m, REAL x)
Associated Legendre polynomial.
Definition mathutil.hxx:1142
bool lessEqualAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point less-or-equal.
Definition mathutil.hxx:1672
bool even(int t)
Check if an integer is even.
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition tinyvector.hxx:2073
void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22, T *r0, T *r1, T *r2)
Compute the eigenvalues of a 3x3 real symmetric matrix.
Definition mathutil.hxx:754
UInt32 floorPower2(UInt32 x)
Round down to the nearest power of 2.
Definition mathutil.hxx:317
Int32 log2i(UInt32 x)
Compute the base-2 logarithm of an integer.
Definition mathutil.hxx:361
double loggamma(double x)
The natural logarithm of the gamma function.
Definition mathutil.hxx:1603
T sign(T t)
The sign function.
Definition mathutil.hxx:591
double noncentralChi2CDFApprox(unsigned int degreesOfFreedom, double noncentrality, double arg)
Cumulative non-central chi square distribution (approximate).
Definition mathutil.hxx:1111
void symmetric2x2Eigenvalues(T a00, T a01, T a11, T *r0, T *r1)
Compute the eigenvalues of a 2x2 real symmetric matrix.
Definition mathutil.hxx:734
double noncentralChi2(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Non-central chi square distribution.
Definition mathutil.hxx:1073
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point equality.
Definition mathutil.hxx:1638
double ellipticIntegralE(double x, double k)
The incomplete elliptic integral of the second kind.
Definition mathutil.hxx:879
bool greaterEqualAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point greater-or-equal.
Definition mathutil.hxx:1697
detail::SelectIntegerType< 32, detail::SignedIntTypes >::type Int32
32-bit signed int
Definition sized_int.hxx:175
int signi(T t)
The integer sign function.
Definition mathutil.hxx:608
double noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double accuracy=1e-7)
Cumulative non-central chi square distribution.
Definition mathutil.hxx:1091
FixedPoint16< IntBits, OverflowHandling > hypot(FixedPoint16< IntBits, OverflowHandling > v1, FixedPoint16< IntBits, OverflowHandling > v2)
Length of hypotenuse.
Definition fixedpoint.hxx:1640
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition fixedpoint.hxx:675
bool isPower2(UInt32 x)
Determine whether x is a power of 2 Bit twiddle from https://graphics.stanford.edu/~seander/bithacks....
Definition mathutil.hxx:278
bool odd(int t)
Check if an integer is odd.

© 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.11.2