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

mathutil.hxx VIGRA

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

© 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.0 (Thu Mar 17 2016)