Chris@16: // Boost.Units - A C++ library for zero-overhead dimensional analysis and 
Chris@16: // unit/quantity manipulation and conversion
Chris@16: //
Chris@16: // Copyright (C) 2003-2008 Matthias Christian Schabel
Chris@16: // Copyright (C) 2008 Steven Watanabe
Chris@16: //
Chris@16: // Distributed under the Boost Software License, Version 1.0. (See
Chris@16: // accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP
Chris@16: #define BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP
Chris@16: 
Chris@16: #include <boost/units/static_rational.hpp>
Chris@16: #include <boost/mpl/next.hpp>
Chris@16: #include <boost/mpl/arithmetic.hpp>
Chris@16: #include <boost/mpl/and.hpp>
Chris@16: #include <boost/mpl/assert.hpp>
Chris@16: 
Chris@16: #include <boost/units/dim.hpp>
Chris@16: #include <boost/units/dimensionless_type.hpp>
Chris@16: #include <boost/units/static_rational.hpp>
Chris@16: #include <boost/units/detail/dimension_list.hpp>
Chris@16: #include <boost/units/detail/sort.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: 
Chris@16: namespace units {
Chris@16: 
Chris@16: namespace detail {
Chris@16: 
Chris@16: // typedef list<rational> equation;
Chris@16: 
Chris@16: template<int N>
Chris@16: struct eliminate_from_pair_of_equations_impl;
Chris@16: 
Chris@16: template<class E1, class E2>
Chris@16: struct eliminate_from_pair_of_equations;
Chris@16: 
Chris@16: template<int N>
Chris@16: struct elimination_impl;
Chris@16: 
Chris@16: template<bool is_zero, bool element_is_last>
Chris@16: struct elimination_skip_leading_zeros_impl;
Chris@16: 
Chris@16: template<class Equation, class Vars>
Chris@16: struct substitute;
Chris@16: 
Chris@16: template<int N>
Chris@16: struct substitute_impl;
Chris@16: 
Chris@16: template<bool is_end>
Chris@16: struct solve_impl;
Chris@16: 
Chris@16: template<class T>
Chris@16: struct solve;
Chris@16: 
Chris@16: template<int N>
Chris@16: struct check_extra_equations_impl;
Chris@16: 
Chris@16: template<int N>
Chris@16: struct normalize_units_impl;
Chris@16: 
Chris@16: struct inconsistent {};
Chris@16: 
Chris@16: // generally useful utilies.
Chris@16: 
Chris@16: template<int N>
Chris@16: struct divide_equation {
Chris@16:     template<class Begin, class Divisor>
Chris@16:     struct apply {
Chris@16:         typedef list<typename mpl::divides<typename Begin::item, Divisor>::type, typename divide_equation<N - 1>::template apply<typename Begin::next, Divisor>::type> type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct divide_equation<0> {
Chris@16:     template<class Begin, class Divisor>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // eliminate_from_pair_of_equations takes a pair of
Chris@16: // equations and eliminates the first variable.
Chris@16: //
Chris@16: // equation eliminate_from_pair_of_equations(equation l1, equation l2) {
Chris@16: //     rational x1 = l1.front();
Chris@16: //     rational x2 = l2.front();
Chris@16: //     return(transform(pop_front(l1), pop_front(l2), _1 * x2 - _2 * x1));
Chris@16: // }
Chris@16: 
Chris@16: template<int N>
Chris@16: struct eliminate_from_pair_of_equations_impl {
Chris@16:     template<class Begin1, class Begin2, class X1, class X2>
Chris@16:     struct apply {
Chris@16:         typedef list<
Chris@16:             typename mpl::minus<
Chris@16:                 typename mpl::times<typename Begin1::item, X2>::type,
Chris@16:                 typename mpl::times<typename Begin2::item, X1>::type
Chris@16:             >::type,
Chris@16:             typename eliminate_from_pair_of_equations_impl<N - 1>::template apply<
Chris@16:                 typename Begin1::next,
Chris@16:                 typename Begin2::next,
Chris@16:                 X1,
Chris@16:                 X2
Chris@16:             >::type
Chris@16:         > type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct eliminate_from_pair_of_equations_impl<0> {
Chris@16:     template<class Begin1, class Begin2, class X1, class X2>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<class E1, class E2>
Chris@16: struct eliminate_from_pair_of_equations {
Chris@16:     typedef E1 begin1;
Chris@16:     typedef E2 begin2;
Chris@16:     typedef typename eliminate_from_pair_of_equations_impl<(E1::size::value - 1)>::template apply<
Chris@16:         typename begin1::next,
Chris@16:         typename begin2::next,
Chris@16:         typename begin1::item,
Chris@16:         typename begin2::item
Chris@16:     >::type type;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: // Stage 1.  Determine which dimensions should
Chris@16: // have dummy base units.  For this purpose
Chris@16: // row reduce the matrix.
Chris@16: 
Chris@16: template<int N>
Chris@16: struct make_zero_vector {
Chris@16:     typedef list<static_rational<0>, typename make_zero_vector<N - 1>::type> type;
Chris@16: };
Chris@16: template<>
Chris@16: struct make_zero_vector<0> {
Chris@16:     typedef dimensionless_type type;
Chris@16: };
Chris@16: 
Chris@16: template<int Column, int TotalColumns>
Chris@16: struct create_row_of_identity {
Chris@16:     typedef list<static_rational<0>, typename create_row_of_identity<Column - 1, TotalColumns - 1>::type> type;
Chris@16: };
Chris@16: template<int TotalColumns>
Chris@16: struct create_row_of_identity<0, TotalColumns> {
Chris@16:     typedef list<static_rational<1>, typename make_zero_vector<TotalColumns - 1>::type> type;
Chris@16: };
Chris@16: template<int Column>
Chris@16: struct create_row_of_identity<Column, 0> {
Chris@16:     // error
Chris@16: };
Chris@16: 
Chris@16: template<int RemainingRows>
Chris@16: struct determine_extra_equations_impl;
Chris@16: 
Chris@16: template<bool first_is_zero, bool is_last>
Chris@16: struct determine_extra_equations_skip_zeros_impl;
Chris@16: 
Chris@16: // not the last row and not zero.
Chris@16: template<>
Chris@16: struct determine_extra_equations_skip_zeros_impl<false, false> {
Chris@16:     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
Chris@16:     struct apply {
Chris@16:         // remove the equation being eliminated against from the set of equations.
Chris@16:         typedef typename determine_extra_equations_impl<RemainingRows - 1>::template apply<typename RowsBegin::next, typename RowsBegin::item>::type next_equations;
Chris@16:         // since this column was present, strip it out.
Chris@16:         typedef Result type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // the last row but not zero.
Chris@16: template<>
Chris@16: struct determine_extra_equations_skip_zeros_impl<false, true> {
Chris@16:     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
Chris@16:     struct apply {
Chris@16:         // remove this equation.
Chris@16:         typedef dimensionless_type next_equations;
Chris@16:         // since this column was present, strip it out.
Chris@16:         typedef Result type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: // the first columns is zero but it is not the last column.
Chris@16: // continue with the same loop.
Chris@16: template<>
Chris@16: struct determine_extra_equations_skip_zeros_impl<true, false> {
Chris@16:     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
Chris@16:     struct apply {
Chris@16:         typedef typename RowsBegin::next::item next_row;
Chris@16:         typedef typename determine_extra_equations_skip_zeros_impl<
Chris@16:             next_row::item::Numerator == 0,
Chris@16:             RemainingRows == 2  // the next one will be the last.
Chris@16:         >::template apply<
Chris@16:             typename RowsBegin::next,
Chris@16:             RemainingRows - 1,
Chris@16:             CurrentColumn,
Chris@16:             TotalColumns,
Chris@16:             Result
Chris@16:         > next;
Chris@16:         typedef list<typename RowsBegin::item::next, typename next::next_equations> next_equations;
Chris@16:         typedef typename next::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // all the elements in this column are zero.
Chris@16: template<>
Chris@16: struct determine_extra_equations_skip_zeros_impl<true, true> {
Chris@16:     template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
Chris@16:     struct apply {
Chris@16:         typedef list<typename RowsBegin::item::next, dimensionless_type> next_equations;
Chris@16:         typedef list<typename create_row_of_identity<CurrentColumn, TotalColumns>::type, Result> type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int RemainingRows>
Chris@16: struct determine_extra_equations_impl {
Chris@16:     template<class RowsBegin, class EliminateAgainst>
Chris@16:     struct apply {
Chris@16:         typedef list<
Chris@16:             typename eliminate_from_pair_of_equations<typename RowsBegin::item, EliminateAgainst>::type,
Chris@16:             typename determine_extra_equations_impl<RemainingRows-1>::template apply<typename RowsBegin::next, EliminateAgainst>::type
Chris@16:         > type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct determine_extra_equations_impl<0> {
Chris@16:     template<class RowsBegin, class EliminateAgainst>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int RemainingColumns, bool is_done>
Chris@16: struct determine_extra_equations {
Chris@16:     template<class RowsBegin, int TotalColumns, class Result>
Chris@16:     struct apply {
Chris@16:         typedef typename RowsBegin::item top_row;
Chris@16:         typedef typename determine_extra_equations_skip_zeros_impl<
Chris@16:             top_row::item::Numerator == 0,
Chris@101:             RowsBegin::size::value == 1
Chris@16:         >::template apply<
Chris@16:             RowsBegin,
Chris@16:             RowsBegin::size::value,
Chris@16:             TotalColumns - RemainingColumns,
Chris@16:             TotalColumns,
Chris@16:             Result
Chris@16:         > column_info;
Chris@16:         typedef typename determine_extra_equations<
Chris@16:             RemainingColumns - 1,
Chris@16:             column_info::next_equations::size::value == 0
Chris@16:         >::template apply<
Chris@16:             typename column_info::next_equations,
Chris@16:             TotalColumns,
Chris@16:             typename column_info::type
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int RemainingColumns>
Chris@16: struct determine_extra_equations<RemainingColumns, true> {
Chris@16:     template<class RowsBegin, int TotalColumns, class Result>
Chris@16:     struct apply {
Chris@16:         typedef typename determine_extra_equations<RemainingColumns - 1, true>::template apply<
Chris@16:             RowsBegin,
Chris@16:             TotalColumns,
Chris@16:             list<typename create_row_of_identity<TotalColumns - RemainingColumns, TotalColumns>::type, Result>
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct determine_extra_equations<0, true> {
Chris@16:     template<class RowsBegin, int TotalColumns, class Result>
Chris@16:     struct apply {
Chris@16:         typedef Result type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // Stage 2
Chris@16: // invert the matrix using Gauss-Jordan elimination
Chris@16: 
Chris@16: 
Chris@16: template<bool is_zero, bool is_last>
Chris@16: struct invert_strip_leading_zeroes;
Chris@16: 
Chris@16: template<int N>
Chris@16: struct invert_handle_after_pivot_row;
Chris@16: 
Chris@16: // When processing column N, none of the first N rows 
Chris@16: // can be the pivot column.
Chris@16: template<int N>
Chris@16: struct invert_handle_inital_rows {
Chris@16:     template<class RowsBegin, class IdentityBegin>
Chris@16:     struct apply {
Chris@16:         typedef typename invert_handle_inital_rows<N - 1>::template apply<
Chris@16:             typename RowsBegin::next,
Chris@16:             typename IdentityBegin::next
Chris@16:         > next;
Chris@16:         typedef typename RowsBegin::item current_row;
Chris@16:         typedef typename IdentityBegin::item current_identity_row;
Chris@16:         typedef typename next::pivot_row pivot_row;
Chris@16:         typedef typename next::identity_pivot_row identity_pivot_row;
Chris@16:         typedef list<
Chris@16:             typename eliminate_from_pair_of_equations_impl<(current_row::size::value) - 1>::template apply<
Chris@16:                 typename current_row::next,
Chris@16:                 pivot_row,
Chris@16:                 typename current_row::item,
Chris@16:                 static_rational<1>
Chris@16:             >::type,
Chris@16:             typename next::new_matrix
Chris@16:         > new_matrix;
Chris@16:         typedef list<
Chris@16:             typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply<
Chris@16:                 current_identity_row,
Chris@16:                 identity_pivot_row,
Chris@16:                 typename current_row::item,
Chris@16:                 static_rational<1>
Chris@16:             >::type,
Chris@16:             typename next::identity_result
Chris@16:         > identity_result;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // This handles the switch to searching for a pivot column.
Chris@16: // The pivot row will be propagated up in the typedefs
Chris@16: // pivot_row and identity_pivot_row.  It is inserted here.
Chris@16: template<>
Chris@16: struct invert_handle_inital_rows<0> {
Chris@16:     template<class RowsBegin, class IdentityBegin>
Chris@16:     struct apply {
Chris@16:         typedef typename RowsBegin::item current_row;
Chris@16:         typedef typename invert_strip_leading_zeroes<
Chris@16:             (current_row::item::Numerator == 0),
Chris@16:             (RowsBegin::size::value == 1)
Chris@16:         >::template apply<
Chris@16:             RowsBegin,
Chris@16:             IdentityBegin
Chris@16:         > next;
Chris@16:         // results
Chris@16:         typedef list<typename next::pivot_row, typename next::new_matrix> new_matrix;
Chris@16:         typedef list<typename next::identity_pivot_row, typename next::identity_result> identity_result;
Chris@16:         typedef typename next::pivot_row pivot_row;
Chris@16:         typedef typename next::identity_pivot_row identity_pivot_row;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // The first internal element which is not zero.
Chris@16: template<>
Chris@16: struct invert_strip_leading_zeroes<false, false> {
Chris@16:     template<class RowsBegin, class IdentityBegin>
Chris@16:     struct apply {
Chris@16:         typedef typename RowsBegin::item current_row;
Chris@16:         typedef typename current_row::item current_value;
Chris@16:         typedef typename divide_equation<(current_row::size::value - 1)>::template apply<typename current_row::next, current_value>::type new_equation;
Chris@16:         typedef typename divide_equation<(IdentityBegin::item::size::value)>::template apply<typename IdentityBegin::item, current_value>::type transformed_identity_equation;
Chris@16:         typedef typename invert_handle_after_pivot_row<(RowsBegin::size::value - 1)>::template apply<
Chris@16:             typename RowsBegin::next,
Chris@16:             typename IdentityBegin::next,
Chris@16:             new_equation,
Chris@16:             transformed_identity_equation
Chris@16:         > next;
Chris@16: 
Chris@16:         // results
Chris@16:         // Note that we don't add the pivot row to the
Chris@16:         // results here, because it needs to propagated up
Chris@16:         // to the diagonal.
Chris@16:         typedef typename next::new_matrix new_matrix;
Chris@16:         typedef typename next::identity_result identity_result;
Chris@16:         typedef new_equation pivot_row;
Chris@16:         typedef transformed_identity_equation identity_pivot_row;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // The one and only non-zero element--at the end
Chris@16: template<>
Chris@16: struct invert_strip_leading_zeroes<false, true> {
Chris@16:     template<class RowsBegin, class IdentityBegin>
Chris@16:     struct apply {
Chris@16:         typedef typename RowsBegin::item current_row;
Chris@16:         typedef typename current_row::item current_value;
Chris@16:         typedef typename divide_equation<(current_row::size::value - 1)>::template apply<typename current_row::next, current_value>::type new_equation;
Chris@16:         typedef typename divide_equation<(IdentityBegin::item::size::value)>::template apply<typename IdentityBegin::item, current_value>::type transformed_identity_equation;
Chris@16: 
Chris@16:         // results
Chris@16:         // Note that we don't add the pivot row to the
Chris@16:         // results here, because it needs to propagated up
Chris@16:         // to the diagonal.
Chris@16:         typedef dimensionless_type identity_result;
Chris@16:         typedef dimensionless_type new_matrix;
Chris@16:         typedef new_equation pivot_row;
Chris@16:         typedef transformed_identity_equation identity_pivot_row;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // One of the initial zeroes
Chris@16: template<>
Chris@16: struct invert_strip_leading_zeroes<true, false> {
Chris@16:     template<class RowsBegin, class IdentityBegin>
Chris@16:     struct apply {
Chris@16:         typedef typename RowsBegin::item current_row;
Chris@16:         typedef typename RowsBegin::next::item next_row;
Chris@16:         typedef typename invert_strip_leading_zeroes<
Chris@16:             next_row::item::Numerator == 0,
Chris@16:             RowsBegin::size::value == 2
Chris@16:         >::template apply<
Chris@16:             typename RowsBegin::next,
Chris@16:             typename IdentityBegin::next
Chris@16:         > next;
Chris@16:         typedef typename IdentityBegin::item current_identity_row;
Chris@16:         // these are propagated up.
Chris@16:         typedef typename next::pivot_row pivot_row;
Chris@16:         typedef typename next::identity_pivot_row identity_pivot_row;
Chris@16:         typedef list<
Chris@16:             typename eliminate_from_pair_of_equations_impl<(current_row::size::value - 1)>::template apply<
Chris@16:                 typename current_row::next,
Chris@16:                 pivot_row,
Chris@16:                 typename current_row::item,
Chris@16:                 static_rational<1>
Chris@16:             >::type,
Chris@16:             typename next::new_matrix
Chris@16:         > new_matrix;
Chris@16:         typedef list<
Chris@16:             typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply<
Chris@16:                 current_identity_row,
Chris@16:                 identity_pivot_row,
Chris@16:                 typename current_row::item,
Chris@16:                 static_rational<1>
Chris@16:             >::type,
Chris@16:             typename next::identity_result
Chris@16:         > identity_result;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // the last element, and is zero.
Chris@16: // Should never happen.
Chris@16: template<>
Chris@16: struct invert_strip_leading_zeroes<true, true> {
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct invert_handle_after_pivot_row {
Chris@16:     template<class RowsBegin, class IdentityBegin, class MatrixPivot, class IdentityPivot>
Chris@16:     struct apply {
Chris@16:         typedef typename invert_handle_after_pivot_row<N - 1>::template apply<
Chris@16:             typename RowsBegin::next,
Chris@16:             typename IdentityBegin::next,
Chris@16:             MatrixPivot,
Chris@16:             IdentityPivot
Chris@16:         > next;
Chris@16:         typedef typename RowsBegin::item current_row;
Chris@16:         typedef typename IdentityBegin::item current_identity_row;
Chris@16:         typedef MatrixPivot pivot_row;
Chris@16:         typedef IdentityPivot identity_pivot_row;
Chris@16: 
Chris@16:         // results
Chris@16:         typedef list<
Chris@16:             typename eliminate_from_pair_of_equations_impl<(current_row::size::value - 1)>::template apply<
Chris@16:                 typename current_row::next,
Chris@16:                 pivot_row,
Chris@16:                 typename current_row::item,
Chris@16:                 static_rational<1>
Chris@16:             >::type,
Chris@16:             typename next::new_matrix
Chris@16:         > new_matrix;
Chris@16:         typedef list<
Chris@16:             typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply<
Chris@16:                 current_identity_row,
Chris@16:                 identity_pivot_row,
Chris@16:                 typename current_row::item,
Chris@16:                 static_rational<1>
Chris@16:             >::type,
Chris@16:             typename next::identity_result
Chris@16:         > identity_result;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct invert_handle_after_pivot_row<0> {
Chris@16:     template<class RowsBegin, class IdentityBegin, class MatrixPivot, class IdentityPivot>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type new_matrix;
Chris@16:         typedef dimensionless_type identity_result;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct invert_impl {
Chris@16:     template<class RowsBegin, class IdentityBegin>
Chris@16:     struct apply {
Chris@16:         typedef typename invert_handle_inital_rows<RowsBegin::size::value - N>::template apply<RowsBegin, IdentityBegin> process_column;
Chris@16:         typedef typename invert_impl<N - 1>::template apply<
Chris@16:             typename process_column::new_matrix,
Chris@16:             typename process_column::identity_result
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct invert_impl<0> {
Chris@16:     template<class RowsBegin, class IdentityBegin>
Chris@16:     struct apply {
Chris@16:         typedef IdentityBegin type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct make_identity {
Chris@16:     template<int Size>
Chris@16:     struct apply {
Chris@16:         typedef list<typename create_row_of_identity<Size - N, Size>::type, typename make_identity<N - 1>::template apply<Size>::type> type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct make_identity<0> {
Chris@16:     template<int Size>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<class Matrix>
Chris@16: struct make_square_and_invert {
Chris@16:     typedef typename Matrix::item top_row;
Chris@16:     typedef typename determine_extra_equations<(top_row::size::value), false>::template apply<
Chris@16:         Matrix,                 // RowsBegin
Chris@16:         top_row::size::value,   // TotalColumns
Chris@16:         Matrix                  // Result
Chris@16:     >::type invertible;
Chris@16:     typedef typename invert_impl<invertible::size::value>::template apply<
Chris@16:         invertible,
Chris@16:         typename make_identity<invertible::size::value>::template apply<invertible::size::value>::type
Chris@16:     >::type type;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: // find_base_dimensions takes a list of
Chris@16: // base_units and returns a sorted list
Chris@16: // of all the base_dimensions they use.
Chris@16: //
Chris@16: // list<base_dimension> find_base_dimensions(list<base_unit> l) {
Chris@16: //     set<base_dimension> dimensions;
Chris@16: //     for_each(base_unit unit : l) {
Chris@16: //         for_each(dim d : unit.dimension_type) {
Chris@16: //             dimensions = insert(dimensions, d.tag_type);
Chris@16: //         }
Chris@16: //     }
Chris@16: //     return(sort(dimensions, _1 > _2, front_inserter(list<base_dimension>())));
Chris@16: // }
Chris@16: 
Chris@16: typedef char set_no;
Chris@16: struct set_yes { set_no dummy[2]; };
Chris@16: 
Chris@16: template<class T>
Chris@16: struct wrap {};
Chris@16: 
Chris@16: struct set_end {
Chris@16:     static set_no lookup(...);
Chris@16:     typedef mpl::long_<0> size;
Chris@16: };
Chris@16: 
Chris@16: template<class T, class Next>
Chris@16: struct set : Next {
Chris@16:     using Next::lookup;
Chris@16:     static set_yes lookup(wrap<T>*);
Chris@16:     typedef T item;
Chris@16:     typedef Next next;
Chris@16:     typedef typename mpl::next<typename Next::size>::type size;
Chris@16: };
Chris@16: 
Chris@16: template<bool has_key>
Chris@16: struct set_insert;
Chris@16: 
Chris@16: template<>
Chris@16: struct set_insert<true> {
Chris@16:     template<class Set, class T>
Chris@16:     struct apply {
Chris@16:         typedef Set type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct set_insert<false> {
Chris@16:     template<class Set, class T>
Chris@16:     struct apply {
Chris@16:         typedef set<T, Set> type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<class Set, class T>
Chris@16: struct has_key {
Chris@16:     static const long size = sizeof(Set::lookup((wrap<T>*)0));
Chris@16:     static const bool value = (size == sizeof(set_yes));
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct find_base_dimensions_impl_impl {
Chris@16:     template<class Begin, class S>
Chris@16:     struct apply {
Chris@16:         typedef typename find_base_dimensions_impl_impl<N-1>::template apply<
Chris@16:             typename Begin::next,
Chris@16:             S
Chris@16:         >::type next;
Chris@16: 
Chris@16:         typedef typename set_insert<
Chris@16:             (has_key<next, typename Begin::item::tag_type>::value)
Chris@16:         >::template apply<
Chris@16:             next,
Chris@16:             typename Begin::item::tag_type
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct find_base_dimensions_impl_impl<0> {
Chris@16:     template<class Begin, class S>
Chris@16:     struct apply {
Chris@16:         typedef S type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct find_base_dimensions_impl {
Chris@16:     template<class Begin>
Chris@16:     struct apply {
Chris@16:         typedef typename find_base_dimensions_impl_impl<(Begin::item::dimension_type::size::value)>::template apply<
Chris@16:             typename Begin::item::dimension_type,
Chris@16:             typename find_base_dimensions_impl<N-1>::template apply<typename Begin::next>::type
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct find_base_dimensions_impl<0> {
Chris@16:     template<class Begin>
Chris@16:     struct apply {
Chris@16:         typedef set_end type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<class T>
Chris@16: struct find_base_dimensions {
Chris@16:     typedef typename insertion_sort<
Chris@16:         typename find_base_dimensions_impl<
Chris@16:             (T::size::value)
Chris@16:         >::template apply<T>::type
Chris@16:     >::type type;
Chris@16: };
Chris@16: 
Chris@16: // calculate_base_dimension_coefficients finds
Chris@16: // the coefficients corresponding to the first
Chris@16: // base_dimension in each of the dimension_lists.
Chris@16: // It returns two values.  The first result
Chris@16: // is a list of the coefficients.  The second
Chris@16: // is a list with all the incremented iterators.
Chris@16: // When we encounter a base_dimension that is
Chris@16: // missing from a dimension_list, we do not
Chris@16: // increment the iterator and we set the
Chris@16: // coefficient to zero.
Chris@16: 
Chris@16: template<bool has_dimension>
Chris@16: struct calculate_base_dimension_coefficients_func;
Chris@16: 
Chris@16: template<>
Chris@16: struct calculate_base_dimension_coefficients_func<true> {
Chris@16:     template<class T>
Chris@16:     struct apply {
Chris@16:         typedef typename T::item::value_type type;
Chris@16:         typedef typename T::next next;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct calculate_base_dimension_coefficients_func<false> {
Chris@16:     template<class T>
Chris@16:     struct apply {
Chris@16:         typedef static_rational<0> type;
Chris@16:         typedef T next;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // begins_with_dimension returns true iff its first
Chris@16: // parameter is a valid iterator which yields its
Chris@16: // second parameter when dereferenced.
Chris@16: 
Chris@16: template<class Iterator>
Chris@16: struct begins_with_dimension {
Chris@16:     template<class Dim>
Chris@16:     struct apply : 
Chris@16:         boost::is_same<
Chris@16:             Dim,
Chris@16:             typename Iterator::item::tag_type
Chris@16:         > {};
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct begins_with_dimension<dimensionless_type> {
Chris@16:     template<class Dim>
Chris@16:     struct apply : mpl::false_ {};
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct calculate_base_dimension_coefficients_impl {
Chris@16:     template<class BaseUnitDimensions,class Dim,class T>
Chris@16:     struct apply {
Chris@16:         typedef typename calculate_base_dimension_coefficients_func<
Chris@16:             begins_with_dimension<typename BaseUnitDimensions::item>::template apply<
Chris@16:                 Dim
Chris@16:             >::value
Chris@16:         >::template apply<
Chris@16:             typename BaseUnitDimensions::item
Chris@16:         > result;
Chris@16:         typedef typename calculate_base_dimension_coefficients_impl<N-1>::template apply<
Chris@16:             typename BaseUnitDimensions::next,
Chris@16:             Dim,
Chris@16:             list<typename result::type, T>
Chris@16:         > next_;
Chris@16:         typedef typename next_::type type;
Chris@16:         typedef list<typename result::next, typename next_::next> next;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct calculate_base_dimension_coefficients_impl<0> {
Chris@16:     template<class Begin, class BaseUnitDimensions, class T>
Chris@16:     struct apply {
Chris@16:         typedef T type;
Chris@16:         typedef dimensionless_type next;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // add_zeroes pushs N zeroes onto the
Chris@16: // front of a list.
Chris@16: //
Chris@16: // list<rational> add_zeroes(list<rational> l, int N) {
Chris@16: //     if(N == 0) {
Chris@16: //         return(l);
Chris@16: //     } else {
Chris@16: //         return(push_front(add_zeroes(l, N-1), 0));
Chris@16: //     }
Chris@16: // }
Chris@16: 
Chris@16: template<int N>
Chris@16: struct add_zeroes_impl {
Chris@16:     // If you get an error here and your base units are
Chris@16:     // in fact linearly independent, please report it.
Chris@16:     BOOST_MPL_ASSERT_MSG((N > 0), base_units_are_probably_not_linearly_independent, (void));
Chris@16:     template<class T>
Chris@16:     struct apply {
Chris@16:         typedef list<
Chris@16:             static_rational<0>,
Chris@16:             typename add_zeroes_impl<N-1>::template apply<T>::type
Chris@16:         > type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct add_zeroes_impl<0> {
Chris@16:     template<class T>
Chris@16:     struct apply {
Chris@16:         typedef T type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // expand_dimensions finds the exponents of
Chris@16: // a set of dimensions in a dimension_list.
Chris@16: // the second parameter is assumed to be
Chris@16: // a superset of the base_dimensions of
Chris@16: // the first parameter.
Chris@16: //
Chris@16: // list<rational> expand_dimensions(dimension_list, list<base_dimension>);
Chris@16: 
Chris@16: template<int N>
Chris@16: struct expand_dimensions {
Chris@16:     template<class Begin, class DimensionIterator>
Chris@16:     struct apply {
Chris@16:         typedef typename calculate_base_dimension_coefficients_func<
Chris@16:             begins_with_dimension<DimensionIterator>::template apply<typename Begin::item>::value
Chris@16:         >::template apply<DimensionIterator> result;
Chris@16:         typedef list<
Chris@16:             typename result::type,
Chris@16:             typename expand_dimensions<N-1>::template apply<typename Begin::next, typename result::next>::type
Chris@16:         > type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct expand_dimensions<0> {
Chris@16:     template<class Begin, class DimensionIterator>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct create_unit_matrix {
Chris@16:     template<class Begin, class Dimensions>
Chris@16:     struct apply {
Chris@16:         typedef typename create_unit_matrix<N - 1>::template apply<typename Begin::next, Dimensions>::type next;
Chris@16:         typedef list<typename expand_dimensions<Dimensions::size::value>::template apply<Dimensions, typename Begin::item::dimension_type>::type, next> type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct create_unit_matrix<0> {
Chris@16:     template<class Begin, class Dimensions>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<class T>
Chris@16: struct normalize_units {
Chris@16:     typedef typename find_base_dimensions<T>::type dimensions;
Chris@16:     typedef typename create_unit_matrix<(T::size::value)>::template apply<
Chris@16:         T,
Chris@16:         dimensions
Chris@16:     >::type matrix;
Chris@16:     typedef typename make_square_and_invert<matrix>::type type;
Chris@16:     static const long extra = (type::size::value) - (T::size::value);
Chris@16: };
Chris@16: 
Chris@16: // multiply_add_units computes M x V
Chris@16: // where M is a matrix and V is a horizontal
Chris@16: // vector
Chris@16: //
Chris@16: // list<rational> multiply_add_units(list<list<rational> >, list<rational>);
Chris@16: 
Chris@16: template<int N>
Chris@16: struct multiply_add_units_impl {
Chris@16:     template<class Begin1, class Begin2 ,class X>
Chris@16:     struct apply {
Chris@16:         typedef list<
Chris@16:             typename mpl::plus<
Chris@16:                 typename mpl::times<
Chris@16:                     typename Begin2::item,
Chris@16:                     X
Chris@16:                 >::type,
Chris@16:                 typename Begin1::item
Chris@16:             >::type,
Chris@16:             typename multiply_add_units_impl<N-1>::template apply<
Chris@16:                 typename Begin1::next,
Chris@16:                 typename Begin2::next,
Chris@16:                 X
Chris@16:             >::type
Chris@16:         > type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct multiply_add_units_impl<0> {
Chris@16:     template<class Begin1, class Begin2 ,class X>
Chris@16:     struct apply {
Chris@16:         typedef dimensionless_type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct multiply_add_units {
Chris@16:     template<class Begin1, class Begin2>
Chris@16:     struct apply {
Chris@16:         typedef typename multiply_add_units_impl<
Chris@16:             (Begin2::item::size::value)
Chris@16:         >::template apply<
Chris@16:             typename multiply_add_units<N-1>::template apply<
Chris@16:                 typename Begin1::next,
Chris@16:                 typename Begin2::next
Chris@16:             >::type,
Chris@16:             typename Begin2::item,
Chris@16:             typename Begin1::item
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct multiply_add_units<1> {
Chris@16:     template<class Begin1, class Begin2>
Chris@16:     struct apply {
Chris@16:         typedef typename add_zeroes_impl<
Chris@16:             (Begin2::item::size::value)
Chris@16:         >::template apply<dimensionless_type>::type type1;
Chris@16:         typedef typename multiply_add_units_impl<
Chris@16:             (Begin2::item::size::value)
Chris@16:         >::template apply<
Chris@16:             type1,
Chris@16:             typename Begin2::item,
Chris@16:             typename Begin1::item
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: // strip_zeroes erases the first N elements of a list if
Chris@16: // they are all zero, otherwise returns inconsistent
Chris@16: //
Chris@16: // list strip_zeroes(list l, int N) {
Chris@16: //     if(N == 0) {
Chris@16: //         return(l);
Chris@16: //     } else if(l.front == 0) {
Chris@16: //         return(strip_zeroes(pop_front(l), N-1));
Chris@16: //     } else {
Chris@16: //         return(inconsistent);
Chris@16: //     }
Chris@16: // }
Chris@16: 
Chris@16: template<int N>
Chris@16: struct strip_zeroes_impl;
Chris@16: 
Chris@16: template<class T>
Chris@16: struct strip_zeroes_func {
Chris@16:     template<class L, int N>
Chris@16:     struct apply {
Chris@16:         typedef inconsistent type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct strip_zeroes_func<static_rational<0> > {
Chris@16:     template<class L, int N>
Chris@16:     struct apply {
Chris@16:         typedef typename strip_zeroes_impl<N-1>::template apply<typename L::next>::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct strip_zeroes_impl {
Chris@16:     template<class T>
Chris@16:     struct apply {
Chris@16:         typedef typename strip_zeroes_func<typename T::item>::template apply<T, N>::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct strip_zeroes_impl<0> {
Chris@16:     template<class T>
Chris@16:     struct apply {
Chris@16:         typedef T type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: // Given a list of base_units, computes the
Chris@16: // exponents of each base unit for a given
Chris@16: // dimension.
Chris@16: //
Chris@16: // list<rational> calculate_base_unit_exponents(list<base_unit> units, dimension_list dimensions);
Chris@16: 
Chris@16: template<class T>
Chris@16: struct is_base_dimension_unit {
Chris@16:     typedef mpl::false_ type;
Chris@16:     typedef void base_dimension_type;
Chris@16: };
Chris@16: template<class T>
Chris@16: struct is_base_dimension_unit<list<dim<T, static_rational<1> >, dimensionless_type> > {
Chris@16:     typedef mpl::true_ type;
Chris@16:     typedef T base_dimension_type;
Chris@16: };
Chris@16: 
Chris@16: template<int N>
Chris@16: struct is_simple_system_impl {
Chris@16:     template<class Begin, class Prev>
Chris@16:     struct apply {
Chris@16:         typedef is_base_dimension_unit<typename Begin::item::dimension_type> test;
Chris@16:         typedef mpl::and_<
Chris@16:             typename test::type,
Chris@16:             mpl::less<Prev, typename test::base_dimension_type>,
Chris@16:             typename is_simple_system_impl<N-1>::template apply<
Chris@16:                 typename Begin::next,
Chris@16:                 typename test::base_dimension_type
Chris@16:             >
Chris@16:         > type;
Chris@16:         static const bool value = (type::value);
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct is_simple_system_impl<0> {
Chris@16:     template<class Begin, class Prev>
Chris@16:     struct apply : mpl::true_ {
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<class T>
Chris@16: struct is_simple_system {
Chris@16:     typedef T Begin;
Chris@16:     typedef is_base_dimension_unit<typename Begin::item::dimension_type> test;
Chris@16:     typedef typename mpl::and_<
Chris@16:         typename test::type,
Chris@16:         typename is_simple_system_impl<
Chris@16:             T::size::value - 1
Chris@16:         >::template apply<
Chris@16:             typename Begin::next::type,
Chris@16:             typename test::base_dimension_type
Chris@16:         >
Chris@16:     >::type type;
Chris@16:     static const bool value = type::value;
Chris@16: };
Chris@16: 
Chris@16: template<bool>
Chris@16: struct calculate_base_unit_exponents_impl;
Chris@16: 
Chris@16: template<>
Chris@16: struct calculate_base_unit_exponents_impl<true> {
Chris@16:     template<class T, class Dimensions>
Chris@16:     struct apply {
Chris@16:         typedef typename expand_dimensions<(T::size::value)>::template apply<
Chris@16:             typename find_base_dimensions<T>::type,
Chris@16:             Dimensions
Chris@16:         >::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<>
Chris@16: struct calculate_base_unit_exponents_impl<false> {
Chris@16:     template<class T, class Dimensions>
Chris@16:     struct apply {
Chris@16:         // find the units that correspond to each base dimension
Chris@16:         typedef normalize_units<T> base_solutions;
Chris@16:         // pad the dimension with zeroes so it can just be a
Chris@16:         // list of numbers, making the multiplication easy
Chris@16:         // e.g. if the arguments are list<pound, foot> and
Chris@16:         // list<mass,time^-2> then this step will
Chris@16:         // yield list<0,1,-2>
Chris@16:         typedef typename expand_dimensions<(base_solutions::dimensions::size::value)>::template apply<
Chris@16:             typename base_solutions::dimensions,
Chris@16:             Dimensions
Chris@16:         >::type dimensions;
Chris@16:         // take the unit corresponding to each base unit
Chris@16:         // multiply each of its exponents by the exponent
Chris@16:         // of the base_dimension in the result and sum.
Chris@16:         typedef typename multiply_add_units<dimensions::size::value>::template apply<
Chris@16:             dimensions,
Chris@16:             typename base_solutions::type
Chris@16:         >::type units;
Chris@16:         // Now, verify that the dummy units really
Chris@16:         // cancel out and remove them.
Chris@16:         typedef typename strip_zeroes_impl<base_solutions::extra>::template apply<units>::type type;
Chris@16:     };
Chris@16: };
Chris@16: 
Chris@16: template<class T, class Dimensions>
Chris@16: struct calculate_base_unit_exponents {
Chris@16:     typedef typename calculate_base_unit_exponents_impl<is_simple_system<T>::value>::template apply<T, Dimensions>::type type;
Chris@16: };
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: } // namespace units
Chris@16: 
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif