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