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

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Tue Sep 24 2013)