ceph/src/boost/libs/units/include/boost/units/detail/linear_algebra.hpp

   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     static const long size = sizeof(Set::lookup((wrap<T>*)0));
 602     static const 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     static const 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         static const 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     static const 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