units/detail/linear_algebra.hpp

0001 // Boost.Units - A C++ library for zero-overhead dimensional analysis and
0002 // unit/quantity manipulation and conversion
0003 //
0004 // Copyright (C) 2003-2008 Matthias Christian Schabel
0005 // Copyright (C) 2008 Steven Watanabe
0006 //
0007 // Distributed under the Boost Software License, Version 1.0. (See
0008 // accompanying file LICENSE_1_0.txt or copy at
0009 // http://www.boost.org/LICENSE_1_0.txt)
0010
0011 #ifndef BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP
0012 #define BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP
0013
0014 #include <boost/units/static_rational.hpp>
0015 #include <boost/mpl/next.hpp>
0016 #include <boost/mpl/arithmetic.hpp>
0017 #include <boost/mpl/and.hpp>
0018 #include <boost/mpl/assert.hpp>
0019
0020 #include <boost/units/dim.hpp>
0021 #include <boost/units/dimensionless_type.hpp>
0022 #include <boost/units/static_rational.hpp>
0023 #include <boost/units/detail/dimension_list.hpp>
0024 #include <boost/units/detail/sort.hpp>
0025
0026 namespace boost {
0027
0028 namespace units {
0029
0030 namespace detail {
0031
0032 // typedef list<rational> equation;
0033
0034 template<int N>
0035 struct eliminate_from_pair_of_equations_impl;
0036
0037 template<class E1, class E2>
0038 struct eliminate_from_pair_of_equations;
0039
0040 template<int N>
0041 struct elimination_impl;
0042
0043 template<bool is_zero, bool element_is_last>
0044 struct elimination_skip_leading_zeros_impl;
0045
0046 template<class Equation, class Vars>
0047 struct substitute;
0048
0049 template<int N>
0050 struct substitute_impl;
0051
0052 template<bool is_end>
0053 struct solve_impl;
0054
0055 template<class T>
0056 struct solve;
0057
0058 template<int N>
0059 struct check_extra_equations_impl;
0060
0061 template<int N>
0062 struct normalize_units_impl;
0063
0064 struct inconsistent {};
0065
0066 // generally useful utilies.
0067
0068 template<int N>
0069 struct divide_equation {
0070     template<class Begin, class Divisor>
0071     struct apply {
0072         typedef list<typename mpl::divides<typename Begin::item, Divisor>::type, typename divide_equation<N - 1>::template apply<typename Begin::next, Divisor>::type> type;
0073     };
0074 };
0075
0076 template<>
0077 struct divide_equation<0> {
0078     template<class Begin, class Divisor>
0079     struct apply {
0080         typedef dimensionless_type type;
0081     };
0082 };
0083
0084 // eliminate_from_pair_of_equations takes a pair of
0085 // equations and eliminates the first variable.
0086 //
0087 // equation eliminate_from_pair_of_equations(equation l1, equation l2) {
0088 //     rational x1 = l1.front();
0089 //     rational x2 = l2.front();
0090 //     return(transform(pop_front(l1), pop_front(l2), _1 * x2 - _2 * x1));
0091 // }
0092
0093 template<int N>
0094 struct eliminate_from_pair_of_equations_impl {
0095     template<class Begin1, class Begin2, class X1, class X2>
0096     struct apply {
0097         typedef list<
0098             typename mpl::minus<
0099                 typename mpl::times<typename Begin1::item, X2>::type,
0100                 typename mpl::times<typename Begin2::item, X1>::type
0101             >::type,
0102             typename eliminate_from_pair_of_equations_impl<N - 1>::template apply<
0103                 typename Begin1::next,
0104                 typename Begin2::next,
0105                 X1,
0106                 X2
0107             >::type
0108         > type;
0109     };
0110 };
0111
0112 template<>
0113 struct eliminate_from_pair_of_equations_impl<0> {
0114     template<class Begin1, class Begin2, class X1, class X2>
0115     struct apply {
0116         typedef dimensionless_type type;
0117     };
0118 };
0119
0120 template<class E1, class E2>
0121 struct eliminate_from_pair_of_equations {
0122     typedef E1 begin1;
0123     typedef E2 begin2;
0124     typedef typename eliminate_from_pair_of_equations_impl<(E1::size::value - 1)>::template apply<
0125         typename begin1::next,
0126         typename begin2::next,
0127         typename begin1::item,
0128         typename begin2::item
0129     >::type type;
0130 };
0131
0132
0133
0134 // Stage 1.  Determine which dimensions should
0135 // have dummy base units.  For this purpose
0136 // row reduce the matrix.
0137
0138 template<int N>
0139 struct make_zero_vector {
0140     typedef list<static_rational<0>, typename make_zero_vector<N - 1>::type> type;
0141 };
0142 template<>
0143 struct make_zero_vector<0> {
0144     typedef dimensionless_type type;
0145 };
0146
0147 template<int Column, int TotalColumns>
0148 struct create_row_of_identity {
0149     typedef list<static_rational<0>, typename create_row_of_identity<Column - 1, TotalColumns - 1>::type> type;
0150 };
0151 template<int TotalColumns>
0152 struct create_row_of_identity<0, TotalColumns> {
0153     typedef list<static_rational<1>, typename make_zero_vector<TotalColumns - 1>::type> type;
0154 };
0155 template<int Column>
0156 struct create_row_of_identity<Column, 0> {
0157     // error
0158 };
0159
0160 template<int RemainingRows>
0161 struct determine_extra_equations_impl;
0162
0163 template<bool first_is_zero, bool is_last>
0164 struct determine_extra_equations_skip_zeros_impl;
0165
0166 // not the last row and not zero.
0167 template<>
0168 struct determine_extra_equations_skip_zeros_impl<false, false> {
0169     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
0170     struct apply {
0171         // remove the equation being eliminated against from the set of equations.
0172         typedef typename determine_extra_equations_impl<RemainingRows - 1>::template apply<typename RowsBegin::next, typename RowsBegin::item>::type next_equations;
0173         // since this column was present, strip it out.
0174         typedef Result type;
0175     };
0176 };
0177
0178 // the last row but not zero.
0179 template<>
0180 struct determine_extra_equations_skip_zeros_impl<false, true> {
0181     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
0182     struct apply {
0183         // remove this equation.
0184         typedef dimensionless_type next_equations;
0185         // since this column was present, strip it out.
0186         typedef Result type;
0187     };
0188 };
0189
0190
0191 // the first columns is zero but it is not the last column.
0192 // continue with the same loop.
0193 template<>
0194 struct determine_extra_equations_skip_zeros_impl<true, false> {
0195     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
0196     struct apply {
0197         typedef typename RowsBegin::next::item next_row;
0198         typedef typename determine_extra_equations_skip_zeros_impl<
0199             next_row::item::Numerator == 0,
0200             RemainingRows == 2  // the next one will be the last.
0201         >::template apply<
0202             typename RowsBegin::next,
0203             RemainingRows - 1,
0204             CurrentColumn,
0205             TotalColumns,
0206             Result
0207         > next;
0208         typedef list<typename RowsBegin::item::next, typename next::next_equations> next_equations;
0209         typedef typename next::type type;
0210     };
0211 };
0212
0213 // all the elements in this column are zero.
0214 template<>
0215 struct determine_extra_equations_skip_zeros_impl<true, true> {
0216     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
0217     struct apply {
0218         typedef list<typename RowsBegin::item::next, dimensionless_type> next_equations;
0219         typedef list<typename create_row_of_identity<CurrentColumn, TotalColumns>::type, Result> type;
0220     };
0221 };
0222
0223 template<int RemainingRows>
0224 struct determine_extra_equations_impl {
0225     template<class RowsBegin, class EliminateAgainst>
0226     struct apply {
0227         typedef list<
0228             typename eliminate_from_pair_of_equations<typename RowsBegin::item, EliminateAgainst>::type,
0229             typename determine_extra_equations_impl<RemainingRows-1>::template apply<typename RowsBegin::next, EliminateAgainst>::type
0230         > type;
0231     };
0232 };
0233
0234 template<>
0235 struct determine_extra_equations_impl<0> {
0236     template<class RowsBegin, class EliminateAgainst>
0237     struct apply {
0238         typedef dimensionless_type type;
0239     };
0240 };
0241
0242 template<int RemainingColumns, bool is_done>
0243 struct determine_extra_equations {
0244     template<class RowsBegin, int TotalColumns, class Result>
0245     struct apply {
0246         typedef typename RowsBegin::item top_row;
0247         typedef typename determine_extra_equations_skip_zeros_impl<
0248             top_row::item::Numerator == 0,
0249             RowsBegin::size::value == 1
0250         >::template apply<
0251             RowsBegin,
0252             RowsBegin::size::value,
0253             TotalColumns - RemainingColumns,
0254             TotalColumns,
0255             Result
0256         > column_info;
0257         typedef typename determine_extra_equations<
0258             RemainingColumns - 1,
0259             column_info::next_equations::size::value == 0
0260         >::template apply<
0261             typename column_info::next_equations,
0262             TotalColumns,
0263             typename column_info::type
0264         >::type type;
0265     };
0266 };
0267
0268 template<int RemainingColumns>
0269 struct determine_extra_equations<RemainingColumns, true> {
0270     template<class RowsBegin, int TotalColumns, class Result>
0271     struct apply {
0272         typedef typename determine_extra_equations<RemainingColumns - 1, true>::template apply<
0273             RowsBegin,
0274             TotalColumns,
0275             list<typename create_row_of_identity<TotalColumns - RemainingColumns, TotalColumns>::type, Result>
0276         >::type type;
0277     };
0278 };
0279
0280 template<>
0281 struct determine_extra_equations<0, true> {
0282     template<class RowsBegin, int TotalColumns, class Result>
0283     struct apply {
0284         typedef Result type;
0285     };
0286 };
0287
0288 // Stage 2
0289 // invert the matrix using Gauss-Jordan elimination
0290
0291
0292 template<bool is_zero, bool is_last>
0293 struct invert_strip_leading_zeroes;
0294
0295 template<int N>
0296 struct invert_handle_after_pivot_row;
0297
0298 // When processing column N, none of the first N rows
0299 // can be the pivot column.
0300 template<int N>
0301 struct invert_handle_inital_rows {
0302     template<class RowsBegin, class IdentityBegin>
0303     struct apply {
0304         typedef typename invert_handle_inital_rows<N - 1>::template apply<
0305             typename RowsBegin::next,
0306             typename IdentityBegin::next
0307         > next;
0308         typedef typename RowsBegin::item current_row;
0309         typedef typename IdentityBegin::item current_identity_row;
0310         typedef typename next::pivot_row pivot_row;
0311         typedef typename next::identity_pivot_row identity_pivot_row;
0312         typedef list<
0313             typename eliminate_from_pair_of_equations_impl<(current_row::size::value) - 1>::template apply<
0314                 typename current_row::next,
0315                 pivot_row,
0316                 typename current_row::item,
0317                 static_rational<1>
0318             >::type,
0319             typename next::new_matrix
0320         > new_matrix;
0321         typedef list<
0322             typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply<
0323                 current_identity_row,
0324                 identity_pivot_row,
0325                 typename current_row::item,
0326                 static_rational<1>
0327             >::type,
0328             typename next::identity_result
0329         > identity_result;
0330     };
0331 };
0332
0333 // This handles the switch to searching for a pivot column.
0334 // The pivot row will be propagated up in the typedefs
0335 // pivot_row and identity_pivot_row.  It is inserted here.
0336 template<>
0337 struct invert_handle_inital_rows<0> {
0338     template<class RowsBegin, class IdentityBegin>
0339     struct apply {
0340         typedef typename RowsBegin::item current_row;
0341         typedef typename invert_strip_leading_zeroes<
0342             (current_row::item::Numerator == 0),
0343             (RowsBegin::size::value == 1)
0344         >::template apply<
0345             RowsBegin,
0346             IdentityBegin
0347         > next;
0348         // results
0349         typedef list<typename next::pivot_row, typename next::new_matrix> new_matrix;
0350         typedef list<typename next::identity_pivot_row, typename next::identity_result> identity_result;
0351         typedef typename next::pivot_row pivot_row;
0352         typedef typename next::identity_pivot_row identity_pivot_row;
0353     };
0354 };
0355
0356 // The first internal element which is not zero.
0357 template<>
0358 struct invert_strip_leading_zeroes<false, false> {
0359     template<class RowsBegin, class IdentityBegin>
0360     struct apply {
0361         typedef typename RowsBegin::item current_row;
0362         typedef typename current_row::item current_value;
0363         typedef typename divide_equation<(current_row::size::value - 1)>::template apply<typename current_row::next, current_value>::type new_equation;
0364         typedef typename divide_equation<(IdentityBegin::item::size::value)>::template apply<typename IdentityBegin::item, current_value>::type transformed_identity_equation;
0365         typedef typename invert_handle_after_pivot_row<(RowsBegin::size::value - 1)>::template apply<
0366             typename RowsBegin::next,
0367             typename IdentityBegin::next,
0368             new_equation,
0369             transformed_identity_equation
0370         > next;
0371
0372         // results
0373         // Note that we don't add the pivot row to the
0374         // results here, because it needs to propagated up
0375         // to the diagonal.
0376         typedef typename next::new_matrix new_matrix;
0377         typedef typename next::identity_result identity_result;
0378         typedef new_equation pivot_row;
0379         typedef transformed_identity_equation identity_pivot_row;
0380     };
0381 };
0382
0383 // The one and only non-zero element--at the end
0384 template<>
0385 struct invert_strip_leading_zeroes<false, true> {
0386     template<class RowsBegin, class IdentityBegin>
0387     struct apply {
0388         typedef typename RowsBegin::item current_row;
0389         typedef typename current_row::item current_value;
0390         typedef typename divide_equation<(current_row::size::value - 1)>::template apply<typename current_row::next, current_value>::type new_equation;
0391         typedef typename divide_equation<(IdentityBegin::item::size::value)>::template apply<typename IdentityBegin::item, current_value>::type transformed_identity_equation;
0392
0393         // results
0394         // Note that we don't add the pivot row to the
0395         // results here, because it needs to propagated up
0396         // to the diagonal.
0397         typedef dimensionless_type identity_result;
0398         typedef dimensionless_type new_matrix;
0399         typedef new_equation pivot_row;
0400         typedef transformed_identity_equation identity_pivot_row;
0401     };
0402 };
0403
0404 // One of the initial zeroes
0405 template<>
0406 struct invert_strip_leading_zeroes<true, false> {
0407     template<class RowsBegin, class IdentityBegin>
0408     struct apply {
0409         typedef typename RowsBegin::item current_row;
0410         typedef typename RowsBegin::next::item next_row;
0411         typedef typename invert_strip_leading_zeroes<
0412             next_row::item::Numerator == 0,
0413             RowsBegin::size::value == 2
0414         >::template apply<
0415             typename RowsBegin::next,
0416             typename IdentityBegin::next
0417         > next;
0418         typedef typename IdentityBegin::item current_identity_row;
0419         // these are propagated up.
0420         typedef typename next::pivot_row pivot_row;
0421         typedef typename next::identity_pivot_row identity_pivot_row;
0422         typedef list<
0423             typename eliminate_from_pair_of_equations_impl<(current_row::size::value - 1)>::template apply<
0424                 typename current_row::next,
0425                 pivot_row,
0426                 typename current_row::item,
0427                 static_rational<1>
0428             >::type,
0429             typename next::new_matrix
0430         > new_matrix;
0431         typedef list<
0432             typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply<
0433                 current_identity_row,
0434                 identity_pivot_row,
0435                 typename current_row::item,
0436                 static_rational<1>
0437             >::type,
0438             typename next::identity_result
0439         > identity_result;
0440     };
0441 };
0442
0443 // the last element, and is zero.
0444 // Should never happen.
0445 template<>
0446 struct invert_strip_leading_zeroes<true, true> {
0447 };
0448
0449 template<int N>
0450 struct invert_handle_after_pivot_row {
0451     template<class RowsBegin, class IdentityBegin, class MatrixPivot, class IdentityPivot>
0452     struct apply {
0453         typedef typename invert_handle_after_pivot_row<N - 1>::template apply<
0454             typename RowsBegin::next,
0455             typename IdentityBegin::next,
0456             MatrixPivot,
0457             IdentityPivot
0458         > next;
0459         typedef typename RowsBegin::item current_row;
0460         typedef typename IdentityBegin::item current_identity_row;
0461         typedef MatrixPivot pivot_row;
0462         typedef IdentityPivot identity_pivot_row;
0463
0464         // results
0465         typedef list<
0466             typename eliminate_from_pair_of_equations_impl<(current_row::size::value - 1)>::template apply<
0467                 typename current_row::next,
0468                 pivot_row,
0469                 typename current_row::item,
0470                 static_rational<1>
0471             >::type,
0472             typename next::new_matrix
0473         > new_matrix;
0474         typedef list<
0475             typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply<
0476                 current_identity_row,
0477                 identity_pivot_row,
0478                 typename current_row::item,
0479                 static_rational<1>
0480             >::type,
0481             typename next::identity_result
0482         > identity_result;
0483     };
0484 };
0485
0486 template<>
0487 struct invert_handle_after_pivot_row<0> {
0488     template<class RowsBegin, class IdentityBegin, class MatrixPivot, class IdentityPivot>
0489     struct apply {
0490         typedef dimensionless_type new_matrix;
0491         typedef dimensionless_type identity_result;
0492     };
0493 };
0494
0495 template<int N>
0496 struct invert_impl {
0497     template<class RowsBegin, class IdentityBegin>
0498     struct apply {
0499         typedef typename invert_handle_inital_rows<RowsBegin::size::value - N>::template apply<RowsBegin, IdentityBegin> process_column;
0500         typedef typename invert_impl<N - 1>::template apply<
0501             typename process_column::new_matrix,
0502             typename process_column::identity_result
0503         >::type type;
0504     };
0505 };
0506
0507 template<>
0508 struct invert_impl<0> {
0509     template<class RowsBegin, class IdentityBegin>
0510     struct apply {
0511         typedef IdentityBegin type;
0512     };
0513 };
0514
0515 template<int N>
0516 struct make_identity {
0517     template<int Size>
0518     struct apply {
0519         typedef list<typename create_row_of_identity<Size - N, Size>::type, typename make_identity<N - 1>::template apply<Size>::type> type;
0520     };
0521 };
0522
0523 template<>
0524 struct make_identity<0> {
0525     template<int Size>
0526     struct apply {
0527         typedef dimensionless_type type;
0528     };
0529 };
0530
0531 template<class Matrix>
0532 struct make_square_and_invert {
0533     typedef typename Matrix::item top_row;
0534     typedef typename determine_extra_equations<(top_row::size::value), false>::template apply<
0535         Matrix,                 // RowsBegin
0536         top_row::size::value,   // TotalColumns
0537         Matrix                  // Result
0538     >::type invertible;
0539     typedef typename invert_impl<invertible::size::value>::template apply<
0540         invertible,
0541         typename make_identity<invertible::size::value>::template apply<invertible::size::value>::type
0542     >::type type;
0543 };
0544
0545
0546 // find_base_dimensions takes a list of
0547 // base_units and returns a sorted list
0548 // of all the base_dimensions they use.
0549 //
0550 // list<base_dimension> find_base_dimensions(list<base_unit> l) {
0551 //     set<base_dimension> dimensions;
0552 //     for_each(base_unit unit : l) {
0553 //         for_each(dim d : unit.dimension_type) {
0554 //             dimensions = insert(dimensions, d.tag_type);
0555 //         }
0556 //     }
0557 //     return(sort(dimensions, _1 > _2, front_inserter(list<base_dimension>())));
0558 // }
0559
0560 typedef char set_no;
0561 struct set_yes { set_no dummy[2]; };
0562
0563 template<class T>
0564 struct wrap {};
0565
0566 struct set_end {
0567     static set_no lookup(...);
0568     typedef mpl::long_<0> size;
0569 };
0570
0571 template<class T, class Next>
0572 struct set : Next {
0573     using Next::lookup;
0574     static set_yes lookup(wrap<T>*);
0575     typedef T item;
0576     typedef Next next;
0577     typedef typename mpl::next<typename Next::size>::type size;
0578 };
0579
0580 template<bool has_key>
0581 struct set_insert;
0582
0583 template<>
0584 struct set_insert<true> {
0585     template<class Set, class T>
0586     struct apply {
0587         typedef Set type;
0588     };
0589 };
0590
0591 template<>
0592 struct set_insert<false> {
0593     template<class Set, class T>
0594     struct apply {
0595         typedef set<T, Set> type;
0596     };
0597 };
0598
0599 template<class Set, class T>
0600 struct has_key {
0601     BOOST_STATIC_CONSTEXPR long size = sizeof(Set::lookup((wrap<T>*)0));
0602     BOOST_STATIC_CONSTEXPR bool value = (size == sizeof(set_yes));
0603 };
0604
0605 template<int N>
0606 struct find_base_dimensions_impl_impl {
0607     template<class Begin, class S>
0608     struct apply {
0609         typedef typename find_base_dimensions_impl_impl<N-1>::template apply<
0610             typename Begin::next,
0611             S
0612         >::type next;
0613
0614         typedef typename set_insert<
0615             (has_key<next, typename Begin::item::tag_type>::value)
0616         >::template apply<
0617             next,
0618             typename Begin::item::tag_type
0619         >::type type;
0620     };
0621 };
0622
0623 template<>
0624 struct find_base_dimensions_impl_impl<0> {
0625     template<class Begin, class S>
0626     struct apply {
0627         typedef S type;
0628     };
0629 };
0630
0631 template<int N>
0632 struct find_base_dimensions_impl {
0633     template<class Begin>
0634     struct apply {
0635         typedef typename find_base_dimensions_impl_impl<(Begin::item::dimension_type::size::value)>::template apply<
0636             typename Begin::item::dimension_type,
0637             typename find_base_dimensions_impl<N-1>::template apply<typename Begin::next>::type
0638         >::type type;
0639     };
0640 };
0641
0642 template<>
0643 struct find_base_dimensions_impl<0> {
0644     template<class Begin>
0645     struct apply {
0646         typedef set_end type;
0647     };
0648 };
0649
0650 template<class T>
0651 struct find_base_dimensions {
0652     typedef typename insertion_sort<
0653         typename find_base_dimensions_impl<
0654             (T::size::value)
0655         >::template apply<T>::type
0656     >::type type;
0657 };
0658
0659 // calculate_base_dimension_coefficients finds
0660 // the coefficients corresponding to the first
0661 // base_dimension in each of the dimension_lists.
0662 // It returns two values.  The first result
0663 // is a list of the coefficients.  The second
0664 // is a list with all the incremented iterators.
0665 // When we encounter a base_dimension that is
0666 // missing from a dimension_list, we do not
0667 // increment the iterator and we set the
0668 // coefficient to zero.
0669
0670 template<bool has_dimension>
0671 struct calculate_base_dimension_coefficients_func;
0672
0673 template<>
0674 struct calculate_base_dimension_coefficients_func<true> {
0675     template<class T>
0676     struct apply {
0677         typedef typename T::item::value_type type;
0678         typedef typename T::next next;
0679     };
0680 };
0681
0682 template<>
0683 struct calculate_base_dimension_coefficients_func<false> {
0684     template<class T>
0685     struct apply {
0686         typedef static_rational<0> type;
0687         typedef T next;
0688     };
0689 };
0690
0691 // begins_with_dimension returns true iff its first
0692 // parameter is a valid iterator which yields its
0693 // second parameter when dereferenced.
0694
0695 template<class Iterator>
0696 struct begins_with_dimension {
0697     template<class Dim>
0698     struct apply :
0699         boost::is_same<
0700             Dim,
0701             typename Iterator::item::tag_type
0702         > {};
0703 };
0704
0705 template<>
0706 struct begins_with_dimension<dimensionless_type> {
0707     template<class Dim>
0708     struct apply : mpl::false_ {};
0709 };
0710
0711 template<int N>
0712 struct calculate_base_dimension_coefficients_impl {
0713     template<class BaseUnitDimensions,class Dim,class T>
0714     struct apply {
0715         typedef typename calculate_base_dimension_coefficients_func<
0716             begins_with_dimension<typename BaseUnitDimensions::item>::template apply<
0717                 Dim
0718             >::value
0719         >::template apply<
0720             typename BaseUnitDimensions::item
0721         > result;
0722         typedef typename calculate_base_dimension_coefficients_impl<N-1>::template apply<
0723             typename BaseUnitDimensions::next,
0724             Dim,
0725             list<typename result::type, T>
0726         > next_;
0727         typedef typename next_::type type;
0728         typedef list<typename result::next, typename next_::next> next;
0729     };
0730 };
0731
0732 template<>
0733 struct calculate_base_dimension_coefficients_impl<0> {
0734     template<class Begin, class BaseUnitDimensions, class T>
0735     struct apply {
0736         typedef T type;
0737         typedef dimensionless_type next;
0738     };
0739 };
0740
0741 // add_zeroes pushs N zeroes onto the
0742 // front of a list.
0743 //
0744 // list<rational> add_zeroes(list<rational> l, int N) {
0745 //     if(N == 0) {
0746 //         return(l);
0747 //     } else {
0748 //         return(push_front(add_zeroes(l, N-1), 0));
0749 //     }
0750 // }
0751
0752 template<int N>
0753 struct add_zeroes_impl {
0754     // If you get an error here and your base units are
0755     // in fact linearly independent, please report it.
0756     BOOST_MPL_ASSERT_MSG((N > 0), base_units_are_probably_not_linearly_independent, (void));
0757     template<class T>
0758     struct apply {
0759         typedef list<
0760             static_rational<0>,
0761             typename add_zeroes_impl<N-1>::template apply<T>::type
0762         > type;
0763     };
0764 };
0765
0766 template<>
0767 struct add_zeroes_impl<0> {
0768     template<class T>
0769     struct apply {
0770         typedef T type;
0771     };
0772 };
0773
0774 // expand_dimensions finds the exponents of
0775 // a set of dimensions in a dimension_list.
0776 // the second parameter is assumed to be
0777 // a superset of the base_dimensions of
0778 // the first parameter.
0779 //
0780 // list<rational> expand_dimensions(dimension_list, list<base_dimension>);
0781
0782 template<int N>
0783 struct expand_dimensions {
0784     template<class Begin, class DimensionIterator>
0785     struct apply {
0786         typedef typename calculate_base_dimension_coefficients_func<
0787             begins_with_dimension<DimensionIterator>::template apply<typename Begin::item>::value
0788         >::template apply<DimensionIterator> result;
0789         typedef list<
0790             typename result::type,
0791             typename expand_dimensions<N-1>::template apply<typename Begin::next, typename result::next>::type
0792         > type;
0793     };
0794 };
0795
0796 template<>
0797 struct expand_dimensions<0> {
0798     template<class Begin, class DimensionIterator>
0799     struct apply {
0800         typedef dimensionless_type type;
0801     };
0802 };
0803
0804 template<int N>
0805 struct create_unit_matrix {
0806     template<class Begin, class Dimensions>
0807     struct apply {
0808         typedef typename create_unit_matrix<N - 1>::template apply<typename Begin::next, Dimensions>::type next;
0809         typedef list<typename expand_dimensions<Dimensions::size::value>::template apply<Dimensions, typename Begin::item::dimension_type>::type, next> type;
0810     };
0811 };
0812
0813 template<>
0814 struct create_unit_matrix<0> {
0815     template<class Begin, class Dimensions>
0816     struct apply {
0817         typedef dimensionless_type type;
0818     };
0819 };
0820
0821 template<class T>
0822 struct normalize_units {
0823     typedef typename find_base_dimensions<T>::type dimensions;
0824     typedef typename create_unit_matrix<(T::size::value)>::template apply<
0825         T,
0826         dimensions
0827     >::type matrix;
0828     typedef typename make_square_and_invert<matrix>::type type;
0829     BOOST_STATIC_CONSTEXPR long extra = (type::size::value) - (T::size::value);
0830 };
0831
0832 // multiply_add_units computes M x V
0833 // where M is a matrix and V is a horizontal
0834 // vector
0835 //
0836 // list<rational> multiply_add_units(list<list<rational> >, list<rational>);
0837
0838 template<int N>
0839 struct multiply_add_units_impl {
0840     template<class Begin1, class Begin2 ,class X>
0841     struct apply {
0842         typedef list<
0843             typename mpl::plus<
0844                 typename mpl::times<
0845                     typename Begin2::item,
0846                     X
0847                 >::type,
0848                 typename Begin1::item
0849             >::type,
0850             typename multiply_add_units_impl<N-1>::template apply<
0851                 typename Begin1::next,
0852                 typename Begin2::next,
0853                 X
0854             >::type
0855         > type;
0856     };
0857 };
0858
0859 template<>
0860 struct multiply_add_units_impl<0> {
0861     template<class Begin1, class Begin2 ,class X>
0862     struct apply {
0863         typedef dimensionless_type type;
0864     };
0865 };
0866
0867 template<int N>
0868 struct multiply_add_units {
0869     template<class Begin1, class Begin2>
0870     struct apply {
0871         typedef typename multiply_add_units_impl<
0872             (Begin2::item::size::value)
0873         >::template apply<
0874             typename multiply_add_units<N-1>::template apply<
0875                 typename Begin1::next,
0876                 typename Begin2::next
0877             >::type,
0878             typename Begin2::item,
0879             typename Begin1::item
0880         >::type type;
0881     };
0882 };
0883
0884 template<>
0885 struct multiply_add_units<1> {
0886     template<class Begin1, class Begin2>
0887     struct apply {
0888         typedef typename add_zeroes_impl<
0889             (Begin2::item::size::value)
0890         >::template apply<dimensionless_type>::type type1;
0891         typedef typename multiply_add_units_impl<
0892             (Begin2::item::size::value)
0893         >::template apply<
0894             type1,
0895             typename Begin2::item,
0896             typename Begin1::item
0897         >::type type;
0898     };
0899 };
0900
0901
0902 // strip_zeroes erases the first N elements of a list if
0903 // they are all zero, otherwise returns inconsistent
0904 //
0905 // list strip_zeroes(list l, int N) {
0906 //     if(N == 0) {
0907 //         return(l);
0908 //     } else if(l.front == 0) {
0909 //         return(strip_zeroes(pop_front(l), N-1));
0910 //     } else {
0911 //         return(inconsistent);
0912 //     }
0913 // }
0914
0915 template<int N>
0916 struct strip_zeroes_impl;
0917
0918 template<class T>
0919 struct strip_zeroes_func {
0920     template<class L, int N>
0921     struct apply {
0922         typedef inconsistent type;
0923     };
0924 };
0925
0926 template<>
0927 struct strip_zeroes_func<static_rational<0> > {
0928     template<class L, int N>
0929     struct apply {
0930         typedef typename strip_zeroes_impl<N-1>::template apply<typename L::next>::type type;
0931     };
0932 };
0933
0934 template<int N>
0935 struct strip_zeroes_impl {
0936     template<class T>
0937     struct apply {
0938         typedef typename strip_zeroes_func<typename T::item>::template apply<T, N>::type type;
0939     };
0940 };
0941
0942 template<>
0943 struct strip_zeroes_impl<0> {
0944     template<class T>
0945     struct apply {
0946         typedef T type;
0947     };
0948 };
0949
0950 // Given a list of base_units, computes the
0951 // exponents of each base unit for a given
0952 // dimension.
0953 //
0954 // list<rational> calculate_base_unit_exponents(list<base_unit> units, dimension_list dimensions);
0955
0956 template<class T>
0957 struct is_base_dimension_unit {
0958     typedef mpl::false_ type;
0959     typedef void base_dimension_type;
0960 };
0961 template<class T>
0962 struct is_base_dimension_unit<list<dim<T, static_rational<1> >, dimensionless_type> > {
0963     typedef mpl::true_ type;
0964     typedef T base_dimension_type;
0965 };
0966
0967 template<int N>
0968 struct is_simple_system_impl {
0969     template<class Begin, class Prev>
0970     struct apply {
0971         typedef is_base_dimension_unit<typename Begin::item::dimension_type> test;
0972         typedef mpl::and_<
0973             typename test::type,
0974             mpl::less<Prev, typename test::base_dimension_type>,
0975             typename is_simple_system_impl<N-1>::template apply<
0976                 typename Begin::next,
0977                 typename test::base_dimension_type
0978             >
0979         > type;
0980         BOOST_STATIC_CONSTEXPR bool value = (type::value);
0981     };
0982 };
0983
0984 template<>
0985 struct is_simple_system_impl<0> {
0986     template<class Begin, class Prev>
0987     struct apply : mpl::true_ {
0988     };
0989 };
0990
0991 template<class T>
0992 struct is_simple_system {
0993     typedef T Begin;
0994     typedef is_base_dimension_unit<typename Begin::item::dimension_type> test;
0995     typedef typename mpl::and_<
0996         typename test::type,
0997         typename is_simple_system_impl<
0998             T::size::value - 1
0999         >::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