2 // Copyright 2013 Francisco Jerez
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:
11 // The above copyright notice and this permission notice shall be included in
12 // all copies or substantial portions of the Software.
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.
23 #ifndef CLOVER_UTIL_FACTOR_HPP
24 #define CLOVER_UTIL_FACTOR_HPP
26 #include "util/range.hpp"
31 /// Calculate all prime integer factors of \p x.
33 /// If \p limit is non-zero, terminate early as soon as enough
34 /// factors have been collected to reach the product \p limit.
38 find_integer_prime_factors(T x, T limit = 0)
40 const T max_d = (limit > 0 && limit < x ? limit : x);
41 const T min_x = x / max_d;
42 std::vector<T> factors;
44 for (T d = 2; d <= max_d && x > min_x; d++) {
46 for (; x % d == 0; x /= d);
56 /// Walk the power set of prime factors of the n-dimensional
57 /// integer array \p grid subject to the constraints given by
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()) {
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);
79 // Try the next power of this factor.
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.
90 // Overflow, back off to the zeroth power,
91 while (block.second[d] % f == 0) {
96 // ...and carry to the next factor.
97 return next_grid_factor(limits, grid, factors,
98 std::move(block), d, i + 1);
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.
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)));
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)