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