]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |