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