Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. //
  2. // Copyright 2013 Francisco Jerez
  3. //
  4. // Permission is hereby granted, free of charge, to any person obtaining a
  5. // copy of this software and associated documentation files (the "Software"),
  6. // to deal in the Software without restriction, including without limitation
  7. // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  8. // and/or sell copies of the Software, and to permit persons to whom the
  9. // Software is furnished to do so, subject to the following conditions:
  10. //
  11. // The above copyright notice and this permission notice shall be included in
  12. // all copies or substantial portions of the Software.
  13. //
  14. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  15. // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  16. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  17. // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  18. // OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  19. // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  20. // OTHER DEALINGS IN THE SOFTWARE.
  21. //
  22.  
  23. #ifndef CLOVER_UTIL_FACTOR_HPP
  24. #define CLOVER_UTIL_FACTOR_HPP
  25.  
  26. #include "util/range.hpp"
  27.  
  28. namespace clover {
  29.    namespace factor {
  30.       ///
  31.       /// Calculate all prime integer factors of \p x.
  32.       ///
  33.       /// If \p limit is non-zero, terminate early as soon as enough
  34.       /// factors have been collected to reach the product \p limit.
  35.       ///
  36.       template<typename T>
  37.       std::vector<T>
  38.       find_integer_prime_factors(T x, T limit = 0)
  39.       {
  40.          const T max_d = (limit > 0 && limit < x ? limit : x);
  41.          const T min_x = x / max_d;
  42.          std::vector<T> factors;
  43.  
  44.          for (T d = 2; d <= max_d && x > min_x; d++) {
  45.             if (x % d == 0) {
  46.                for (; x % d == 0; x /= d);
  47.                factors.push_back(d);
  48.             }
  49.          }
  50.  
  51.          return factors;
  52.       }
  53.  
  54.       namespace detail {
  55.          ///
  56.          /// Walk the power set of prime factors of the n-dimensional
  57.          /// integer array \p grid subject to the constraints given by
  58.          /// \p limits.
  59.          ///
  60.          template<typename T>
  61.          std::pair<T, std::vector<T>>
  62.          next_grid_factor(const std::pair<T, std::vector<T>> &limits,
  63.                           const std::vector<T> &grid,
  64.                           const std::vector<std::vector<T>> &factors,
  65.                           std::pair<T, std::vector<T>> block,
  66.                           unsigned d = 0, unsigned i = 0) {
  67.             if (d >= factors.size()) {
  68.                // We're done.
  69.                return {};
  70.  
  71.             } else if (i >= factors[d].size()) {
  72.                // We're done with this grid dimension, try the next.
  73.                return next_grid_factor(limits, grid, factors,
  74.                                        std::move(block), d + 1, 0);
  75.  
  76.             } else {
  77.                T f = factors[d][i];
  78.  
  79.                // Try the next power of this factor.
  80.                block.first *= f;
  81.                block.second[d] *= f;
  82.  
  83.                if (block.first <= limits.first &&
  84.                    block.second[d] <= limits.second[d] &&
  85.                    grid[d] % block.second[d] == 0) {
  86.                   // We've found a valid grid divisor.
  87.                   return block;
  88.  
  89.                } else {
  90.                   // Overflow, back off to the zeroth power,
  91.                   while (block.second[d] % f == 0) {
  92.                      block.second[d] /= f;
  93.                      block.first /= f;
  94.                   }
  95.  
  96.                   // ...and carry to the next factor.
  97.                   return next_grid_factor(limits, grid, factors,
  98.                                           std::move(block), d, i + 1);
  99.                }
  100.             }
  101.          }
  102.       }
  103.  
  104.       ///
  105.       /// Find the divisor of the integer array \p grid that gives the
  106.       /// highest possible product not greater than \p product_limit
  107.       /// subject to the constraints given by \p coord_limit.
  108.       ///
  109.       template<typename T>
  110.       std::vector<T>
  111.       find_grid_optimal_factor(T product_limit,
  112.                                const std::vector<T> &coord_limit,
  113.                                const std::vector<T> &grid) {
  114.          const std::vector<std::vector<T>> factors =
  115.             map(find_integer_prime_factors<T>, grid, coord_limit);
  116.          const auto limits = std::make_pair(product_limit, coord_limit);
  117.          auto best = std::make_pair(T(1), std::vector<T>(grid.size(), T(1)));
  118.  
  119.          for (auto block = best;
  120.               block.first != 0 && best.first != product_limit;
  121.               block = detail::next_grid_factor(limits, grid, factors, block)) {
  122.             if (block.first > best.first)
  123.                best = block;
  124.          }
  125.  
  126.          return best.second;
  127.       }
  128.    }
  129. }
  130.  
  131. #endif
  132.