src/gallium/state_trackers/clover/util/factor.hpp

   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