1 // boost quaternion.hpp header file 2 3 // (C) Copyright Hubert Holin 2001. 4 // Distributed under the Boost Software License, Version 1.0. (See 5 // accompanying file LICENSE_1_0.txt or copy at 6 // http://www.boost.org/LICENSE_1_0.txt) 7 8 // See http://www.boost.org for updates, documentation, and revision history. 9 10 #ifndef BOOST_QUATERNION_HPP 11 #define BOOST_QUATERNION_HPP 12 13 14 #include <complex> 15 #include <iosfwd> // for the "<<" and ">>" operators 16 #include <sstream> // for the "<<" operator 17 18 #include <boost/config.hpp> // for BOOST_NO_STD_LOCALE 19 #include <boost/detail/workaround.hpp> 20 #ifndef BOOST_NO_STD_LOCALE 21 #include <locale> // for the "<<" operator 22 #endif /* BOOST_NO_STD_LOCALE */ 23 24 #include <valarray> 25 26 27 28 #include <boost/math/special_functions/sinc.hpp> // for the Sinus cardinal 29 #include <boost/math/special_functions/sinhc.hpp> // for the Hyperbolic Sinus cardinal 30 31 32 namespace boost 33 { 34 namespace math 35 { 36 37 #define BOOST_QUATERNION_ACCESSOR_GENERATOR(type) \ 38 type real() const \ 39 { \ 40 return(a); \ 41 } \ 42 \ 43 quaternion<type> unreal() const \ 44 { \ 45 return(quaternion<type>(static_cast<type>(0),b,c,d)); \ 46 } \ 47 \ 48 type R_component_1() const \ 49 { \ 50 return(a); \ 51 } \ 52 \ 53 type R_component_2() const \ 54 { \ 55 return(b); \ 56 } \ 57 \ 58 type R_component_3() const \ 59 { \ 60 return(c); \ 61 } \ 62 \ 63 type R_component_4() const \ 64 { \ 65 return(d); \ 66 } \ 67 \ 68 ::std::complex<type> C_component_1() const \ 69 { \ 70 return(::std::complex<type>(a,b)); \ 71 } \ 72 \ 73 ::std::complex<type> C_component_2() const \ 74 { \ 75 return(::std::complex<type>(c,d)); \ 76 } 77 78 79 #define BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(type) \ 80 template<typename X> \ 81 quaternion<type> & operator = (quaternion<X> const & a_affecter) \ 82 { \ 83 a = static_cast<type>(a_affecter.R_component_1()); \ 84 b = static_cast<type>(a_affecter.R_component_2()); \ 85 c = static_cast<type>(a_affecter.R_component_3()); \ 86 d = static_cast<type>(a_affecter.R_component_4()); \ 87 \ 88 return(*this); \ 89 } \ 90 \ 91 quaternion<type> & operator = (quaternion<type> const & a_affecter) \ 92 { \ 93 a = a_affecter.a; \ 94 b = a_affecter.b; \ 95 c = a_affecter.c; \ 96 d = a_affecter.d; \ 97 \ 98 return(*this); \ 99 } \ 100 \ 101 quaternion<type> & operator = (type const & a_affecter) \ 102 { \ 103 a = a_affecter; \ 104 \ 105 b = c = d = static_cast<type>(0); \ 106 \ 107 return(*this); \ 108 } \ 109 \ 110 quaternion<type> & operator = (::std::complex<type> const & a_affecter) \ 111 { \ 112 a = a_affecter.real(); \ 113 b = a_affecter.imag(); \ 114 \ 115 c = d = static_cast<type>(0); \ 116 \ 117 return(*this); \ 118 } 119 120 121 #define BOOST_QUATERNION_MEMBER_DATA_GENERATOR(type) \ 122 type a; \ 123 type b; \ 124 type c; \ 125 type d; 126 127 128 template<typename T> 129 class quaternion 130 { 131 public: 132 133 typedef T value_type; 134 135 136 // constructor for H seen as R^4 137 // (also default constructor) 138 quaternion(T const & requested_a=T (),T const & requested_b=T (),T const & requested_c=T (),T const & requested_d=T ())139 explicit quaternion( T const & requested_a = T(), 140 T const & requested_b = T(), 141 T const & requested_c = T(), 142 T const & requested_d = T()) 143 : a(requested_a), 144 b(requested_b), 145 c(requested_c), 146 d(requested_d) 147 { 148 // nothing to do! 149 } 150 151 152 // constructor for H seen as C^2 153 quaternion(::std::complex<T> const & z0,::std::complex<T> const & z1=::std::complex<T> ())154 explicit quaternion( ::std::complex<T> const & z0, 155 ::std::complex<T> const & z1 = ::std::complex<T>()) 156 : a(z0.real()), 157 b(z0.imag()), 158 c(z1.real()), 159 d(z1.imag()) 160 { 161 // nothing to do! 162 } 163 164 165 // UNtemplated copy constructor 166 // (this is taken care of by the compiler itself) 167 168 169 // templated copy constructor 170 171 template<typename X> quaternion(quaternion<X> const & a_recopier)172 explicit quaternion(quaternion<X> const & a_recopier) 173 : a(static_cast<T>(a_recopier.R_component_1())), 174 b(static_cast<T>(a_recopier.R_component_2())), 175 c(static_cast<T>(a_recopier.R_component_3())), 176 d(static_cast<T>(a_recopier.R_component_4())) 177 { 178 // nothing to do! 179 } 180 181 182 // destructor 183 // (this is taken care of by the compiler itself) 184 185 186 // accessors 187 // 188 // Note: Like complex number, quaternions do have a meaningful notion of "real part", 189 // but unlike them there is no meaningful notion of "imaginary part". 190 // Instead there is an "unreal part" which itself is a quaternion, and usually 191 // nothing simpler (as opposed to the complex number case). 192 // However, for practicallity, there are accessors for the other components 193 // (these are necessary for the templated copy constructor, for instance). 194 195 BOOST_QUATERNION_ACCESSOR_GENERATOR(T) 196 197 // assignment operators 198 199 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(T) 200 201 // other assignment-related operators 202 // 203 // NOTE: Quaternion multiplication is *NOT* commutative; 204 // symbolically, "q *= rhs;" means "q = q * rhs;" 205 // and "q /= rhs;" means "q = q * inverse_of(rhs);" 206 207 quaternion<T> & operator += (T const & rhs) 208 { 209 T at = a + rhs; // exception guard 210 211 a = at; 212 213 return(*this); 214 } 215 216 operator +=(::std::complex<T> const & rhs)217 quaternion<T> & operator += (::std::complex<T> const & rhs) 218 { 219 T at = a + rhs.real(); // exception guard 220 T bt = b + rhs.imag(); // exception guard 221 222 a = at; 223 b = bt; 224 225 return(*this); 226 } 227 228 229 template<typename X> operator +=(quaternion<X> const & rhs)230 quaternion<T> & operator += (quaternion<X> const & rhs) 231 { 232 T at = a + static_cast<T>(rhs.R_component_1()); // exception guard 233 T bt = b + static_cast<T>(rhs.R_component_2()); // exception guard 234 T ct = c + static_cast<T>(rhs.R_component_3()); // exception guard 235 T dt = d + static_cast<T>(rhs.R_component_4()); // exception guard 236 237 a = at; 238 b = bt; 239 c = ct; 240 d = dt; 241 242 return(*this); 243 } 244 245 246 operator -=(T const & rhs)247 quaternion<T> & operator -= (T const & rhs) 248 { 249 T at = a - rhs; // exception guard 250 251 a = at; 252 253 return(*this); 254 } 255 256 operator -=(::std::complex<T> const & rhs)257 quaternion<T> & operator -= (::std::complex<T> const & rhs) 258 { 259 T at = a - rhs.real(); // exception guard 260 T bt = b - rhs.imag(); // exception guard 261 262 a = at; 263 b = bt; 264 265 return(*this); 266 } 267 268 269 template<typename X> operator -=(quaternion<X> const & rhs)270 quaternion<T> & operator -= (quaternion<X> const & rhs) 271 { 272 T at = a - static_cast<T>(rhs.R_component_1()); // exception guard 273 T bt = b - static_cast<T>(rhs.R_component_2()); // exception guard 274 T ct = c - static_cast<T>(rhs.R_component_3()); // exception guard 275 T dt = d - static_cast<T>(rhs.R_component_4()); // exception guard 276 277 a = at; 278 b = bt; 279 c = ct; 280 d = dt; 281 282 return(*this); 283 } 284 285 operator *=(T const & rhs)286 quaternion<T> & operator *= (T const & rhs) 287 { 288 T at = a * rhs; // exception guard 289 T bt = b * rhs; // exception guard 290 T ct = c * rhs; // exception guard 291 T dt = d * rhs; // exception guard 292 293 a = at; 294 b = bt; 295 c = ct; 296 d = dt; 297 298 return(*this); 299 } 300 301 operator *=(::std::complex<T> const & rhs)302 quaternion<T> & operator *= (::std::complex<T> const & rhs) 303 { 304 T ar = rhs.real(); 305 T br = rhs.imag(); 306 307 T at = +a*ar-b*br; 308 T bt = +a*br+b*ar; 309 T ct = +c*ar+d*br; 310 T dt = -c*br+d*ar; 311 312 a = at; 313 b = bt; 314 c = ct; 315 d = dt; 316 317 return(*this); 318 } 319 320 321 template<typename X> operator *=(quaternion<X> const & rhs)322 quaternion<T> & operator *= (quaternion<X> const & rhs) 323 { 324 T ar = static_cast<T>(rhs.R_component_1()); 325 T br = static_cast<T>(rhs.R_component_2()); 326 T cr = static_cast<T>(rhs.R_component_3()); 327 T dr = static_cast<T>(rhs.R_component_4()); 328 329 T at = +a*ar-b*br-c*cr-d*dr; 330 T bt = +a*br+b*ar+c*dr-d*cr; //(a*br+ar*b)+(c*dr-cr*d); 331 T ct = +a*cr-b*dr+c*ar+d*br; //(a*cr+ar*c)+(d*br-dr*b); 332 T dt = +a*dr+b*cr-c*br+d*ar; //(a*dr+ar*d)+(b*cr-br*c); 333 334 a = at; 335 b = bt; 336 c = ct; 337 d = dt; 338 339 return(*this); 340 } 341 342 343 operator /=(T const & rhs)344 quaternion<T> & operator /= (T const & rhs) 345 { 346 T at = a / rhs; // exception guard 347 T bt = b / rhs; // exception guard 348 T ct = c / rhs; // exception guard 349 T dt = d / rhs; // exception guard 350 351 a = at; 352 b = bt; 353 c = ct; 354 d = dt; 355 356 return(*this); 357 } 358 359 operator /=(::std::complex<T> const & rhs)360 quaternion<T> & operator /= (::std::complex<T> const & rhs) 361 { 362 T ar = rhs.real(); 363 T br = rhs.imag(); 364 365 T denominator = ar*ar+br*br; 366 367 T at = (+a*ar+b*br)/denominator; //(a*ar+b*br)/denominator; 368 T bt = (-a*br+b*ar)/denominator; //(ar*b-a*br)/denominator; 369 T ct = (+c*ar-d*br)/denominator; //(ar*c-d*br)/denominator; 370 T dt = (+c*br+d*ar)/denominator; //(ar*d+br*c)/denominator; 371 372 a = at; 373 b = bt; 374 c = ct; 375 d = dt; 376 377 return(*this); 378 } 379 380 381 template<typename X> operator /=(quaternion<X> const & rhs)382 quaternion<T> & operator /= (quaternion<X> const & rhs) 383 { 384 T ar = static_cast<T>(rhs.R_component_1()); 385 T br = static_cast<T>(rhs.R_component_2()); 386 T cr = static_cast<T>(rhs.R_component_3()); 387 T dr = static_cast<T>(rhs.R_component_4()); 388 389 T denominator = ar*ar+br*br+cr*cr+dr*dr; 390 391 T at = (+a*ar+b*br+c*cr+d*dr)/denominator; //(a*ar+b*br+c*cr+d*dr)/denominator; 392 T bt = (-a*br+b*ar-c*dr+d*cr)/denominator; //((ar*b-a*br)+(cr*d-c*dr))/denominator; 393 T ct = (-a*cr+b*dr+c*ar-d*br)/denominator; //((ar*c-a*cr)+(dr*b-d*br))/denominator; 394 T dt = (-a*dr-b*cr+c*br+d*ar)/denominator; //((ar*d-a*dr)+(br*c-b*cr))/denominator; 395 396 a = at; 397 b = bt; 398 c = ct; 399 d = dt; 400 401 return(*this); 402 } 403 404 405 protected: 406 407 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(T) 408 409 410 private: 411 412 }; 413 414 415 // declaration of quaternion specialization 416 417 template<> class quaternion<float>; 418 template<> class quaternion<double>; 419 template<> class quaternion<long double>; 420 421 422 // helper templates for converting copy constructors (declaration) 423 424 namespace detail 425 { 426 427 template< typename T, 428 typename U 429 > 430 quaternion<T> quaternion_type_converter(quaternion<U> const & rhs); 431 } 432 433 434 // implementation of quaternion specialization 435 436 437 #define BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(type) \ 438 explicit quaternion( type const & requested_a = static_cast<type>(0), \ 439 type const & requested_b = static_cast<type>(0), \ 440 type const & requested_c = static_cast<type>(0), \ 441 type const & requested_d = static_cast<type>(0)) \ 442 : a(requested_a), \ 443 b(requested_b), \ 444 c(requested_c), \ 445 d(requested_d) \ 446 { \ 447 } \ 448 \ 449 explicit quaternion( ::std::complex<type> const & z0, \ 450 ::std::complex<type> const & z1 = ::std::complex<type>()) \ 451 : a(z0.real()), \ 452 b(z0.imag()), \ 453 c(z1.real()), \ 454 d(z1.imag()) \ 455 { \ 456 } 457 458 459 #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type) \ 460 quaternion<type> & operator += (type const & rhs) \ 461 { \ 462 a += rhs; \ 463 \ 464 return(*this); \ 465 } 466 467 #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type) \ 468 quaternion<type> & operator += (::std::complex<type> const & rhs) \ 469 { \ 470 a += rhs.real(); \ 471 b += rhs.imag(); \ 472 \ 473 return(*this); \ 474 } 475 476 #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type) \ 477 template<typename X> \ 478 quaternion<type> & operator += (quaternion<X> const & rhs) \ 479 { \ 480 a += static_cast<type>(rhs.R_component_1()); \ 481 b += static_cast<type>(rhs.R_component_2()); \ 482 c += static_cast<type>(rhs.R_component_3()); \ 483 d += static_cast<type>(rhs.R_component_4()); \ 484 \ 485 return(*this); \ 486 } 487 488 #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type) \ 489 quaternion<type> & operator -= (type const & rhs) \ 490 { \ 491 a -= rhs; \ 492 \ 493 return(*this); \ 494 } 495 496 #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type) \ 497 quaternion<type> & operator -= (::std::complex<type> const & rhs) \ 498 { \ 499 a -= rhs.real(); \ 500 b -= rhs.imag(); \ 501 \ 502 return(*this); \ 503 } 504 505 #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type) \ 506 template<typename X> \ 507 quaternion<type> & operator -= (quaternion<X> const & rhs) \ 508 { \ 509 a -= static_cast<type>(rhs.R_component_1()); \ 510 b -= static_cast<type>(rhs.R_component_2()); \ 511 c -= static_cast<type>(rhs.R_component_3()); \ 512 d -= static_cast<type>(rhs.R_component_4()); \ 513 \ 514 return(*this); \ 515 } 516 517 #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1(type) \ 518 quaternion<type> & operator *= (type const & rhs) \ 519 { \ 520 a *= rhs; \ 521 b *= rhs; \ 522 c *= rhs; \ 523 d *= rhs; \ 524 \ 525 return(*this); \ 526 } 527 528 #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2(type) \ 529 quaternion<type> & operator *= (::std::complex<type> const & rhs) \ 530 { \ 531 type ar = rhs.real(); \ 532 type br = rhs.imag(); \ 533 \ 534 type at = +a*ar-b*br; \ 535 type bt = +a*br+b*ar; \ 536 type ct = +c*ar+d*br; \ 537 type dt = -c*br+d*ar; \ 538 \ 539 a = at; \ 540 b = bt; \ 541 c = ct; \ 542 d = dt; \ 543 \ 544 return(*this); \ 545 } 546 547 #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3(type) \ 548 template<typename X> \ 549 quaternion<type> & operator *= (quaternion<X> const & rhs) \ 550 { \ 551 type ar = static_cast<type>(rhs.R_component_1()); \ 552 type br = static_cast<type>(rhs.R_component_2()); \ 553 type cr = static_cast<type>(rhs.R_component_3()); \ 554 type dr = static_cast<type>(rhs.R_component_4()); \ 555 \ 556 type at = +a*ar-b*br-c*cr-d*dr; \ 557 type bt = +a*br+b*ar+c*dr-d*cr; \ 558 type ct = +a*cr-b*dr+c*ar+d*br; \ 559 type dt = +a*dr+b*cr-c*br+d*ar; \ 560 \ 561 a = at; \ 562 b = bt; \ 563 c = ct; \ 564 d = dt; \ 565 \ 566 return(*this); \ 567 } 568 569 // There is quite a lot of repetition in the code below. This is intentional. 570 // The last conditional block is the normal form, and the others merely 571 // consist of workarounds for various compiler deficiencies. Hopefuly, when 572 // more compilers are conformant and we can retire support for those that are 573 // not, we will be able to remove the clutter. This is makes the situation 574 // (painfully) explicit. 575 576 #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1(type) \ 577 quaternion<type> & operator /= (type const & rhs) \ 578 { \ 579 a /= rhs; \ 580 b /= rhs; \ 581 c /= rhs; \ 582 d /= rhs; \ 583 \ 584 return(*this); \ 585 } 586 587 #if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP) 588 #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ 589 quaternion<type> & operator /= (::std::complex<type> const & rhs) \ 590 { \ 591 using ::std::valarray; \ 592 using ::std::abs; \ 593 \ 594 valarray<type> tr(2); \ 595 \ 596 tr[0] = rhs.real(); \ 597 tr[1] = rhs.imag(); \ 598 \ 599 type mixam = static_cast<type>(1)/(abs(tr).max)(); \ 600 \ 601 tr *= mixam; \ 602 \ 603 valarray<type> tt(4); \ 604 \ 605 tt[0] = +a*tr[0]+b*tr[1]; \ 606 tt[1] = -a*tr[1]+b*tr[0]; \ 607 tt[2] = +c*tr[0]-d*tr[1]; \ 608 tt[3] = +c*tr[1]+d*tr[0]; \ 609 \ 610 tr *= tr; \ 611 \ 612 tt *= (mixam/tr.sum()); \ 613 \ 614 a = tt[0]; \ 615 b = tt[1]; \ 616 c = tt[2]; \ 617 d = tt[3]; \ 618 \ 619 return(*this); \ 620 } 621 #else 622 #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ 623 quaternion<type> & operator /= (::std::complex<type> const & rhs) \ 624 { \ 625 using ::std::valarray; \ 626 \ 627 valarray<type> tr(2); \ 628 \ 629 tr[0] = rhs.real(); \ 630 tr[1] = rhs.imag(); \ 631 \ 632 type mixam = static_cast<type>(1)/(abs(tr).max)(); \ 633 \ 634 tr *= mixam; \ 635 \ 636 valarray<type> tt(4); \ 637 \ 638 tt[0] = +a*tr[0]+b*tr[1]; \ 639 tt[1] = -a*tr[1]+b*tr[0]; \ 640 tt[2] = +c*tr[0]-d*tr[1]; \ 641 tt[3] = +c*tr[1]+d*tr[0]; \ 642 \ 643 tr *= tr; \ 644 \ 645 tt *= (mixam/tr.sum()); \ 646 \ 647 a = tt[0]; \ 648 b = tt[1]; \ 649 c = tt[2]; \ 650 d = tt[3]; \ 651 \ 652 return(*this); \ 653 } 654 #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ 655 656 #if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP) 657 #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) \ 658 template<typename X> \ 659 quaternion<type> & operator /= (quaternion<X> const & rhs) \ 660 { \ 661 using ::std::valarray; \ 662 using ::std::abs; \ 663 \ 664 valarray<type> tr(4); \ 665 \ 666 tr[0] = static_cast<type>(rhs.R_component_1()); \ 667 tr[1] = static_cast<type>(rhs.R_component_2()); \ 668 tr[2] = static_cast<type>(rhs.R_component_3()); \ 669 tr[3] = static_cast<type>(rhs.R_component_4()); \ 670 \ 671 type mixam = static_cast<type>(1)/(abs(tr).max)(); \ 672 \ 673 tr *= mixam; \ 674 \ 675 valarray<type> tt(4); \ 676 \ 677 tt[0] = +a*tr[0]+b*tr[1]+c*tr[2]+d*tr[3]; \ 678 tt[1] = -a*tr[1]+b*tr[0]-c*tr[3]+d*tr[2]; \ 679 tt[2] = -a*tr[2]+b*tr[3]+c*tr[0]-d*tr[1]; \ 680 tt[3] = -a*tr[3]-b*tr[2]+c*tr[1]+d*tr[0]; \ 681 \ 682 tr *= tr; \ 683 \ 684 tt *= (mixam/tr.sum()); \ 685 \ 686 a = tt[0]; \ 687 b = tt[1]; \ 688 c = tt[2]; \ 689 d = tt[3]; \ 690 \ 691 return(*this); \ 692 } 693 #else 694 #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) \ 695 template<typename X> \ 696 quaternion<type> & operator /= (quaternion<X> const & rhs) \ 697 { \ 698 using ::std::valarray; \ 699 \ 700 valarray<type> tr(4); \ 701 \ 702 tr[0] = static_cast<type>(rhs.R_component_1()); \ 703 tr[1] = static_cast<type>(rhs.R_component_2()); \ 704 tr[2] = static_cast<type>(rhs.R_component_3()); \ 705 tr[3] = static_cast<type>(rhs.R_component_4()); \ 706 \ 707 type mixam = static_cast<type>(1)/(abs(tr).max)(); \ 708 \ 709 tr *= mixam; \ 710 \ 711 valarray<type> tt(4); \ 712 \ 713 tt[0] = +a*tr[0]+b*tr[1]+c*tr[2]+d*tr[3]; \ 714 tt[1] = -a*tr[1]+b*tr[0]-c*tr[3]+d*tr[2]; \ 715 tt[2] = -a*tr[2]+b*tr[3]+c*tr[0]-d*tr[1]; \ 716 tt[3] = -a*tr[3]-b*tr[2]+c*tr[1]+d*tr[0]; \ 717 \ 718 tr *= tr; \ 719 \ 720 tt *= (mixam/tr.sum()); \ 721 \ 722 a = tt[0]; \ 723 b = tt[1]; \ 724 c = tt[2]; \ 725 d = tt[3]; \ 726 \ 727 return(*this); \ 728 } 729 #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ 730 731 #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR(type) \ 732 BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type) \ 733 BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type) \ 734 BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type) 735 736 #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR(type) \ 737 BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type) \ 738 BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type) \ 739 BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type) 740 741 #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR(type) \ 742 BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1(type) \ 743 BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2(type) \ 744 BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3(type) 745 746 #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR(type) \ 747 BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1(type) \ 748 BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ 749 BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) 750 751 #define BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(type) \ 752 BOOST_QUATERNION_MEMBER_ADD_GENERATOR(type) \ 753 BOOST_QUATERNION_MEMBER_SUB_GENERATOR(type) \ 754 BOOST_QUATERNION_MEMBER_MUL_GENERATOR(type) \ 755 BOOST_QUATERNION_MEMBER_DIV_GENERATOR(type) 756 757 758 template<> 759 class quaternion<float> 760 { 761 public: 762 763 typedef float value_type; 764 765 BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(float) 766 767 // UNtemplated copy constructor 768 // (this is taken care of by the compiler itself) 769 770 // explicit copy constructors (precision-loosing converters) 771 quaternion(quaternion<double> const & a_recopier)772 explicit quaternion(quaternion<double> const & a_recopier) 773 { 774 *this = detail::quaternion_type_converter<float, double>(a_recopier); 775 } 776 quaternion(quaternion<long double> const & a_recopier)777 explicit quaternion(quaternion<long double> const & a_recopier) 778 { 779 *this = detail::quaternion_type_converter<float, long double>(a_recopier); 780 } 781 782 // destructor 783 // (this is taken care of by the compiler itself) 784 785 // accessors 786 // 787 // Note: Like complex number, quaternions do have a meaningful notion of "real part", 788 // but unlike them there is no meaningful notion of "imaginary part". 789 // Instead there is an "unreal part" which itself is a quaternion, and usually 790 // nothing simpler (as opposed to the complex number case). 791 // However, for practicallity, there are accessors for the other components 792 // (these are necessary for the templated copy constructor, for instance). 793 794 BOOST_QUATERNION_ACCESSOR_GENERATOR(float) 795 796 // assignment operators 797 798 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(float) 799 800 // other assignment-related operators 801 // 802 // NOTE: Quaternion multiplication is *NOT* commutative; 803 // symbolically, "q *= rhs;" means "q = q * rhs;" 804 // and "q /= rhs;" means "q = q * inverse_of(rhs);" 805 806 BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(float) 807 808 809 protected: 810 811 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(float) 812 813 814 private: 815 816 }; 817 818 819 template<> 820 class quaternion<double> 821 { 822 public: 823 824 typedef double value_type; 825 826 BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(double) 827 828 // UNtemplated copy constructor 829 // (this is taken care of by the compiler itself) 830 831 // converting copy constructor 832 quaternion(quaternion<float> const & a_recopier)833 explicit quaternion(quaternion<float> const & a_recopier) 834 { 835 *this = detail::quaternion_type_converter<double, float>(a_recopier); 836 } 837 838 // explicit copy constructors (precision-loosing converters) 839 quaternion(quaternion<long double> const & a_recopier)840 explicit quaternion(quaternion<long double> const & a_recopier) 841 { 842 *this = detail::quaternion_type_converter<double, long double>(a_recopier); 843 } 844 845 // destructor 846 // (this is taken care of by the compiler itself) 847 848 // accessors 849 // 850 // Note: Like complex number, quaternions do have a meaningful notion of "real part", 851 // but unlike them there is no meaningful notion of "imaginary part". 852 // Instead there is an "unreal part" which itself is a quaternion, and usually 853 // nothing simpler (as opposed to the complex number case). 854 // However, for practicallity, there are accessors for the other components 855 // (these are necessary for the templated copy constructor, for instance). 856 857 BOOST_QUATERNION_ACCESSOR_GENERATOR(double) 858 859 // assignment operators 860 861 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(double) 862 863 // other assignment-related operators 864 // 865 // NOTE: Quaternion multiplication is *NOT* commutative; 866 // symbolically, "q *= rhs;" means "q = q * rhs;" 867 // and "q /= rhs;" means "q = q * inverse_of(rhs);" 868 869 BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(double) 870 871 872 protected: 873 874 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(double) 875 876 877 private: 878 879 }; 880 881 882 template<> 883 class quaternion<long double> 884 { 885 public: 886 887 typedef long double value_type; 888 889 BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(long double) 890 891 // UNtemplated copy constructor 892 // (this is taken care of by the compiler itself) 893 894 // converting copy constructors 895 quaternion(quaternion<float> const & a_recopier)896 explicit quaternion(quaternion<float> const & a_recopier) 897 { 898 *this = detail::quaternion_type_converter<long double, float>(a_recopier); 899 } 900 quaternion(quaternion<double> const & a_recopier)901 explicit quaternion(quaternion<double> const & a_recopier) 902 { 903 *this = detail::quaternion_type_converter<long double, double>(a_recopier); 904 } 905 906 // destructor 907 // (this is taken care of by the compiler itself) 908 909 // accessors 910 // 911 // Note: Like complex number, quaternions do have a meaningful notion of "real part", 912 // but unlike them there is no meaningful notion of "imaginary part". 913 // Instead there is an "unreal part" which itself is a quaternion, and usually 914 // nothing simpler (as opposed to the complex number case). 915 // However, for practicallity, there are accessors for the other components 916 // (these are necessary for the templated copy constructor, for instance). 917 918 BOOST_QUATERNION_ACCESSOR_GENERATOR(long double) 919 920 // assignment operators 921 922 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(long double) 923 924 // other assignment-related operators 925 // 926 // NOTE: Quaternion multiplication is *NOT* commutative; 927 // symbolically, "q *= rhs;" means "q = q * rhs;" 928 // and "q /= rhs;" means "q = q * inverse_of(rhs);" 929 930 BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(long double) 931 932 933 protected: 934 935 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(long double) 936 937 938 private: 939 940 }; 941 942 943 #undef BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR 944 #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR 945 #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR 946 #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR 947 #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR 948 #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1 949 #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2 950 #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3 951 #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1 952 #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2 953 #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3 954 #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1 955 #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2 956 #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3 957 #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1 958 #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2 959 #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3 960 961 #undef BOOST_QUATERNION_CONSTRUCTOR_GENERATOR 962 963 964 #undef BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR 965 966 #undef BOOST_QUATERNION_MEMBER_DATA_GENERATOR 967 968 #undef BOOST_QUATERNION_ACCESSOR_GENERATOR 969 970 971 // operators 972 973 #define BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) \ 974 { \ 975 quaternion<T> res(lhs); \ 976 res op##= rhs; \ 977 return(res); \ 978 } 979 980 #define BOOST_QUATERNION_OPERATOR_GENERATOR_1_L(op) \ 981 template<typename T> \ 982 inline quaternion<T> operator op (T const & lhs, quaternion<T> const & rhs) \ 983 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) 984 985 #define BOOST_QUATERNION_OPERATOR_GENERATOR_1_R(op) \ 986 template<typename T> \ 987 inline quaternion<T> operator op (quaternion<T> const & lhs, T const & rhs) \ 988 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) 989 990 #define BOOST_QUATERNION_OPERATOR_GENERATOR_2_L(op) \ 991 template<typename T> \ 992 inline quaternion<T> operator op (::std::complex<T> const & lhs, quaternion<T> const & rhs) \ 993 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) 994 995 #define BOOST_QUATERNION_OPERATOR_GENERATOR_2_R(op) \ 996 template<typename T> \ 997 inline quaternion<T> operator op (quaternion<T> const & lhs, ::std::complex<T> const & rhs) \ 998 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) 999 1000 #define BOOST_QUATERNION_OPERATOR_GENERATOR_3(op) \ 1001 template<typename T> \ 1002 inline quaternion<T> operator op (quaternion<T> const & lhs, quaternion<T> const & rhs) \ 1003 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) 1004 1005 #define BOOST_QUATERNION_OPERATOR_GENERATOR(op) \ 1006 BOOST_QUATERNION_OPERATOR_GENERATOR_1_L(op) \ 1007 BOOST_QUATERNION_OPERATOR_GENERATOR_1_R(op) \ 1008 BOOST_QUATERNION_OPERATOR_GENERATOR_2_L(op) \ 1009 BOOST_QUATERNION_OPERATOR_GENERATOR_2_R(op) \ 1010 BOOST_QUATERNION_OPERATOR_GENERATOR_3(op) 1011 1012 1013 BOOST_QUATERNION_OPERATOR_GENERATOR(+) 1014 BOOST_QUATERNION_OPERATOR_GENERATOR(-) 1015 BOOST_QUATERNION_OPERATOR_GENERATOR(*) 1016 BOOST_QUATERNION_OPERATOR_GENERATOR(/) 1017 1018 1019 #undef BOOST_QUATERNION_OPERATOR_GENERATOR 1020 1021 #undef BOOST_QUATERNION_OPERATOR_GENERATOR_1_L 1022 #undef BOOST_QUATERNION_OPERATOR_GENERATOR_1_R 1023 #undef BOOST_QUATERNION_OPERATOR_GENERATOR_2_L 1024 #undef BOOST_QUATERNION_OPERATOR_GENERATOR_2_R 1025 #undef BOOST_QUATERNION_OPERATOR_GENERATOR_3 1026 1027 #undef BOOST_QUATERNION_OPERATOR_GENERATOR_BODY 1028 1029 1030 template<typename T> operator +(quaternion<T> const & q)1031 inline quaternion<T> operator + (quaternion<T> const & q) 1032 { 1033 return(q); 1034 } 1035 1036 1037 template<typename T> operator -(quaternion<T> const & q)1038 inline quaternion<T> operator - (quaternion<T> const & q) 1039 { 1040 return(quaternion<T>(-q.R_component_1(),-q.R_component_2(),-q.R_component_3(),-q.R_component_4())); 1041 } 1042 1043 1044 template<typename T> operator ==(T const & lhs,quaternion<T> const & rhs)1045 inline bool operator == (T const & lhs, quaternion<T> const & rhs) 1046 { 1047 return ( 1048 (rhs.R_component_1() == lhs)&& 1049 (rhs.R_component_2() == static_cast<T>(0))&& 1050 (rhs.R_component_3() == static_cast<T>(0))&& 1051 (rhs.R_component_4() == static_cast<T>(0)) 1052 ); 1053 } 1054 1055 1056 template<typename T> operator ==(quaternion<T> const & lhs,T const & rhs)1057 inline bool operator == (quaternion<T> const & lhs, T const & rhs) 1058 { 1059 return ( 1060 (lhs.R_component_1() == rhs)&& 1061 (lhs.R_component_2() == static_cast<T>(0))&& 1062 (lhs.R_component_3() == static_cast<T>(0))&& 1063 (lhs.R_component_4() == static_cast<T>(0)) 1064 ); 1065 } 1066 1067 1068 template<typename T> operator ==(::std::complex<T> const & lhs,quaternion<T> const & rhs)1069 inline bool operator == (::std::complex<T> const & lhs, quaternion<T> const & rhs) 1070 { 1071 return ( 1072 (rhs.R_component_1() == lhs.real())&& 1073 (rhs.R_component_2() == lhs.imag())&& 1074 (rhs.R_component_3() == static_cast<T>(0))&& 1075 (rhs.R_component_4() == static_cast<T>(0)) 1076 ); 1077 } 1078 1079 1080 template<typename T> operator ==(quaternion<T> const & lhs,::std::complex<T> const & rhs)1081 inline bool operator == (quaternion<T> const & lhs, ::std::complex<T> const & rhs) 1082 { 1083 return ( 1084 (lhs.R_component_1() == rhs.real())&& 1085 (lhs.R_component_2() == rhs.imag())&& 1086 (lhs.R_component_3() == static_cast<T>(0))&& 1087 (lhs.R_component_4() == static_cast<T>(0)) 1088 ); 1089 } 1090 1091 1092 template<typename T> operator ==(quaternion<T> const & lhs,quaternion<T> const & rhs)1093 inline bool operator == (quaternion<T> const & lhs, quaternion<T> const & rhs) 1094 { 1095 return ( 1096 (rhs.R_component_1() == lhs.R_component_1())&& 1097 (rhs.R_component_2() == lhs.R_component_2())&& 1098 (rhs.R_component_3() == lhs.R_component_3())&& 1099 (rhs.R_component_4() == lhs.R_component_4()) 1100 ); 1101 } 1102 1103 1104 #define BOOST_QUATERNION_NOT_EQUAL_GENERATOR \ 1105 { \ 1106 return(!(lhs == rhs)); \ 1107 } 1108 1109 template<typename T> 1110 inline bool operator != (T const & lhs, quaternion<T> const & rhs) 1111 BOOST_QUATERNION_NOT_EQUAL_GENERATOR 1112 1113 template<typename T> 1114 inline bool operator != (quaternion<T> const & lhs, T const & rhs) 1115 BOOST_QUATERNION_NOT_EQUAL_GENERATOR 1116 1117 template<typename T> 1118 inline bool operator != (::std::complex<T> const & lhs, quaternion<T> const & rhs) 1119 BOOST_QUATERNION_NOT_EQUAL_GENERATOR 1120 1121 template<typename T> 1122 inline bool operator != (quaternion<T> const & lhs, ::std::complex<T> const & rhs) 1123 BOOST_QUATERNION_NOT_EQUAL_GENERATOR 1124 1125 template<typename T> 1126 inline bool operator != (quaternion<T> const & lhs, quaternion<T> const & rhs) 1127 BOOST_QUATERNION_NOT_EQUAL_GENERATOR 1128 1129 #undef BOOST_QUATERNION_NOT_EQUAL_GENERATOR 1130 1131 1132 // Note: we allow the following formats, whith a, b, c, and d reals 1133 // a 1134 // (a), (a,b), (a,b,c), (a,b,c,d) 1135 // (a,(c)), (a,(c,d)), ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b),(c,d)) 1136 template<typename T, typename charT, class traits> 1137 ::std::basic_istream<charT,traits> & operator >> ( ::std::basic_istream<charT,traits> & is, 1138 quaternion<T> & q) 1139 { 1140 1141 #ifdef BOOST_NO_STD_LOCALE 1142 #else 1143 const ::std::ctype<charT> & ct = ::std::use_facet< ::std::ctype<charT> >(is.getloc()); 1144 #endif /* BOOST_NO_STD_LOCALE */ 1145 1146 T a = T(); 1147 T b = T(); 1148 T c = T(); 1149 T d = T(); 1150 1151 ::std::complex<T> u = ::std::complex<T>(); 1152 ::std::complex<T> v = ::std::complex<T>(); 1153 1154 charT ch = charT(); 1155 char cc; 1156 1157 is >> ch; // get the first lexeme 1158 1159 if (!is.good()) goto finish; 1160 1161 #ifdef BOOST_NO_STD_LOCALE 1162 cc = ch; 1163 #else 1164 cc = ct.narrow(ch, char()); 1165 #endif /* BOOST_NO_STD_LOCALE */ 1166 1167 if (cc == '(') // read "(", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d)), ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,)) 1168 { 1169 is >> ch; // get the second lexeme 1170 1171 if (!is.good()) goto finish; 1172 1173 #ifdef BOOST_NO_STD_LOCALE 1174 cc = ch; 1175 #else 1176 cc = ct.narrow(ch, char()); 1177 #endif /* BOOST_NO_STD_LOCALE */ 1178 1179 if (cc == '(') // read "((", possible: ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,)) 1180 { 1181 is.putback(ch); 1182 1183 is >> u; // we extract the first and second components 1184 a = u.real(); 1185 b = u.imag(); 1186 1187 if (!is.good()) goto finish; 1188 1189 is >> ch; // get the next lexeme 1190 1191 if (!is.good()) goto finish; 1192 1193 #ifdef BOOST_NO_STD_LOCALE 1194 cc = ch; 1195 #else 1196 cc = ct.narrow(ch, char()); 1197 #endif /* BOOST_NO_STD_LOCALE */ 1198 1199 if (cc == ')') // format: ((a)) or ((a,b)) 1200 { 1201 q = quaternion<T>(a,b); 1202 } 1203 else if (cc == ',') // read "((a)," or "((a,b),", possible: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,)) 1204 { 1205 is >> v; // we extract the third and fourth components 1206 c = v.real(); 1207 d = v.imag(); 1208 1209 if (!is.good()) goto finish; 1210 1211 is >> ch; // get the last lexeme 1212 1213 if (!is.good()) goto finish; 1214 1215 #ifdef BOOST_NO_STD_LOCALE 1216 cc = ch; 1217 #else 1218 cc = ct.narrow(ch, char()); 1219 #endif /* BOOST_NO_STD_LOCALE */ 1220 1221 if (cc == ')') // format: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)) or ((a,b,),(c,d,)) 1222 { 1223 q = quaternion<T>(a,b,c,d); 1224 } 1225 else // error 1226 { 1227 is.setstate(::std::ios_base::failbit); 1228 } 1229 } 1230 else // error 1231 { 1232 is.setstate(::std::ios_base::failbit); 1233 } 1234 } 1235 else // read "(a", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d)) 1236 { 1237 is.putback(ch); 1238 1239 is >> a; // we extract the first component 1240 1241 if (!is.good()) goto finish; 1242 1243 is >> ch; // get the third lexeme 1244 1245 if (!is.good()) goto finish; 1246 1247 #ifdef BOOST_NO_STD_LOCALE 1248 cc = ch; 1249 #else 1250 cc = ct.narrow(ch, char()); 1251 #endif /* BOOST_NO_STD_LOCALE */ 1252 1253 if (cc == ')') // format: (a) 1254 { 1255 q = quaternion<T>(a); 1256 } 1257 else if (cc == ',') // read "(a,", possible: (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d)) 1258 { 1259 is >> ch; // get the fourth lexeme 1260 1261 if (!is.good()) goto finish; 1262 1263 #ifdef BOOST_NO_STD_LOCALE 1264 cc = ch; 1265 #else 1266 cc = ct.narrow(ch, char()); 1267 #endif /* BOOST_NO_STD_LOCALE */ 1268 1269 if (cc == '(') // read "(a,(", possible: (a,(c)), (a,(c,d)) 1270 { 1271 is.putback(ch); 1272 1273 is >> v; // we extract the third and fourth component 1274 1275 c = v.real(); 1276 d = v.imag(); 1277 1278 if (!is.good()) goto finish; 1279 1280 is >> ch; // get the ninth lexeme 1281 1282 if (!is.good()) goto finish; 1283 1284 #ifdef BOOST_NO_STD_LOCALE 1285 cc = ch; 1286 #else 1287 cc = ct.narrow(ch, char()); 1288 #endif /* BOOST_NO_STD_LOCALE */ 1289 1290 if (cc == ')') // format: (a,(c)) or (a,(c,d)) 1291 { 1292 q = quaternion<T>(a,b,c,d); 1293 } 1294 else // error 1295 { 1296 is.setstate(::std::ios_base::failbit); 1297 } 1298 } 1299 else // read "(a,b", possible: (a,b), (a,b,c), (a,b,c,d) 1300 { 1301 is.putback(ch); 1302 1303 is >> b; // we extract the second component 1304 1305 if (!is.good()) goto finish; 1306 1307 is >> ch; // get the fifth lexeme 1308 1309 if (!is.good()) goto finish; 1310 1311 #ifdef BOOST_NO_STD_LOCALE 1312 cc = ch; 1313 #else 1314 cc = ct.narrow(ch, char()); 1315 #endif /* BOOST_NO_STD_LOCALE */ 1316 1317 if (cc == ')') // format: (a,b) 1318 { 1319 q = quaternion<T>(a,b); 1320 } 1321 else if (cc == ',') // read "(a,b,", possible: (a,b,c), (a,b,c,d) 1322 { 1323 is >> c; // we extract the third component 1324 1325 if (!is.good()) goto finish; 1326 1327 is >> ch; // get the seventh lexeme 1328 1329 if (!is.good()) goto finish; 1330 1331 #ifdef BOOST_NO_STD_LOCALE 1332 cc = ch; 1333 #else 1334 cc = ct.narrow(ch, char()); 1335 #endif /* BOOST_NO_STD_LOCALE */ 1336 1337 if (cc == ')') // format: (a,b,c) 1338 { 1339 q = quaternion<T>(a,b,c); 1340 } 1341 else if (cc == ',') // read "(a,b,c,", possible: (a,b,c,d) 1342 { 1343 is >> d; // we extract the fourth component 1344 1345 if (!is.good()) goto finish; 1346 1347 is >> ch; // get the ninth lexeme 1348 1349 if (!is.good()) goto finish; 1350 1351 #ifdef BOOST_NO_STD_LOCALE 1352 cc = ch; 1353 #else 1354 cc = ct.narrow(ch, char()); 1355 #endif /* BOOST_NO_STD_LOCALE */ 1356 1357 if (cc == ')') // format: (a,b,c,d) 1358 { 1359 q = quaternion<T>(a,b,c,d); 1360 } 1361 else // error 1362 { 1363 is.setstate(::std::ios_base::failbit); 1364 } 1365 } 1366 else // error 1367 { 1368 is.setstate(::std::ios_base::failbit); 1369 } 1370 } 1371 else // error 1372 { 1373 is.setstate(::std::ios_base::failbit); 1374 } 1375 } 1376 } 1377 else // error 1378 { 1379 is.setstate(::std::ios_base::failbit); 1380 } 1381 } 1382 } 1383 else // format: a 1384 { 1385 is.putback(ch); 1386 1387 is >> a; // we extract the first component 1388 1389 if (!is.good()) goto finish; 1390 1391 q = quaternion<T>(a); 1392 } 1393 1394 finish: 1395 return(is); 1396 } 1397 1398 1399 template<typename T, typename charT, class traits> operator <<(::std::basic_ostream<charT,traits> & os,quaternion<T> const & q)1400 ::std::basic_ostream<charT,traits> & operator << ( ::std::basic_ostream<charT,traits> & os, 1401 quaternion<T> const & q) 1402 { 1403 ::std::basic_ostringstream<charT,traits> s; 1404 1405 s.flags(os.flags()); 1406 #ifdef BOOST_NO_STD_LOCALE 1407 #else 1408 s.imbue(os.getloc()); 1409 #endif /* BOOST_NO_STD_LOCALE */ 1410 s.precision(os.precision()); 1411 1412 s << '(' << q.R_component_1() << ',' 1413 << q.R_component_2() << ',' 1414 << q.R_component_3() << ',' 1415 << q.R_component_4() << ')'; 1416 1417 return os << s.str(); 1418 } 1419 1420 1421 // values 1422 1423 template<typename T> real(quaternion<T> const & q)1424 inline T real(quaternion<T> const & q) 1425 { 1426 return(q.real()); 1427 } 1428 1429 1430 template<typename T> unreal(quaternion<T> const & q)1431 inline quaternion<T> unreal(quaternion<T> const & q) 1432 { 1433 return(q.unreal()); 1434 } 1435 1436 1437 #define BOOST_QUATERNION_VALARRAY_LOADER \ 1438 using ::std::valarray; \ 1439 \ 1440 valarray<T> temp(4); \ 1441 \ 1442 temp[0] = q.R_component_1(); \ 1443 temp[1] = q.R_component_2(); \ 1444 temp[2] = q.R_component_3(); \ 1445 temp[3] = q.R_component_4(); 1446 1447 1448 template<typename T> sup(quaternion<T> const & q)1449 inline T sup(quaternion<T> const & q) 1450 { 1451 #ifdef BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP 1452 using ::std::abs; 1453 #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ 1454 1455 BOOST_QUATERNION_VALARRAY_LOADER 1456 1457 return((abs(temp).max)()); 1458 } 1459 1460 1461 template<typename T> l1(quaternion<T> const & q)1462 inline T l1(quaternion<T> const & q) 1463 { 1464 #ifdef BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP 1465 using ::std::abs; 1466 #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ 1467 1468 BOOST_QUATERNION_VALARRAY_LOADER 1469 1470 return(abs(temp).sum()); 1471 } 1472 1473 1474 template<typename T> abs(quaternion<T> const & q)1475 inline T abs(quaternion<T> const & q) 1476 { 1477 #ifdef BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP 1478 using ::std::abs; 1479 #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ 1480 1481 using ::std::sqrt; 1482 1483 BOOST_QUATERNION_VALARRAY_LOADER 1484 1485 T maxim = (abs(temp).max)(); // overflow protection 1486 1487 if (maxim == static_cast<T>(0)) 1488 { 1489 return(maxim); 1490 } 1491 else 1492 { 1493 T mixam = static_cast<T>(1)/maxim; // prefer multiplications over divisions 1494 1495 temp *= mixam; 1496 1497 temp *= temp; 1498 1499 return(maxim*sqrt(temp.sum())); 1500 } 1501 1502 //return(sqrt(norm(q))); 1503 } 1504 1505 1506 #undef BOOST_QUATERNION_VALARRAY_LOADER 1507 1508 1509 // Note: This is the Cayley norm, not the Euclidian norm... 1510 1511 template<typename T> norm(quaternion<T> const & q)1512 inline T norm(quaternion<T>const & q) 1513 { 1514 return(real(q*conj(q))); 1515 } 1516 1517 1518 template<typename T> conj(quaternion<T> const & q)1519 inline quaternion<T> conj(quaternion<T> const & q) 1520 { 1521 return(quaternion<T>( +q.R_component_1(), 1522 -q.R_component_2(), 1523 -q.R_component_3(), 1524 -q.R_component_4())); 1525 } 1526 1527 1528 template<typename T> spherical(T const & rho,T const & theta,T const & phi1,T const & phi2)1529 inline quaternion<T> spherical( T const & rho, 1530 T const & theta, 1531 T const & phi1, 1532 T const & phi2) 1533 { 1534 using ::std::cos; 1535 using ::std::sin; 1536 1537 //T a = cos(theta)*cos(phi1)*cos(phi2); 1538 //T b = sin(theta)*cos(phi1)*cos(phi2); 1539 //T c = sin(phi1)*cos(phi2); 1540 //T d = sin(phi2); 1541 1542 T courrant = static_cast<T>(1); 1543 1544 T d = sin(phi2); 1545 1546 courrant *= cos(phi2); 1547 1548 T c = sin(phi1)*courrant; 1549 1550 courrant *= cos(phi1); 1551 1552 T b = sin(theta)*courrant; 1553 T a = cos(theta)*courrant; 1554 1555 return(rho*quaternion<T>(a,b,c,d)); 1556 } 1557 1558 1559 template<typename T> semipolar(T const & rho,T const & alpha,T const & theta1,T const & theta2)1560 inline quaternion<T> semipolar( T const & rho, 1561 T const & alpha, 1562 T const & theta1, 1563 T const & theta2) 1564 { 1565 using ::std::cos; 1566 using ::std::sin; 1567 1568 T a = cos(alpha)*cos(theta1); 1569 T b = cos(alpha)*sin(theta1); 1570 T c = sin(alpha)*cos(theta2); 1571 T d = sin(alpha)*sin(theta2); 1572 1573 return(rho*quaternion<T>(a,b,c,d)); 1574 } 1575 1576 1577 template<typename T> multipolar(T const & rho1,T const & theta1,T const & rho2,T const & theta2)1578 inline quaternion<T> multipolar( T const & rho1, 1579 T const & theta1, 1580 T const & rho2, 1581 T const & theta2) 1582 { 1583 using ::std::cos; 1584 using ::std::sin; 1585 1586 T a = rho1*cos(theta1); 1587 T b = rho1*sin(theta1); 1588 T c = rho2*cos(theta2); 1589 T d = rho2*sin(theta2); 1590 1591 return(quaternion<T>(a,b,c,d)); 1592 } 1593 1594 1595 template<typename T> cylindrospherical(T const & t,T const & radius,T const & longitude,T const & latitude)1596 inline quaternion<T> cylindrospherical( T const & t, 1597 T const & radius, 1598 T const & longitude, 1599 T const & latitude) 1600 { 1601 using ::std::cos; 1602 using ::std::sin; 1603 1604 1605 1606 T b = radius*cos(longitude)*cos(latitude); 1607 T c = radius*sin(longitude)*cos(latitude); 1608 T d = radius*sin(latitude); 1609 1610 return(quaternion<T>(t,b,c,d)); 1611 } 1612 1613 1614 template<typename T> cylindrical(T const & r,T const & angle,T const & h1,T const & h2)1615 inline quaternion<T> cylindrical(T const & r, 1616 T const & angle, 1617 T const & h1, 1618 T const & h2) 1619 { 1620 using ::std::cos; 1621 using ::std::sin; 1622 1623 T a = r*cos(angle); 1624 T b = r*sin(angle); 1625 1626 return(quaternion<T>(a,b,h1,h2)); 1627 } 1628 1629 1630 // transcendentals 1631 // (please see the documentation) 1632 1633 1634 template<typename T> exp(quaternion<T> const & q)1635 inline quaternion<T> exp(quaternion<T> const & q) 1636 { 1637 using ::std::exp; 1638 using ::std::cos; 1639 1640 using ::boost::math::sinc_pi; 1641 1642 T u = exp(real(q)); 1643 1644 T z = abs(unreal(q)); 1645 1646 T w = sinc_pi(z); 1647 1648 return(u*quaternion<T>(cos(z), 1649 w*q.R_component_2(), w*q.R_component_3(), 1650 w*q.R_component_4())); 1651 } 1652 1653 1654 template<typename T> cos(quaternion<T> const & q)1655 inline quaternion<T> cos(quaternion<T> const & q) 1656 { 1657 using ::std::sin; 1658 using ::std::cos; 1659 using ::std::cosh; 1660 1661 using ::boost::math::sinhc_pi; 1662 1663 T z = abs(unreal(q)); 1664 1665 T w = -sin(q.real())*sinhc_pi(z); 1666 1667 return(quaternion<T>(cos(q.real())*cosh(z), 1668 w*q.R_component_2(), w*q.R_component_3(), 1669 w*q.R_component_4())); 1670 } 1671 1672 1673 template<typename T> sin(quaternion<T> const & q)1674 inline quaternion<T> sin(quaternion<T> const & q) 1675 { 1676 using ::std::sin; 1677 using ::std::cos; 1678 using ::std::cosh; 1679 1680 using ::boost::math::sinhc_pi; 1681 1682 T z = abs(unreal(q)); 1683 1684 T w = +cos(q.real())*sinhc_pi(z); 1685 1686 return(quaternion<T>(sin(q.real())*cosh(z), 1687 w*q.R_component_2(), w*q.R_component_3(), 1688 w*q.R_component_4())); 1689 } 1690 1691 1692 template<typename T> tan(quaternion<T> const & q)1693 inline quaternion<T> tan(quaternion<T> const & q) 1694 { 1695 return(sin(q)/cos(q)); 1696 } 1697 1698 1699 template<typename T> cosh(quaternion<T> const & q)1700 inline quaternion<T> cosh(quaternion<T> const & q) 1701 { 1702 return((exp(+q)+exp(-q))/static_cast<T>(2)); 1703 } 1704 1705 1706 template<typename T> sinh(quaternion<T> const & q)1707 inline quaternion<T> sinh(quaternion<T> const & q) 1708 { 1709 return((exp(+q)-exp(-q))/static_cast<T>(2)); 1710 } 1711 1712 1713 template<typename T> tanh(quaternion<T> const & q)1714 inline quaternion<T> tanh(quaternion<T> const & q) 1715 { 1716 return(sinh(q)/cosh(q)); 1717 } 1718 1719 1720 template<typename T> pow(quaternion<T> const & q,int n)1721 quaternion<T> pow(quaternion<T> const & q, 1722 int n) 1723 { 1724 if (n > 1) 1725 { 1726 int m = n>>1; 1727 1728 quaternion<T> result = pow(q, m); 1729 1730 result *= result; 1731 1732 if (n != (m<<1)) 1733 { 1734 result *= q; // n odd 1735 } 1736 1737 return(result); 1738 } 1739 else if (n == 1) 1740 { 1741 return(q); 1742 } 1743 else if (n == 0) 1744 { 1745 return(quaternion<T>(static_cast<T>(1))); 1746 } 1747 else /* n < 0 */ 1748 { 1749 return(pow(quaternion<T>(static_cast<T>(1))/q,-n)); 1750 } 1751 } 1752 1753 1754 // helper templates for converting copy constructors (definition) 1755 1756 namespace detail 1757 { 1758 1759 template< typename T, 1760 typename U 1761 > quaternion_type_converter(quaternion<U> const & rhs)1762 quaternion<T> quaternion_type_converter(quaternion<U> const & rhs) 1763 { 1764 return(quaternion<T>( static_cast<T>(rhs.R_component_1()), 1765 static_cast<T>(rhs.R_component_2()), 1766 static_cast<T>(rhs.R_component_3()), 1767 static_cast<T>(rhs.R_component_4()))); 1768 } 1769 } 1770 } 1771 } 1772 1773 #endif /* BOOST_QUATERNION_HPP */ 1774