1 /* 2 * Tiny Vector Matrix Library 3 * Dense Vector Matrix Libary of Tiny size using Expression Templates 4 * 5 * Copyright (C) 2001 - 2007 Olaf Petzold <opetzold@users.sourceforge.net> 6 * 7 * This library is free software; you can redistribute it and/or 8 * modify it under the terms of the GNU Lesser General Public 9 * License as published by the Free Software Foundation; either 10 * version 2.1 of the License, or (at your option) any later version. 11 * 12 * This library is distributed in the hope that it will be useful, 13 * but WITHOUT ANY WARRANTY; without even the implied warranty of 14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15 * Lesser General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General Public 18 * License along with this library; if not, write to the Free Software 19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 20 * 21 * $Id: NumericTraits.h,v 1.16 2007-06-23 15:58:58 opetzold Exp $ 22 */ 23 24 #ifndef TVMET_NUMERIC_TRAITS_H 25 #define TVMET_NUMERIC_TRAITS_H 26 27 #if defined(TVMET_HAVE_COMPLEX) 28 # include <complex> 29 #endif 30 #include <cmath> 31 #include <limits> 32 33 #include <tvmet/CompileTimeError.h> 34 35 36 namespace tvmet { 37 38 39 /** 40 * \class NumericTraits NumericTraits.h "tvmet/NumericTraits.h" 41 * \brief Traits for integral types for operations. 42 * 43 * For each type we have to specialize this traits. 44 * 45 * \note Keep in mind that the long types long long and long double doesn't 46 * have traits. This is due to the sum_type. We can't give a guarantee 47 * that there is a type of holding the sum. Therefore using this traits 48 * is only safe if you have long long resp. long double types by 49 * working on long ints and doubles. Otherwise you will get not expected 50 * result for some circumstances. Anyway, you can use big integer/float 51 * libraries and specialize the traits by your own. 52 * 53 * \todo The abs function of complex<non_float_type> can have an 54 * overrun due to numeric computation. Solve it (someone 55 * using value_type=long here?) 56 */ 57 template<class T> 58 struct NumericTraits { 59 typedef T base_type; 60 typedef T value_type; 61 typedef value_type sum_type; 62 typedef value_type diff_type; 63 typedef value_type float_type; 64 typedef value_type signed_type; 65 66 typedef NumericTraits<value_type> traits_type; 67 typedef const value_type& argument_type; 68 69 static inline 70 base_type real(argument_type x); 71 72 static inline 73 base_type imag(argument_type x); 74 75 static inline 76 value_type conj(argument_type x); 77 78 static inline 79 base_type abs(argument_type x); 80 81 static inline 82 value_type sqrt(argument_type x); 83 84 static inline norm_1NumericTraits85 base_type norm_1(argument_type x) { 86 return NumericTraits<base_type>::abs(traits_type::real(x)) 87 + NumericTraits<base_type>::abs(traits_type::imag(x)); 88 } 89 90 static inline norm_2NumericTraits91 base_type norm_2(argument_type x) { return traits_type::abs(x); } 92 93 static inline norm_infNumericTraits94 base_type norm_inf(argument_type x) { 95 return std::max(NumericTraits<base_type>::abs(traits_type::real(x)), 96 NumericTraits<base_type>::abs(traits_type::imag(x))); 97 } 98 99 static inline equalsNumericTraits100 bool equals(argument_type lhs, argument_type rhs) { 101 static base_type sqrt_epsilon( 102 NumericTraits<base_type>::sqrt( 103 std::numeric_limits<base_type>::epsilon())); 104 105 return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon * 106 std::max(std::max(traits_type::norm_inf(lhs), 107 traits_type::norm_inf(rhs)), 108 std::numeric_limits<base_type>::min()); 109 } 110 }; 111 112 113 /* 114 * numeric traits for standard types 115 */ 116 117 118 /** 119 * \class NumericTraits<char> NumericTraits.h "tvmet/NumericTraits.h" 120 * \brief Traits specialized for char. 121 */ 122 template<> 123 struct NumericTraits<char> { 124 typedef char value_type; 125 typedef value_type base_type; 126 typedef long sum_type; 127 typedef int diff_type; 128 typedef float float_type; 129 typedef char signed_type; 130 131 typedef NumericTraits<value_type> traits_type; 132 typedef value_type argument_type; 133 134 static inline 135 base_type real(argument_type x) { return x; } 136 137 static inline 138 base_type imag(argument_type) { return 0; } 139 140 static inline 141 value_type conj(argument_type x) { return x; } 142 143 static inline 144 base_type abs(argument_type x) { return std::abs(x); } 145 146 static inline 147 value_type sqrt(argument_type x) { 148 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 149 } 150 151 static inline 152 base_type norm_1(argument_type x) { return traits_type::abs(x); } 153 154 static inline 155 base_type norm_2(argument_type x) { return traits_type::abs(x); } 156 157 static inline 158 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 159 160 static inline 161 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 162 163 enum { is_complex = false }; 164 165 /** Complexity on operations. */ 166 enum { 167 ops_plus = 1, /**< Complexity on plus/minus ops. */ 168 ops_muls = 1 /**< Complexity on multiplications. */ 169 }; 170 }; 171 172 173 /** 174 * \class NumericTraits<unsigned char> NumericTraits.h "tvmet/NumericTraits.h" 175 * \brief Traits specialized for unsigned char. 176 * 177 * \note Normally it doesn't make sense to call <tt>conj</tt> 178 * for an unsigned type! An unary minus operator 179 * applied to unsigned type will result unsigned. Therefore 180 * this function is missing here. 181 */ 182 template<> 183 struct NumericTraits<unsigned char> { 184 typedef unsigned char value_type; 185 typedef value_type base_type; 186 typedef unsigned long sum_type; 187 typedef int diff_type; 188 typedef float float_type; 189 typedef int signed_type; 190 191 typedef NumericTraits<value_type> traits_type; 192 typedef value_type argument_type; 193 194 static inline 195 base_type real(argument_type x) { return x; } 196 197 static inline 198 base_type imag(argument_type) { return 0; } 199 200 static inline 201 base_type abs(argument_type x) { return std::abs(x); } 202 203 static inline 204 value_type sqrt(argument_type x) { 205 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 206 } 207 208 static inline 209 base_type norm_1(argument_type x) { return traits_type::abs(x); } 210 211 static inline 212 base_type norm_2(argument_type x) { return traits_type::abs(x); } 213 214 static inline 215 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 216 217 static inline 218 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 219 220 enum { is_complex = false }; 221 222 /** Complexity on operations. */ 223 enum { 224 ops_plus = 1, /**< Complexity on plus/minus ops. */ 225 ops_muls = 1 /**< Complexity on multiplications. */ 226 }; 227 }; 228 229 230 /** 231 * \class NumericTraits<short int> NumericTraits.h "tvmet/NumericTraits.h" 232 * \brief Traits specialized for short int. 233 */ 234 template<> 235 struct NumericTraits<short int> { 236 typedef short int value_type; 237 typedef value_type base_type; 238 #if defined(TVMET_HAVE_LONG_LONG) 239 typedef long long sum_type; 240 #else 241 typedef long sum_type; 242 #endif 243 typedef int diff_type; 244 typedef float float_type; 245 typedef short int signed_type; 246 247 typedef NumericTraits<value_type> traits_type; 248 typedef value_type argument_type; 249 250 static inline 251 base_type real(argument_type x) { return x; } 252 253 static inline 254 base_type imag(argument_type) { return 0; } 255 256 static inline 257 value_type conj(argument_type x) { return x; } 258 259 static inline 260 base_type abs(argument_type x) { return std::abs(x); } 261 262 static inline 263 value_type sqrt(argument_type x) { 264 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 265 } 266 267 static inline 268 base_type norm_1(argument_type x) { return traits_type::abs(x); } 269 270 static inline 271 base_type norm_2(argument_type x) { return traits_type::abs(x); } 272 273 static inline 274 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 275 276 static inline 277 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 278 279 enum { is_complex = false }; 280 281 /** Complexity on operations. */ 282 enum { 283 ops_plus = 1, /**< Complexity on plus/minus ops. */ 284 ops_muls = 1 /**< Complexity on multiplications. */ 285 }; 286 }; 287 288 289 /** 290 * \class NumericTraits<short unsigned int> NumericTraits.h "tvmet/NumericTraits.h" 291 * \brief Traits specialized for short unsigned int. 292 * 293 * \note Normally it doesn't make sense to call <tt>conj</tt> 294 * for an unsigned type! An unary minus operator 295 * applied to unsigned type will result unsigned. Therefore 296 * this function is missing here. 297 */ 298 template<> 299 struct NumericTraits<short unsigned int> { 300 typedef short unsigned int value_type; 301 typedef value_type base_type; 302 #if defined(TVMET_HAVE_LONG_LONG) 303 typedef unsigned long long sum_type; 304 #else 305 typedef unsigned long sum_type; 306 #endif 307 typedef int diff_type; 308 typedef float float_type; 309 typedef int signed_type; 310 311 typedef NumericTraits<value_type> traits_type; 312 typedef value_type argument_type; 313 314 static inline 315 base_type real(argument_type x) { return x; } 316 317 static inline 318 base_type imag(argument_type) { return 0; } 319 320 static inline 321 base_type abs(argument_type x) { return std::abs(x); } 322 323 static inline 324 value_type sqrt(argument_type x) { 325 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 326 } 327 328 static inline 329 base_type norm_1(argument_type x) { return traits_type::abs(x); } 330 331 static inline 332 base_type norm_2(argument_type x) { return traits_type::abs(x); } 333 334 static inline 335 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 336 337 static inline 338 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 339 340 enum { is_complex = false }; 341 342 /** Complexity on operations. */ 343 enum { 344 ops_plus = 1, /**< Complexity on plus/minus ops. */ 345 ops_muls = 1 /**< Complexity on multiplications. */ 346 }; 347 }; 348 349 350 /** 351 * \class NumericTraits<int> NumericTraits.h "tvmet/NumericTraits.h" 352 * \brief Traits specialized for int. 353 */ 354 template<> 355 struct NumericTraits<int> { 356 typedef int value_type; 357 typedef value_type base_type; 358 #if defined(TVMET_HAVE_LONG_LONG) 359 typedef long long sum_type; 360 #else 361 typedef long sum_type; 362 #endif 363 typedef int diff_type; 364 typedef double float_type; 365 typedef int signed_type; 366 367 typedef NumericTraits<value_type> traits_type; 368 typedef value_type argument_type; 369 370 static inline 371 base_type real(argument_type x) { return x; } 372 373 static inline 374 base_type imag(argument_type) { return 0; } 375 376 static inline 377 value_type conj(argument_type x) { return x; } 378 379 static inline 380 base_type abs(argument_type x) { return std::abs(x); } 381 382 static inline 383 value_type sqrt(argument_type x) { 384 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 385 } 386 387 static inline 388 base_type norm_1(argument_type x) { return traits_type::abs(x); } 389 390 static inline 391 base_type norm_2(argument_type x) { return traits_type::abs(x); } 392 393 static inline 394 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 395 396 static inline 397 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 398 399 enum { is_complex = false }; 400 401 /** Complexity on operations. */ 402 enum { 403 ops_plus = 1, /**< Complexity on plus/minus ops. */ 404 ops_muls = 1 /**< Complexity on multiplications. */ 405 }; 406 }; 407 408 409 /** 410 * \class NumericTraits<unsigned int> NumericTraits.h "tvmet/NumericTraits.h" 411 * \brief Traits specialized for unsigned int. 412 * 413 * \note Normally it doesn't make sense to call <tt>conj</tt> 414 * for an unsigned type! An unary minus operator 415 * applied to unsigned type will result unsigned. Therefore 416 * this function is missing here. 417 */ 418 template<> 419 struct NumericTraits<unsigned int> { 420 typedef unsigned int value_type; 421 typedef value_type base_type; 422 #if defined(TVMET_HAVE_LONG_LONG) 423 typedef unsigned long long sum_type; 424 #else 425 typedef unsigned long sum_type; 426 #endif 427 typedef int diff_type; 428 typedef double float_type; 429 typedef long signed_type; 430 431 typedef NumericTraits<value_type> traits_type; 432 typedef value_type argument_type; 433 434 static inline 435 base_type real(argument_type x) { return x; } 436 437 static inline 438 base_type imag(argument_type) { return 0; } 439 440 static inline 441 base_type abs(argument_type x) { return x; } 442 443 static inline 444 value_type sqrt(argument_type x) { 445 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 446 } 447 448 static inline 449 base_type norm_1(argument_type x) { return traits_type::abs(x); } 450 451 static inline 452 base_type norm_2(argument_type x) { return traits_type::abs(x); } 453 454 static inline 455 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 456 457 static inline 458 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 459 460 enum { is_complex = false }; 461 462 /** Complexity on operations. */ 463 enum { 464 ops_plus = 1, /**< Complexity on plus/minus ops. */ 465 ops_muls = 1 /**< Complexity on multiplications. */ 466 }; 467 }; 468 469 470 /** 471 * \class NumericTraits<long> NumericTraits.h "tvmet/NumericTraits.h" 472 * \brief Traits specialized for long. 473 */ 474 template<> 475 struct NumericTraits<long> { 476 typedef long value_type; 477 typedef value_type base_type; 478 #if defined(TVMET_HAVE_LONG_LONG) 479 typedef long long sum_type; 480 #else 481 typedef long sum_type; 482 #endif 483 typedef long diff_type; 484 typedef double float_type; 485 typedef long signed_type; 486 487 typedef NumericTraits<value_type> traits_type; 488 typedef value_type argument_type; 489 490 static inline 491 base_type real(argument_type x) { return x; } 492 493 static inline 494 base_type imag(argument_type) { return 0; } 495 496 static inline 497 value_type conj(argument_type x) { return x; } 498 499 static inline 500 base_type abs(argument_type x) { return std::abs(x); } 501 502 static inline 503 value_type sqrt(argument_type x) { 504 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 505 } 506 507 static inline 508 base_type norm_1(argument_type x) { return traits_type::abs(x); } 509 510 static inline 511 base_type norm_2(argument_type x) { return traits_type::abs(x); } 512 513 static inline 514 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 515 516 static inline 517 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 518 519 enum { is_complex = false }; 520 521 /** Complexity on operations. */ 522 enum { 523 ops_plus = 1, /**< Complexity on plus/minus ops. */ 524 ops_muls = 1 /**< Complexity on multiplications. */ 525 }; 526 }; 527 528 529 /** 530 * \class NumericTraits<unsigned long> NumericTraits.h "tvmet/NumericTraits.h" 531 * \brief Traits specialized for unsigned long. 532 * 533 * \note Normally it doesn't make sense to call <tt>conj</tt> 534 * for an unsigned type! An unary minus operator 535 * applied to unsigned type will result unsigned. Therefore 536 * this function is missing here. 537 */ 538 template<> 539 struct NumericTraits<unsigned long> { 540 typedef unsigned long value_type; 541 typedef value_type base_type; 542 #if defined(TVMET_HAVE_LONG_LONG) 543 typedef unsigned long long sum_type; 544 #else 545 typedef unsigned long sum_type; 546 #endif 547 typedef unsigned long diff_type; 548 typedef double float_type; 549 typedef long signed_type; 550 551 typedef NumericTraits<value_type> traits_type; 552 typedef value_type argument_type; 553 554 static inline 555 base_type real(argument_type x) { return x; } 556 557 static inline 558 base_type imag(argument_type) { return 0; } 559 560 static inline 561 base_type abs(argument_type x) { return x; } 562 563 static inline 564 value_type sqrt(argument_type x) { 565 return static_cast<value_type>(std::sqrt(static_cast<float_type>(x))); 566 } 567 568 static inline 569 base_type norm_1(argument_type x) { return traits_type::abs(x); } 570 571 static inline 572 base_type norm_2(argument_type x) { return traits_type::abs(x); } 573 574 static inline 575 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 576 577 static inline 578 bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; } 579 580 enum { is_complex = false }; 581 582 /** Complexity on operations. */ 583 enum { 584 ops_plus = 1, /**< Complexity on plus/minus ops. */ 585 ops_muls = 1 /**< Complexity on multiplications. */ 586 }; 587 }; 588 589 590 /** 591 * \class NumericTraits<float> NumericTraits.h "tvmet/NumericTraits.h" 592 * \brief Traits specialized for float. 593 */ 594 template<> 595 struct NumericTraits<float> { 596 typedef float value_type; 597 typedef value_type base_type; 598 typedef double sum_type; 599 typedef float diff_type; 600 typedef float float_type; 601 typedef float signed_type; 602 603 typedef NumericTraits<value_type> traits_type; 604 typedef value_type argument_type; 605 606 static inline 607 base_type real(argument_type x) { return x; } 608 609 static inline 610 base_type imag(argument_type) { return 0; } 611 612 static inline 613 value_type conj(argument_type x) { return x; } 614 615 static inline 616 base_type abs(argument_type x) { return std::abs(x); } 617 618 static inline 619 value_type sqrt(argument_type x) { return std::sqrt(x); } 620 621 static inline 622 base_type norm_1(argument_type x) { return traits_type::abs(x); } 623 624 static inline 625 base_type norm_2(argument_type x) { return traits_type::abs(x); } 626 627 static inline 628 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 629 630 static inline 631 bool equals(argument_type lhs, argument_type rhs) { 632 static base_type sqrt_epsilon( 633 NumericTraits<base_type>::sqrt( 634 std::numeric_limits<base_type>::epsilon())); 635 636 return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon * 637 std::max(std::max(traits_type::norm_inf(lhs), 638 traits_type::norm_inf(rhs)), 639 std::numeric_limits<base_type>::min()); 640 } 641 642 enum { is_complex = false }; 643 644 /** Complexity on operations. */ 645 enum { 646 ops_plus = 1, /**< Complexity on plus/minus ops. */ 647 ops_muls = 1 /**< Complexity on multiplications. */ 648 }; 649 }; 650 651 652 /** 653 * \class NumericTraits<double> NumericTraits.h "tvmet/NumericTraits.h" 654 * \brief Traits specialized for double. 655 */ 656 template<> 657 struct NumericTraits<double> { 658 typedef double value_type; 659 typedef value_type base_type; 660 #if defined(TVMET_HAVE_LONG_DOUBLE) 661 typedef long double sum_type; 662 #else 663 typedef double sum_type; 664 #endif 665 typedef double diff_type; 666 typedef double float_type; 667 typedef double signed_type; 668 669 typedef NumericTraits<value_type> traits_type; 670 typedef value_type argument_type; 671 672 static inline 673 base_type real(argument_type x) { return x; } 674 675 static inline 676 base_type imag(argument_type) { return 0; } 677 678 static inline 679 value_type conj(argument_type x) { return x; } 680 681 static inline 682 base_type abs(argument_type x) { return std::abs(x); } 683 684 static inline 685 value_type sqrt(argument_type x) { return std::sqrt(x); } 686 687 static inline 688 base_type norm_1(argument_type x) { return traits_type::abs(x); } 689 690 static inline 691 base_type norm_2(argument_type x) { return traits_type::abs(x); } 692 693 static inline 694 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 695 696 static inline 697 bool equals(argument_type lhs, argument_type rhs) { 698 static base_type sqrt_epsilon( 699 NumericTraits<base_type>::sqrt( 700 std::numeric_limits<base_type>::epsilon())); 701 702 return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon * 703 std::max(std::max(traits_type::norm_inf(lhs), 704 traits_type::norm_inf(rhs)), 705 std::numeric_limits<base_type>::min()); 706 } 707 708 enum { is_complex = false }; 709 710 /** Complexity on operations. */ 711 enum { 712 ops_plus = 1, /**< Complexity on plus/minus ops. */ 713 ops_muls = 1 /**< Complexity on multiplications. */ 714 }; 715 }; 716 717 718 #if defined(TVMET_HAVE_LONG_DOUBLE) 719 /** 720 * \class NumericTraits<long double> NumericTraits.h "tvmet/NumericTraits.h" 721 * \brief Traits specialized for long double. 722 */ 723 template<> 724 struct NumericTraits<long double> { 725 typedef long double value_type; 726 typedef value_type base_type; 727 typedef long double sum_type; 728 typedef long double diff_type; 729 typedef long double float_type; 730 typedef long double signed_type; 731 732 typedef NumericTraits<value_type> traits_type; 733 typedef value_type argument_type; 734 735 static inline 736 base_type real(argument_type x) { return x; } 737 738 static inline 739 base_type imag(argument_type) { return 0; } 740 741 static inline 742 value_type conj(argument_type x) { return x; } 743 744 static inline 745 base_type abs(argument_type x) { return std::abs(x); } 746 747 static inline 748 value_type sqrt(argument_type x) { return std::sqrt(x); } 749 750 static inline 751 base_type norm_1(argument_type x) { return traits_type::abs(x); } 752 753 static inline 754 base_type norm_2(argument_type x) { return traits_type::abs(x); } 755 756 static inline 757 base_type norm_inf(argument_type x) { return traits_type::abs(x); } 758 759 static inline 760 bool equals(argument_type lhs, argument_type rhs) { 761 static base_type sqrt_epsilon( 762 NumericTraits<base_type>::sqrt( 763 std::numeric_limits<base_type>::epsilon())); 764 765 return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon * 766 std::max(std::max(traits_type::norm_inf(lhs), 767 traits_type::norm_inf(rhs)), 768 std::numeric_limits<base_type>::min()); 769 } 770 771 enum { is_complex = false }; 772 773 /** Complexity on operations. */ 774 enum { 775 ops_plus = 1, /**< Complexity on plus/minus ops. */ 776 ops_muls = 1 /**< Complexity on multiplications. */ 777 }; 778 }; 779 #endif // TVMET_HAVE_LONG_DOUBLE 780 781 782 /* 783 * numeric traits for complex types 784 */ 785 #if defined(TVMET_HAVE_COMPLEX) 786 787 /** 788 * \class NumericTraits< std::complex<int> > NumericTraits.h "tvmet/NumericTraits.h" 789 * \brief Traits specialized for std::complex<int>. 790 */ 791 template<> 792 struct NumericTraits< std::complex<int> > { 793 typedef int base_type; 794 typedef std::complex<int> value_type; 795 typedef std::complex<long> sum_type; 796 typedef std::complex<int> diff_type; 797 typedef std::complex<float> float_type; 798 typedef std::complex<int> signed_type; 799 800 typedef NumericTraits<value_type> traits_type; 801 typedef const value_type& argument_type; 802 803 static inline 804 base_type real(argument_type z) { return std::real(z); } 805 806 static inline 807 base_type imag(argument_type z) { return std::imag(z); } 808 809 static inline 810 value_type conj(argument_type z) { return std::conj(z); } 811 812 static inline 813 base_type abs(argument_type z) { 814 base_type x = z.real(); 815 base_type y = z.imag(); 816 817 // XXX probably case of overrun; header complex uses scaling 818 return static_cast<base_type>(NumericTraits<base_type>::sqrt(x * x + y * y)); 819 } 820 821 static /* inline */ 822 value_type sqrt(argument_type z) { 823 // borrowed and adapted from header complex 824 base_type x = z.real(); 825 base_type y = z.imag(); 826 827 if(x == base_type()) { 828 base_type t = NumericTraits<base_type>::sqrt( 829 NumericTraits<base_type>::abs(y) / 2); 830 return value_type(t, y < base_type() ? -t : t); 831 } 832 else { 833 base_type t = NumericTraits<base_type>::sqrt( 834 2 * (traits_type::abs(z) 835 + NumericTraits<base_type>::abs(x))); 836 base_type u = t / 2; 837 return x > base_type() 838 ? value_type(u, y / t) 839 : value_type(NumericTraits<base_type>::abs(y) / t, y < base_type() ? -u : u); 840 } 841 } 842 843 static inline 844 base_type norm_1(argument_type z) { 845 return NumericTraits<base_type>::abs((traits_type::real(z))) 846 + NumericTraits<base_type>::abs((traits_type::imag(z))); 847 } 848 849 static inline 850 base_type norm_2(argument_type z) { return traits_type::abs(z); } 851 852 static inline 853 base_type norm_inf(argument_type z) { 854 return std::max(NumericTraits<base_type>::abs(traits_type::real(z)), 855 NumericTraits<base_type>::abs(traits_type::imag(z))); 856 } 857 858 static inline 859 bool equals(argument_type lhs, argument_type rhs) { 860 return (traits_type::real(lhs) == traits_type::real(rhs)) 861 && (traits_type::imag(lhs) == traits_type::imag(rhs)); 862 } 863 864 enum { is_complex = true }; 865 866 /** Complexity on operations. */ 867 enum { 868 ops_plus = 2, /**< Complexity on plus/minus ops. */ 869 ops_muls = 6 /**< Complexity on multiplications. */ 870 }; 871 }; 872 873 874 /** 875 * \class NumericTraits< std::complex<unsigned int> > NumericTraits.h "tvmet/NumericTraits.h" 876 * \brief Traits specialized for std::complex<unsigned int>. 877 * 878 * \note Normally it doesn't make sense to call <tt>conj</tt> 879 * for an unsigned type! An unary minus operator 880 * applied to unsigned type will result unsigned. Therefore 881 * this function is missing here. 882 */ 883 template<> 884 struct NumericTraits< std::complex<unsigned int> > { 885 typedef unsigned int base_type; 886 typedef std::complex<unsigned int> value_type; 887 typedef std::complex<unsigned long> sum_type; 888 typedef std::complex<int> diff_type; 889 typedef std::complex<float> float_type; 890 typedef std::complex<int> signed_type; 891 892 typedef NumericTraits<value_type> traits_type; 893 typedef const value_type& argument_type; 894 895 static inline 896 base_type real(argument_type z) { return std::real(z); } 897 898 static inline 899 base_type imag(argument_type z) { return std::imag(z); } 900 901 static inline 902 base_type abs(argument_type z) { 903 base_type x = z.real(); 904 base_type y = z.imag(); 905 906 // XXX probably case of overrun; header complex uses scaling 907 return static_cast<base_type>(NumericTraits<base_type>::sqrt(x * x + y * y)); 908 } 909 910 static /* inline */ 911 value_type sqrt(argument_type z) { 912 // borrowed and adapted from header complex 913 base_type x = z.real(); 914 base_type y = z.imag(); 915 916 if(x == base_type()) { 917 base_type t = NumericTraits<base_type>::sqrt( 918 NumericTraits<base_type>::abs(y) / 2); 919 return value_type(t, t); 920 } 921 else { 922 base_type t = NumericTraits<base_type>::sqrt( 923 2 * (traits_type::abs(z) 924 + NumericTraits<base_type>::abs(x))); 925 return value_type(t / 2, y / t); 926 } 927 } 928 929 static inline 930 base_type norm_1(argument_type z) { 931 return NumericTraits<base_type>::abs((traits_type::real(z))) 932 + NumericTraits<base_type>::abs((traits_type::imag(z))); 933 } 934 935 static inline 936 base_type norm_2(argument_type z) { return traits_type::abs(z); } 937 938 static inline 939 base_type norm_inf(argument_type z) { 940 return std::max(NumericTraits<base_type>::abs(traits_type::real(z)), 941 NumericTraits<base_type>::abs(traits_type::imag(z))); 942 } 943 944 static inline 945 bool equals(argument_type lhs, argument_type rhs) { 946 return (traits_type::real(lhs) == traits_type::real(rhs)) 947 && (traits_type::imag(lhs) == traits_type::imag(rhs)); 948 } 949 950 enum { is_complex = true }; 951 952 /** Complexity on operations. */ 953 enum { 954 ops_plus = 2, /**< Complexity on plus/minus ops. */ 955 ops_muls = 6 /**< Complexity on multiplications. */ 956 }; 957 }; 958 959 960 /** 961 * \class NumericTraits< std::complex<long> > NumericTraits.h "tvmet/NumericTraits.h" 962 * \brief Traits specialized for std::complex<long>. 963 */ 964 template<> 965 struct NumericTraits< std::complex<long> > { 966 typedef long base_type; 967 typedef std::complex<long> value_type; 968 #if defined(TVMET_HAVE_LONG_LONG) 969 typedef std::complex<long long> sum_type; 970 #else 971 typedef std::complex<long> sum_type; 972 #endif 973 typedef std::complex<int> diff_type; 974 typedef std::complex<float> float_type; 975 typedef std::complex<int> signed_type; 976 977 typedef NumericTraits<value_type> traits_type; 978 typedef const value_type& argument_type; 979 980 static inline 981 base_type real(argument_type z) { return std::real(z); } 982 983 static inline 984 base_type imag(argument_type z) { return std::imag(z); } 985 986 static inline 987 value_type conj(argument_type z) { return std::conj(z); } 988 989 static inline 990 base_type abs(argument_type z) { 991 base_type x = z.real(); 992 base_type y = z.imag(); 993 994 // XXX probably case of overrun; header complex uses scaling 995 return static_cast<base_type>(NumericTraits<base_type>::sqrt(x * x + y * y)); 996 } 997 998 static /* inline */ 999 value_type sqrt(argument_type z) { 1000 // borrowed and adapted from header complex 1001 base_type x = z.real(); 1002 base_type y = z.imag(); 1003 1004 if(x == base_type()) { 1005 base_type t = NumericTraits<base_type>::sqrt( 1006 NumericTraits<base_type>::abs(y) / 2); 1007 return value_type(t, y < base_type() ? -t : t); 1008 } 1009 else { 1010 base_type t = NumericTraits<base_type>::sqrt( 1011 2 * (traits_type::abs(z) 1012 + NumericTraits<base_type>::abs(x))); 1013 base_type u = t / 2; 1014 return x > base_type() 1015 ? value_type(u, y / t) 1016 : value_type(NumericTraits<base_type>::abs(y) / t, y < base_type() ? -u : u); 1017 } 1018 } 1019 1020 static inline 1021 base_type norm_1(argument_type z) { 1022 return NumericTraits<base_type>::abs((traits_type::real(z))) 1023 + NumericTraits<base_type>::abs((traits_type::imag(z))); 1024 } 1025 1026 static inline 1027 base_type norm_2(argument_type z) { return traits_type::abs(z); } 1028 1029 static inline 1030 base_type norm_inf(argument_type z) { 1031 return std::max(NumericTraits<base_type>::abs(traits_type::real(z)), 1032 NumericTraits<base_type>::abs(traits_type::imag(z))); 1033 } 1034 1035 static inline 1036 bool equals(argument_type lhs, argument_type rhs) { 1037 return (traits_type::real(lhs) == traits_type::real(rhs)) 1038 && (traits_type::imag(lhs) == traits_type::imag(rhs)); 1039 } 1040 1041 enum { is_complex = true }; 1042 1043 /** Complexity on operations. */ 1044 enum { 1045 ops_plus = 2, /**< Complexity on plus/minus ops. */ 1046 ops_muls = 6 /**< Complexity on multiplications. */ 1047 }; 1048 }; 1049 1050 1051 /** 1052 * \class NumericTraits< std::complex<unsigned long> > NumericTraits.h "tvmet/NumericTraits.h" 1053 * \brief Traits specialized for std::complex<unsigned long>. 1054 * 1055 * \note Normally it doesn't make sense to call <tt>conj</tt> 1056 * for an unsigned type! An unary minus operator 1057 * applied to unsigned type will result unsigned. Therefore 1058 * this function is missing here. 1059 */ 1060 template<> 1061 struct NumericTraits< std::complex<unsigned long> > { 1062 typedef unsigned long base_type; 1063 typedef std::complex<unsigned long> value_type; 1064 #if defined(TVMET_HAVE_LONG_LONG) 1065 typedef std::complex<unsigned long long> sum_type; 1066 #else 1067 typedef std::complex<unsigned long> sum_type; 1068 #endif 1069 typedef std::complex<long> diff_type; 1070 typedef std::complex<float> float_type; 1071 typedef std::complex<long> signed_type; 1072 1073 typedef NumericTraits<value_type> traits_type; 1074 typedef const value_type& argument_type; 1075 1076 static inline 1077 base_type real(argument_type z) { return std::real(z); } 1078 1079 static inline 1080 base_type imag(argument_type z) { return std::imag(z); } 1081 1082 static inline 1083 base_type abs(argument_type z) { 1084 base_type x = z.real(); 1085 base_type y = z.imag(); 1086 1087 // XXX probably case of overrun; header complex uses scaling 1088 return static_cast<base_type>(NumericTraits<base_type>::sqrt(x * x + y * y)); 1089 } 1090 1091 static /* inline */ 1092 value_type sqrt(argument_type z) { 1093 // borrowed and adapted from header complex 1094 base_type x = z.real(); 1095 base_type y = z.imag(); 1096 1097 if(x == base_type()) { 1098 base_type t = NumericTraits<base_type>::sqrt( 1099 NumericTraits<base_type>::abs(y) / 2); 1100 return value_type(t, t); 1101 } 1102 else { 1103 base_type t = NumericTraits<base_type>::sqrt( 1104 2 * (traits_type::abs(z) 1105 + NumericTraits<base_type>::abs(x))); 1106 return value_type(t / 2, y / t); 1107 } 1108 } 1109 1110 static inline 1111 base_type norm_1(argument_type z) { 1112 return NumericTraits<base_type>::abs((traits_type::real(z))) 1113 + NumericTraits<base_type>::abs((traits_type::imag(z))); 1114 } 1115 1116 static inline 1117 base_type norm_2(argument_type z) { return traits_type::abs(z); } 1118 1119 static inline 1120 base_type norm_inf(argument_type z) { 1121 return std::max(NumericTraits<base_type>::abs(traits_type::real(z)), 1122 NumericTraits<base_type>::abs(traits_type::imag(z))); 1123 } 1124 1125 static inline 1126 bool equals(argument_type lhs, argument_type rhs) { 1127 return (traits_type::real(lhs) == traits_type::real(rhs)) 1128 && (traits_type::imag(lhs) == traits_type::imag(rhs)); 1129 } 1130 1131 enum { is_complex = true }; 1132 1133 /** Complexity on operations.*/ 1134 enum { 1135 ops_plus = 2, /**< Complexity on plus/minus ops. */ 1136 ops_muls = 6 /**< Complexity on multiplications. */ 1137 }; 1138 }; 1139 1140 1141 /** 1142 * \class NumericTraits< std::complex<float> > NumericTraits.h "tvmet/NumericTraits.h" 1143 * \brief Traits specialized for std::complex<float>. 1144 */ 1145 template<> 1146 struct NumericTraits< std::complex<float> > { 1147 typedef float base_type; 1148 typedef std::complex<float> value_type; 1149 typedef std::complex<double> sum_type; 1150 typedef std::complex<float> diff_type; 1151 typedef std::complex<float> float_type; 1152 typedef std::complex<float> signed_type; 1153 1154 typedef NumericTraits<value_type> traits_type; 1155 typedef const value_type& argument_type; 1156 1157 static inline 1158 base_type real(argument_type z) { return std::real(z); } 1159 1160 static inline 1161 base_type imag(argument_type z) { return std::imag(z); } 1162 1163 static inline 1164 value_type conj(argument_type z) { return std::conj(z); } 1165 1166 static inline 1167 base_type abs(argument_type z) { return std::abs(z); } 1168 1169 static inline 1170 value_type sqrt(argument_type z) { return std::sqrt(z); } 1171 1172 static inline 1173 base_type norm_1(argument_type z) { 1174 return NumericTraits<base_type>::abs((traits_type::real(z))) 1175 + NumericTraits<base_type>::abs((traits_type::imag(z))); 1176 } 1177 1178 static inline 1179 base_type norm_2(argument_type z) { return traits_type::abs(z); } 1180 1181 static inline 1182 base_type norm_inf(argument_type z) { 1183 return std::max(NumericTraits<base_type>::abs(traits_type::real(z)), 1184 NumericTraits<base_type>::abs(traits_type::imag(z))); 1185 } 1186 1187 static inline 1188 bool equals(argument_type lhs, argument_type rhs) { 1189 static base_type sqrt_epsilon( 1190 NumericTraits<base_type>::sqrt( 1191 std::numeric_limits<base_type>::epsilon())); 1192 1193 return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon * 1194 std::max(std::max(traits_type::norm_inf(lhs), 1195 traits_type::norm_inf(rhs)), 1196 std::numeric_limits<base_type>::min()); 1197 } 1198 1199 enum { is_complex = true }; 1200 1201 /** Complexity on operations. */ 1202 enum { 1203 ops_plus = 2, /**< Complexity on plus/minus ops. */ 1204 ops_muls = 6 /**< Complexity on multiplications. */ 1205 }; 1206 }; 1207 1208 1209 /** 1210 * \class NumericTraits< std::complex<double> > NumericTraits.h "tvmet/NumericTraits.h" 1211 * \brief Traits specialized for std::complex<double>. 1212 */ 1213 template<> 1214 struct NumericTraits< std::complex<double> > { 1215 typedef double base_type; 1216 typedef std::complex<double> value_type; 1217 #if defined(TVMET_HAVE_LONG_DOUBLE) 1218 typedef std::complex<long double> sum_type; 1219 #else 1220 typedef std::complex<double> sum_type; 1221 #endif 1222 typedef std::complex<double> diff_type; 1223 typedef std::complex<double> float_type; 1224 typedef std::complex<double> signed_type; 1225 1226 typedef NumericTraits<value_type> traits_type; 1227 typedef const value_type& argument_type; 1228 1229 static inline 1230 base_type real(argument_type z) { return std::real(z); } 1231 1232 static inline 1233 base_type imag(argument_type z) { return std::imag(z); } 1234 1235 static inline 1236 value_type conj(argument_type z) { return std::conj(z); } 1237 1238 static inline 1239 base_type abs(argument_type z) { return std::abs(z); } 1240 1241 static inline 1242 value_type sqrt(argument_type z) { return std::sqrt(z); } 1243 1244 static inline 1245 base_type norm_1(argument_type z) { 1246 return NumericTraits<base_type>::abs((traits_type::real(z))) 1247 + NumericTraits<base_type>::abs((traits_type::imag(z))); 1248 } 1249 1250 static inline 1251 base_type norm_2(argument_type z) { return traits_type::abs(z); } 1252 1253 static inline 1254 base_type norm_inf(argument_type z) { 1255 return std::max(NumericTraits<base_type>::abs(traits_type::real(z)), 1256 NumericTraits<base_type>::abs(traits_type::imag(z))); 1257 } 1258 1259 static inline 1260 bool equals(argument_type lhs, argument_type rhs) { 1261 static base_type sqrt_epsilon( 1262 NumericTraits<base_type>::sqrt( 1263 std::numeric_limits<base_type>::epsilon())); 1264 1265 return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon * 1266 std::max(std::max(traits_type::norm_inf(lhs), 1267 traits_type::norm_inf(rhs)), 1268 std::numeric_limits<base_type>::min()); 1269 } 1270 1271 enum { is_complex = true }; 1272 1273 /** Complexity on operations. */ 1274 enum { 1275 ops_plus = 2, /**< Complexity on plus/minus ops. */ 1276 ops_muls = 6 /**< Complexity on multiplications. */ 1277 }; 1278 }; 1279 1280 1281 #if defined(TVMET_HAVE_LONG_DOUBLE) 1282 /** 1283 * \class NumericTraits< std::complex<long double> > NumericTraits.h "tvmet/NumericTraits.h" 1284 * \brief Traits specialized for std::complex<double>. 1285 */ 1286 template<> 1287 struct NumericTraits< std::complex<long double> > { 1288 typedef long double base_type; 1289 typedef std::complex<long double> value_type; 1290 typedef std::complex<long double> sum_type; 1291 typedef std::complex<long double> diff_type; 1292 typedef std::complex<long double> float_type; 1293 typedef std::complex<long double> signed_type; 1294 1295 typedef NumericTraits<value_type> traits_type; 1296 typedef const value_type& argument_type; 1297 1298 static inline 1299 base_type real(argument_type z) { return std::real(z); } 1300 1301 static inline 1302 base_type imag(argument_type z) { return std::imag(z); } 1303 1304 static inline 1305 value_type conj(argument_type z) { return std::conj(z); } 1306 1307 static inline 1308 base_type abs(argument_type z) { return std::abs(z); } 1309 1310 static inline 1311 value_type sqrt(argument_type z) { return std::sqrt(z); } 1312 1313 static inline 1314 base_type norm_1(argument_type z) { 1315 return NumericTraits<base_type>::abs((traits_type::real(z))) 1316 + NumericTraits<base_type>::abs((traits_type::imag(z))); 1317 } 1318 1319 static inline 1320 base_type norm_2(argument_type z) { return traits_type::abs(z); } 1321 1322 static inline 1323 base_type norm_inf(argument_type z) { 1324 return std::max(NumericTraits<base_type>::abs(traits_type::real(z)), 1325 NumericTraits<base_type>::abs(traits_type::imag(z))); 1326 } 1327 1328 static inline 1329 bool equals(argument_type lhs, argument_type rhs) { 1330 static base_type sqrt_epsilon( 1331 NumericTraits<base_type>::sqrt( 1332 std::numeric_limits<base_type>::epsilon())); 1333 1334 return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon * 1335 std::max(std::max(traits_type::norm_inf(lhs), 1336 traits_type::norm_inf(rhs)), 1337 std::numeric_limits<base_type>::min()); 1338 } 1339 1340 enum { is_complex = true }; 1341 1342 /** Complexity on operations. */ 1343 enum { 1344 ops_plus = 2, /**< Complexity on plus/minus ops. */ 1345 ops_muls = 6 /**< Complexity on multiplications. */ 1346 }; 1347 }; 1348 #endif // defined(TVMET_HAVE_LONG_DOUBLE) 1349 1350 1351 #endif // defined(TVMET_HAVE_COMPLEX) 1352 1353 1354 } // namespace tvmet 1355 1356 1357 #endif // TVMET_NUMERIC_TRAITS_H 1358 1359 1360 // Local Variables: 1361 // mode:C++ 1362 // tab-width:8 1363 // End: 1364