Home | History | Annotate | Download | only in util 
      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