1 // Boost.Units - A C++ library for zero-overhead dimensional analysis and 2 // unit/quantity manipulation and conversion 3 // 4 // Copyright (C) 2003-2008 Matthias Christian Schabel 5 // Copyright (C) 2008 Steven Watanabe 6 // 7 // Distributed under the Boost Software License, Version 1.0. (See 8 // accompanying file LICENSE_1_0.txt or copy at 9 // http://www.boost.org/LICENSE_1_0.txt) 10 11 #ifndef BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP 12 #define BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP 13 14 #include <boost/units/static_rational.hpp> 15 #include <boost/mpl/next.hpp> 16 #include <boost/mpl/arithmetic.hpp> 17 #include <boost/mpl/and.hpp> 18 #include <boost/mpl/assert.hpp> 19 20 #include <boost/units/dim.hpp> 21 #include <boost/units/dimensionless_type.hpp> 22 #include <boost/units/static_rational.hpp> 23 #include <boost/units/detail/dimension_list.hpp> 24 #include <boost/units/detail/sort.hpp> 25 26 namespace boost { 27 28 namespace units { 29 30 namespace detail { 31 32 // typedef list<rational> equation; 33 34 template<int N> 35 struct eliminate_from_pair_of_equations_impl; 36 37 template<class E1, class E2> 38 struct eliminate_from_pair_of_equations; 39 40 template<int N> 41 struct elimination_impl; 42 43 template<bool is_zero, bool element_is_last> 44 struct elimination_skip_leading_zeros_impl; 45 46 template<class Equation, class Vars> 47 struct substitute; 48 49 template<int N> 50 struct substitute_impl; 51 52 template<bool is_end> 53 struct solve_impl; 54 55 template<class T> 56 struct solve; 57 58 template<int N> 59 struct check_extra_equations_impl; 60 61 template<int N> 62 struct normalize_units_impl; 63 64 struct inconsistent {}; 65 66 // generally useful utilies. 67 68 template<int N> 69 struct divide_equation { 70 template<class Begin, class Divisor> 71 struct apply { 72 typedef list<typename mpl::divides<typename Begin::item, Divisor>::type, typename divide_equation<N - 1>::template apply<typename Begin::next, Divisor>::type> type; 73 }; 74 }; 75 76 template<> 77 struct divide_equation<0> { 78 template<class Begin, class Divisor> 79 struct apply { 80 typedef dimensionless_type type; 81 }; 82 }; 83 84 // eliminate_from_pair_of_equations takes a pair of 85 // equations and eliminates the first variable. 86 // 87 // equation eliminate_from_pair_of_equations(equation l1, equation l2) { 88 // rational x1 = l1.front(); 89 // rational x2 = l2.front(); 90 // return(transform(pop_front(l1), pop_front(l2), _1 * x2 - _2 * x1)); 91 // } 92 93 template<int N> 94 struct eliminate_from_pair_of_equations_impl { 95 template<class Begin1, class Begin2, class X1, class X2> 96 struct apply { 97 typedef list< 98 typename mpl::minus< 99 typename mpl::times<typename Begin1::item, X2>::type, 100 typename mpl::times<typename Begin2::item, X1>::type 101 >::type, 102 typename eliminate_from_pair_of_equations_impl<N - 1>::template apply< 103 typename Begin1::next, 104 typename Begin2::next, 105 X1, 106 X2 107 >::type 108 > type; 109 }; 110 }; 111 112 template<> 113 struct eliminate_from_pair_of_equations_impl<0> { 114 template<class Begin1, class Begin2, class X1, class X2> 115 struct apply { 116 typedef dimensionless_type type; 117 }; 118 }; 119 120 template<class E1, class E2> 121 struct eliminate_from_pair_of_equations { 122 typedef E1 begin1; 123 typedef E2 begin2; 124 typedef typename eliminate_from_pair_of_equations_impl<(E1::size::value - 1)>::template apply< 125 typename begin1::next, 126 typename begin2::next, 127 typename begin1::item, 128 typename begin2::item 129 >::type type; 130 }; 131 132 133 134 // Stage 1. Determine which dimensions should 135 // have dummy base units. For this purpose 136 // row reduce the matrix. 137 138 template<int N> 139 struct make_zero_vector { 140 typedef list<static_rational<0>, typename make_zero_vector<N - 1>::type> type; 141 }; 142 template<> 143 struct make_zero_vector<0> { 144 typedef dimensionless_type type; 145 }; 146 147 template<int Column, int TotalColumns> 148 struct create_row_of_identity { 149 typedef list<static_rational<0>, typename create_row_of_identity<Column - 1, TotalColumns - 1>::type> type; 150 }; 151 template<int TotalColumns> 152 struct create_row_of_identity<0, TotalColumns> { 153 typedef list<static_rational<1>, typename make_zero_vector<TotalColumns - 1>::type> type; 154 }; 155 template<int Column> 156 struct create_row_of_identity<Column, 0> { 157 // error 158 }; 159 160 template<int RemainingRows> 161 struct determine_extra_equations_impl; 162 163 template<bool first_is_zero, bool is_last> 164 struct determine_extra_equations_skip_zeros_impl; 165 166 // not the last row and not zero. 167 template<> 168 struct determine_extra_equations_skip_zeros_impl<false, false> { 169 template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result> 170 struct apply { 171 // remove the equation being eliminated against from the set of equations. 172 typedef typename determine_extra_equations_impl<RemainingRows - 1>::template apply<typename RowsBegin::next, typename RowsBegin::item>::type next_equations; 173 // since this column was present, strip it out. 174 typedef Result type; 175 }; 176 }; 177 178 // the last row but not zero. 179 template<> 180 struct determine_extra_equations_skip_zeros_impl<false, true> { 181 template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result> 182 struct apply { 183 // remove this equation. 184 typedef dimensionless_type next_equations; 185 // since this column was present, strip it out. 186 typedef Result type; 187 }; 188 }; 189 190 191 // the first columns is zero but it is not the last column. 192 // continue with the same loop. 193 template<> 194 struct determine_extra_equations_skip_zeros_impl<true, false> { 195 template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result> 196 struct apply { 197 typedef typename RowsBegin::next::item next_row; 198 typedef typename determine_extra_equations_skip_zeros_impl< 199 next_row::item::Numerator == 0, 200 RemainingRows == 2 // the next one will be the last. 201 >::template apply< 202 typename RowsBegin::next, 203 RemainingRows - 1, 204 CurrentColumn, 205 TotalColumns, 206 Result 207 > next; 208 typedef list<typename RowsBegin::item::next, typename next::next_equations> next_equations; 209 typedef typename next::type type; 210 }; 211 }; 212 213 // all the elements in this column are zero. 214 template<> 215 struct determine_extra_equations_skip_zeros_impl<true, true> { 216 template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result> 217 struct apply { 218 typedef list<typename RowsBegin::item::next, dimensionless_type> next_equations; 219 typedef list<typename create_row_of_identity<CurrentColumn, TotalColumns>::type, Result> type; 220 }; 221 }; 222 223 template<int RemainingRows> 224 struct determine_extra_equations_impl { 225 template<class RowsBegin, class EliminateAgainst> 226 struct apply { 227 typedef list< 228 typename eliminate_from_pair_of_equations<typename RowsBegin::item, EliminateAgainst>::type, 229 typename determine_extra_equations_impl<RemainingRows-1>::template apply<typename RowsBegin::next, EliminateAgainst>::type 230 > type; 231 }; 232 }; 233 234 template<> 235 struct determine_extra_equations_impl<0> { 236 template<class RowsBegin, class EliminateAgainst> 237 struct apply { 238 typedef dimensionless_type type; 239 }; 240 }; 241 242 template<int RemainingColumns, bool is_done> 243 struct determine_extra_equations { 244 template<class RowsBegin, int TotalColumns, class Result> 245 struct apply { 246 typedef typename RowsBegin::item top_row; 247 typedef typename determine_extra_equations_skip_zeros_impl< 248 top_row::item::Numerator == 0, 249 RowsBegin::size::value == 1 250 >::template apply< 251 RowsBegin, 252 RowsBegin::size::value, 253 TotalColumns - RemainingColumns, 254 TotalColumns, 255 Result 256 > column_info; 257 typedef typename determine_extra_equations< 258 RemainingColumns - 1, 259 column_info::next_equations::size::value == 0 260 >::template apply< 261 typename column_info::next_equations, 262 TotalColumns, 263 typename column_info::type 264 >::type type; 265 }; 266 }; 267 268 template<int RemainingColumns> 269 struct determine_extra_equations<RemainingColumns, true> { 270 template<class RowsBegin, int TotalColumns, class Result> 271 struct apply { 272 typedef typename determine_extra_equations<RemainingColumns - 1, true>::template apply< 273 RowsBegin, 274 TotalColumns, 275 list<typename create_row_of_identity<TotalColumns - RemainingColumns, TotalColumns>::type, Result> 276 >::type type; 277 }; 278 }; 279 280 template<> 281 struct determine_extra_equations<0, true> { 282 template<class RowsBegin, int TotalColumns, class Result> 283 struct apply { 284 typedef Result type; 285 }; 286 }; 287 288 // Stage 2 289 // invert the matrix using Gauss-Jordan elimination 290 291 292 template<bool is_zero, bool is_last> 293 struct invert_strip_leading_zeroes; 294 295 template<int N> 296 struct invert_handle_after_pivot_row; 297 298 // When processing column N, none of the first N rows 299 // can be the pivot column. 300 template<int N> 301 struct invert_handle_inital_rows { 302 template<class RowsBegin, class IdentityBegin> 303 struct apply { 304 typedef typename invert_handle_inital_rows<N - 1>::template apply< 305 typename RowsBegin::next, 306 typename IdentityBegin::next 307 > next; 308 typedef typename RowsBegin::item current_row; 309 typedef typename IdentityBegin::item current_identity_row; 310 typedef typename next::pivot_row pivot_row; 311 typedef typename next::identity_pivot_row identity_pivot_row; 312 typedef list< 313 typename eliminate_from_pair_of_equations_impl<(current_row::size::value) - 1>::template apply< 314 typename current_row::next, 315 pivot_row, 316 typename current_row::item, 317 static_rational<1> 318 >::type, 319 typename next::new_matrix 320 > new_matrix; 321 typedef list< 322 typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply< 323 current_identity_row, 324 identity_pivot_row, 325 typename current_row::item, 326 static_rational<1> 327 >::type, 328 typename next::identity_result 329 > identity_result; 330 }; 331 }; 332 333 // This handles the switch to searching for a pivot column. 334 // The pivot row will be propagated up in the typedefs 335 // pivot_row and identity_pivot_row. It is inserted here. 336 template<> 337 struct invert_handle_inital_rows<0> { 338 template<class RowsBegin, class IdentityBegin> 339 struct apply { 340 typedef typename RowsBegin::item current_row; 341 typedef typename invert_strip_leading_zeroes< 342 (current_row::item::Numerator == 0), 343 (RowsBegin::size::value == 1) 344 >::template apply< 345 RowsBegin, 346 IdentityBegin 347 > next; 348 // results 349 typedef list<typename next::pivot_row, typename next::new_matrix> new_matrix; 350 typedef list<typename next::identity_pivot_row, typename next::identity_result> identity_result; 351 typedef typename next::pivot_row pivot_row; 352 typedef typename next::identity_pivot_row identity_pivot_row; 353 }; 354 }; 355 356 // The first internal element which is not zero. 357 template<> 358 struct invert_strip_leading_zeroes<false, false> { 359 template<class RowsBegin, class IdentityBegin> 360 struct apply { 361 typedef typename RowsBegin::item current_row; 362 typedef typename current_row::item current_value; 363 typedef typename divide_equation<(current_row::size::value - 1)>::template apply<typename current_row::next, current_value>::type new_equation; 364 typedef typename divide_equation<(IdentityBegin::item::size::value)>::template apply<typename IdentityBegin::item, current_value>::type transformed_identity_equation; 365 typedef typename invert_handle_after_pivot_row<(RowsBegin::size::value - 1)>::template apply< 366 typename RowsBegin::next, 367 typename IdentityBegin::next, 368 new_equation, 369 transformed_identity_equation 370 > next; 371 372 // results 373 // Note that we don't add the pivot row to the 374 // results here, because it needs to propagated up 375 // to the diagonal. 376 typedef typename next::new_matrix new_matrix; 377 typedef typename next::identity_result identity_result; 378 typedef new_equation pivot_row; 379 typedef transformed_identity_equation identity_pivot_row; 380 }; 381 }; 382 383 // The one and only non-zero element--at the end 384 template<> 385 struct invert_strip_leading_zeroes<false, true> { 386 template<class RowsBegin, class IdentityBegin> 387 struct apply { 388 typedef typename RowsBegin::item current_row; 389 typedef typename current_row::item current_value; 390 typedef typename divide_equation<(current_row::size::value - 1)>::template apply<typename current_row::next, current_value>::type new_equation; 391 typedef typename divide_equation<(IdentityBegin::item::size::value)>::template apply<typename IdentityBegin::item, current_value>::type transformed_identity_equation; 392 393 // results 394 // Note that we don't add the pivot row to the 395 // results here, because it needs to propagated up 396 // to the diagonal. 397 typedef dimensionless_type identity_result; 398 typedef dimensionless_type new_matrix; 399 typedef new_equation pivot_row; 400 typedef transformed_identity_equation identity_pivot_row; 401 }; 402 }; 403 404 // One of the initial zeroes 405 template<> 406 struct invert_strip_leading_zeroes<true, false> { 407 template<class RowsBegin, class IdentityBegin> 408 struct apply { 409 typedef typename RowsBegin::item current_row; 410 typedef typename RowsBegin::next::item next_row; 411 typedef typename invert_strip_leading_zeroes< 412 next_row::item::Numerator == 0, 413 RowsBegin::size::value == 2 414 >::template apply< 415 typename RowsBegin::next, 416 typename IdentityBegin::next 417 > next; 418 typedef typename IdentityBegin::item current_identity_row; 419 // these are propagated up. 420 typedef typename next::pivot_row pivot_row; 421 typedef typename next::identity_pivot_row identity_pivot_row; 422 typedef list< 423 typename eliminate_from_pair_of_equations_impl<(current_row::size::value - 1)>::template apply< 424 typename current_row::next, 425 pivot_row, 426 typename current_row::item, 427 static_rational<1> 428 >::type, 429 typename next::new_matrix 430 > new_matrix; 431 typedef list< 432 typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply< 433 current_identity_row, 434 identity_pivot_row, 435 typename current_row::item, 436 static_rational<1> 437 >::type, 438 typename next::identity_result 439 > identity_result; 440 }; 441 }; 442 443 // the last element, and is zero. 444 // Should never happen. 445 template<> 446 struct invert_strip_leading_zeroes<true, true> { 447 }; 448 449 template<int N> 450 struct invert_handle_after_pivot_row { 451 template<class RowsBegin, class IdentityBegin, class MatrixPivot, class IdentityPivot> 452 struct apply { 453 typedef typename invert_handle_after_pivot_row<N - 1>::template apply< 454 typename RowsBegin::next, 455 typename IdentityBegin::next, 456 MatrixPivot, 457 IdentityPivot 458 > next; 459 typedef typename RowsBegin::item current_row; 460 typedef typename IdentityBegin::item current_identity_row; 461 typedef MatrixPivot pivot_row; 462 typedef IdentityPivot identity_pivot_row; 463 464 // results 465 typedef list< 466 typename eliminate_from_pair_of_equations_impl<(current_row::size::value - 1)>::template apply< 467 typename current_row::next, 468 pivot_row, 469 typename current_row::item, 470 static_rational<1> 471 >::type, 472 typename next::new_matrix 473 > new_matrix; 474 typedef list< 475 typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply< 476 current_identity_row, 477 identity_pivot_row, 478 typename current_row::item, 479 static_rational<1> 480 >::type, 481 typename next::identity_result 482 > identity_result; 483 }; 484 }; 485 486 template<> 487 struct invert_handle_after_pivot_row<0> { 488 template<class RowsBegin, class IdentityBegin, class MatrixPivot, class IdentityPivot> 489 struct apply { 490 typedef dimensionless_type new_matrix; 491 typedef dimensionless_type identity_result; 492 }; 493 }; 494 495 template<int N> 496 struct invert_impl { 497 template<class RowsBegin, class IdentityBegin> 498 struct apply { 499 typedef typename invert_handle_inital_rows<RowsBegin::size::value - N>::template apply<RowsBegin, IdentityBegin> process_column; 500 typedef typename invert_impl<N - 1>::template apply< 501 typename process_column::new_matrix, 502 typename process_column::identity_result 503 >::type type; 504 }; 505 }; 506 507 template<> 508 struct invert_impl<0> { 509 template<class RowsBegin, class IdentityBegin> 510 struct apply { 511 typedef IdentityBegin type; 512 }; 513 }; 514 515 template<int N> 516 struct make_identity { 517 template<int Size> 518 struct apply { 519 typedef list<typename create_row_of_identity<Size - N, Size>::type, typename make_identity<N - 1>::template apply<Size>::type> type; 520 }; 521 }; 522 523 template<> 524 struct make_identity<0> { 525 template<int Size> 526 struct apply { 527 typedef dimensionless_type type; 528 }; 529 }; 530 531 template<class Matrix> 532 struct make_square_and_invert { 533 typedef typename Matrix::item top_row; 534 typedef typename determine_extra_equations<(top_row::size::value), false>::template apply< 535 Matrix, // RowsBegin 536 top_row::size::value, // TotalColumns 537 Matrix // Result 538 >::type invertible; 539 typedef typename invert_impl<invertible::size::value>::template apply< 540 invertible, 541 typename make_identity<invertible::size::value>::template apply<invertible::size::value>::type 542 >::type type; 543 }; 544 545 546 // find_base_dimensions takes a list of 547 // base_units and returns a sorted list 548 // of all the base_dimensions they use. 549 // 550 // list<base_dimension> find_base_dimensions(list<base_unit> l) { 551 // set<base_dimension> dimensions; 552 // for_each(base_unit unit : l) { 553 // for_each(dim d : unit.dimension_type) { 554 // dimensions = insert(dimensions, d.tag_type); 555 // } 556 // } 557 // return(sort(dimensions, _1 > _2, front_inserter(list<base_dimension>()))); 558 // } 559 560 typedef char set_no; 561 struct set_yes { set_no dummy[2]; }; 562 563 template<class T> 564 struct wrap {}; 565 566 struct set_end { 567 static set_no lookup(...); 568 typedef mpl::long_<0> size; 569 }; 570 571 template<class T, class Next> 572 struct set : Next { 573 using Next::lookup; 574 static set_yes lookup(wrap<T>*); 575 typedef T item; 576 typedef Next next; 577 typedef typename mpl::next<typename Next::size>::type size; 578 }; 579 580 template<bool has_key> 581 struct set_insert; 582 583 template<> 584 struct set_insert<true> { 585 template<class Set, class T> 586 struct apply { 587 typedef Set type; 588 }; 589 }; 590 591 template<> 592 struct set_insert<false> { 593 template<class Set, class T> 594 struct apply { 595 typedef set<T, Set> type; 596 }; 597 }; 598 599 template<class Set, class T> 600 struct has_key { 601 BOOST_STATIC_CONSTEXPR long size = sizeof(Set::lookup((wrap<T>*)0)); 602 BOOST_STATIC_CONSTEXPR bool value = (size == sizeof(set_yes)); 603 }; 604 605 template<int N> 606 struct find_base_dimensions_impl_impl { 607 template<class Begin, class S> 608 struct apply { 609 typedef typename find_base_dimensions_impl_impl<N-1>::template apply< 610 typename Begin::next, 611 S 612 >::type next; 613 614 typedef typename set_insert< 615 (has_key<next, typename Begin::item::tag_type>::value) 616 >::template apply< 617 next, 618 typename Begin::item::tag_type 619 >::type type; 620 }; 621 }; 622 623 template<> 624 struct find_base_dimensions_impl_impl<0> { 625 template<class Begin, class S> 626 struct apply { 627 typedef S type; 628 }; 629 }; 630 631 template<int N> 632 struct find_base_dimensions_impl { 633 template<class Begin> 634 struct apply { 635 typedef typename find_base_dimensions_impl_impl<(Begin::item::dimension_type::size::value)>::template apply< 636 typename Begin::item::dimension_type, 637 typename find_base_dimensions_impl<N-1>::template apply<typename Begin::next>::type 638 >::type type; 639 }; 640 }; 641 642 template<> 643 struct find_base_dimensions_impl<0> { 644 template<class Begin> 645 struct apply { 646 typedef set_end type; 647 }; 648 }; 649 650 template<class T> 651 struct find_base_dimensions { 652 typedef typename insertion_sort< 653 typename find_base_dimensions_impl< 654 (T::size::value) 655 >::template apply<T>::type 656 >::type type; 657 }; 658 659 // calculate_base_dimension_coefficients finds 660 // the coefficients corresponding to the first 661 // base_dimension in each of the dimension_lists. 662 // It returns two values. The first result 663 // is a list of the coefficients. The second 664 // is a list with all the incremented iterators. 665 // When we encounter a base_dimension that is 666 // missing from a dimension_list, we do not 667 // increment the iterator and we set the 668 // coefficient to zero. 669 670 template<bool has_dimension> 671 struct calculate_base_dimension_coefficients_func; 672 673 template<> 674 struct calculate_base_dimension_coefficients_func<true> { 675 template<class T> 676 struct apply { 677 typedef typename T::item::value_type type; 678 typedef typename T::next next; 679 }; 680 }; 681 682 template<> 683 struct calculate_base_dimension_coefficients_func<false> { 684 template<class T> 685 struct apply { 686 typedef static_rational<0> type; 687 typedef T next; 688 }; 689 }; 690 691 // begins_with_dimension returns true iff its first 692 // parameter is a valid iterator which yields its 693 // second parameter when dereferenced. 694 695 template<class Iterator> 696 struct begins_with_dimension { 697 template<class Dim> 698 struct apply : 699 boost::is_same< 700 Dim, 701 typename Iterator::item::tag_type 702 > {}; 703 }; 704 705 template<> 706 struct begins_with_dimension<dimensionless_type> { 707 template<class Dim> 708 struct apply : mpl::false_ {}; 709 }; 710 711 template<int N> 712 struct calculate_base_dimension_coefficients_impl { 713 template<class BaseUnitDimensions,class Dim,class T> 714 struct apply { 715 typedef typename calculate_base_dimension_coefficients_func< 716 begins_with_dimension<typename BaseUnitDimensions::item>::template apply< 717 Dim 718 >::value 719 >::template apply< 720 typename BaseUnitDimensions::item 721 > result; 722 typedef typename calculate_base_dimension_coefficients_impl<N-1>::template apply< 723 typename BaseUnitDimensions::next, 724 Dim, 725 list<typename result::type, T> 726 > next_; 727 typedef typename next_::type type; 728 typedef list<typename result::next, typename next_::next> next; 729 }; 730 }; 731 732 template<> 733 struct calculate_base_dimension_coefficients_impl<0> { 734 template<class Begin, class BaseUnitDimensions, class T> 735 struct apply { 736 typedef T type; 737 typedef dimensionless_type next; 738 }; 739 }; 740 741 // add_zeroes pushs N zeroes onto the 742 // front of a list. 743 // 744 // list<rational> add_zeroes(list<rational> l, int N) { 745 // if(N == 0) { 746 // return(l); 747 // } else { 748 // return(push_front(add_zeroes(l, N-1), 0)); 749 // } 750 // } 751 752 template<int N> 753 struct add_zeroes_impl { 754 // If you get an error here and your base units are 755 // in fact linearly independent, please report it. 756 BOOST_MPL_ASSERT_MSG((N > 0), base_units_are_probably_not_linearly_independent, (void)); 757 template<class T> 758 struct apply { 759 typedef list< 760 static_rational<0>, 761 typename add_zeroes_impl<N-1>::template apply<T>::type 762 > type; 763 }; 764 }; 765 766 template<> 767 struct add_zeroes_impl<0> { 768 template<class T> 769 struct apply { 770 typedef T type; 771 }; 772 }; 773 774 // expand_dimensions finds the exponents of 775 // a set of dimensions in a dimension_list. 776 // the second parameter is assumed to be 777 // a superset of the base_dimensions of 778 // the first parameter. 779 // 780 // list<rational> expand_dimensions(dimension_list, list<base_dimension>); 781 782 template<int N> 783 struct expand_dimensions { 784 template<class Begin, class DimensionIterator> 785 struct apply { 786 typedef typename calculate_base_dimension_coefficients_func< 787 begins_with_dimension<DimensionIterator>::template apply<typename Begin::item>::value 788 >::template apply<DimensionIterator> result; 789 typedef list< 790 typename result::type, 791 typename expand_dimensions<N-1>::template apply<typename Begin::next, typename result::next>::type 792 > type; 793 }; 794 }; 795 796 template<> 797 struct expand_dimensions<0> { 798 template<class Begin, class DimensionIterator> 799 struct apply { 800 typedef dimensionless_type type; 801 }; 802 }; 803 804 template<int N> 805 struct create_unit_matrix { 806 template<class Begin, class Dimensions> 807 struct apply { 808 typedef typename create_unit_matrix<N - 1>::template apply<typename Begin::next, Dimensions>::type next; 809 typedef list<typename expand_dimensions<Dimensions::size::value>::template apply<Dimensions, typename Begin::item::dimension_type>::type, next> type; 810 }; 811 }; 812 813 template<> 814 struct create_unit_matrix<0> { 815 template<class Begin, class Dimensions> 816 struct apply { 817 typedef dimensionless_type type; 818 }; 819 }; 820 821 template<class T> 822 struct normalize_units { 823 typedef typename find_base_dimensions<T>::type dimensions; 824 typedef typename create_unit_matrix<(T::size::value)>::template apply< 825 T, 826 dimensions 827 >::type matrix; 828 typedef typename make_square_and_invert<matrix>::type type; 829 BOOST_STATIC_CONSTEXPR long extra = (type::size::value) - (T::size::value); 830 }; 831 832 // multiply_add_units computes M x V 833 // where M is a matrix and V is a horizontal 834 // vector 835 // 836 // list<rational> multiply_add_units(list<list<rational> >, list<rational>); 837 838 template<int N> 839 struct multiply_add_units_impl { 840 template<class Begin1, class Begin2 ,class X> 841 struct apply { 842 typedef list< 843 typename mpl::plus< 844 typename mpl::times< 845 typename Begin2::item, 846 X 847 >::type, 848 typename Begin1::item 849 >::type, 850 typename multiply_add_units_impl<N-1>::template apply< 851 typename Begin1::next, 852 typename Begin2::next, 853 X 854 >::type 855 > type; 856 }; 857 }; 858 859 template<> 860 struct multiply_add_units_impl<0> { 861 template<class Begin1, class Begin2 ,class X> 862 struct apply { 863 typedef dimensionless_type type; 864 }; 865 }; 866 867 template<int N> 868 struct multiply_add_units { 869 template<class Begin1, class Begin2> 870 struct apply { 871 typedef typename multiply_add_units_impl< 872 (Begin2::item::size::value) 873 >::template apply< 874 typename multiply_add_units<N-1>::template apply< 875 typename Begin1::next, 876 typename Begin2::next 877 >::type, 878 typename Begin2::item, 879 typename Begin1::item 880 >::type type; 881 }; 882 }; 883 884 template<> 885 struct multiply_add_units<1> { 886 template<class Begin1, class Begin2> 887 struct apply { 888 typedef typename add_zeroes_impl< 889 (Begin2::item::size::value) 890 >::template apply<dimensionless_type>::type type1; 891 typedef typename multiply_add_units_impl< 892 (Begin2::item::size::value) 893 >::template apply< 894 type1, 895 typename Begin2::item, 896 typename Begin1::item 897 >::type type; 898 }; 899 }; 900 901 902 // strip_zeroes erases the first N elements of a list if 903 // they are all zero, otherwise returns inconsistent 904 // 905 // list strip_zeroes(list l, int N) { 906 // if(N == 0) { 907 // return(l); 908 // } else if(l.front == 0) { 909 // return(strip_zeroes(pop_front(l), N-1)); 910 // } else { 911 // return(inconsistent); 912 // } 913 // } 914 915 template<int N> 916 struct strip_zeroes_impl; 917 918 template<class T> 919 struct strip_zeroes_func { 920 template<class L, int N> 921 struct apply { 922 typedef inconsistent type; 923 }; 924 }; 925 926 template<> 927 struct strip_zeroes_func<static_rational<0> > { 928 template<class L, int N> 929 struct apply { 930 typedef typename strip_zeroes_impl<N-1>::template apply<typename L::next>::type type; 931 }; 932 }; 933 934 template<int N> 935 struct strip_zeroes_impl { 936 template<class T> 937 struct apply { 938 typedef typename strip_zeroes_func<typename T::item>::template apply<T, N>::type type; 939 }; 940 }; 941 942 template<> 943 struct strip_zeroes_impl<0> { 944 template<class T> 945 struct apply { 946 typedef T type; 947 }; 948 }; 949 950 // Given a list of base_units, computes the 951 // exponents of each base unit for a given 952 // dimension. 953 // 954 // list<rational> calculate_base_unit_exponents(list<base_unit> units, dimension_list dimensions); 955 956 template<class T> 957 struct is_base_dimension_unit { 958 typedef mpl::false_ type; 959 typedef void base_dimension_type; 960 }; 961 template<class T> 962 struct is_base_dimension_unit<list<dim<T, static_rational<1> >, dimensionless_type> > { 963 typedef mpl::true_ type; 964 typedef T base_dimension_type; 965 }; 966 967 template<int N> 968 struct is_simple_system_impl { 969 template<class Begin, class Prev> 970 struct apply { 971 typedef is_base_dimension_unit<typename Begin::item::dimension_type> test; 972 typedef mpl::and_< 973 typename test::type, 974 mpl::less<Prev, typename test::base_dimension_type>, 975 typename is_simple_system_impl<N-1>::template apply< 976 typename Begin::next, 977 typename test::base_dimension_type 978 > 979 > type; 980 BOOST_STATIC_CONSTEXPR bool value = (type::value); 981 }; 982 }; 983 984 template<> 985 struct is_simple_system_impl<0> { 986 template<class Begin, class Prev> 987 struct apply : mpl::true_ { 988 }; 989 }; 990 991 template<class T> 992 struct is_simple_system { 993 typedef T Begin; 994 typedef is_base_dimension_unit<typename Begin::item::dimension_type> test; 995 typedef typename mpl::and_< 996 typename test::type, 997 typename is_simple_system_impl< 998 T::size::value - 1 999 >::template apply< 1000 typename Begin::next::type, 1001 typename test::base_dimension_type 1002 > 1003 >::type type; 1004 BOOST_STATIC_CONSTEXPR bool value = type::value; 1005 }; 1006 1007 template<bool> 1008 struct calculate_base_unit_exponents_impl; 1009 1010 template<> 1011 struct calculate_base_unit_exponents_impl<true> { 1012 template<class T, class Dimensions> 1013 struct apply { 1014 typedef typename expand_dimensions<(T::size::value)>::template apply< 1015 typename find_base_dimensions<T>::type, 1016 Dimensions 1017 >::type type; 1018 }; 1019 }; 1020 1021 template<> 1022 struct calculate_base_unit_exponents_impl<false> { 1023 template<class T, class Dimensions> 1024 struct apply { 1025 // find the units that correspond to each base dimension 1026 typedef normalize_units<T> base_solutions; 1027 // pad the dimension with zeroes so it can just be a 1028 // list of numbers, making the multiplication easy 1029 // e.g. if the arguments are list<pound, foot> and 1030 // list<mass,time^-2> then this step will 1031 // yield list<0,1,-2> 1032 typedef typename expand_dimensions<(base_solutions::dimensions::size::value)>::template apply< 1033 typename base_solutions::dimensions, 1034 Dimensions 1035 >::type dimensions; 1036 // take the unit corresponding to each base unit 1037 // multiply each of its exponents by the exponent 1038 // of the base_dimension in the result and sum. 1039 typedef typename multiply_add_units<dimensions::size::value>::template apply< 1040 dimensions, 1041 typename base_solutions::type 1042 >::type units; 1043 // Now, verify that the dummy units really 1044 // cancel out and remove them. 1045 typedef typename strip_zeroes_impl<base_solutions::extra>::template apply<units>::type type; 1046 }; 1047 }; 1048 1049 template<class T, class Dimensions> 1050 struct calculate_base_unit_exponents { 1051 typedef typename calculate_base_unit_exponents_impl<is_simple_system<T>::value>::template apply<T, Dimensions>::type type; 1052 }; 1053 1054 } // namespace detail 1055 1056 } // namespace units 1057 1058 } // namespace boost 1059 1060 #endif 1061