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