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> find_integer_prime_factors(T x,T limit=0)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>> next_grid_factor(const std::pair<T,std::vector<T>> & limits,const std::vector<T> & grid,const std::vector<std::vector<T>> & factors,std::pair<T,std::vector<T>> block,unsigned d=0,unsigned i=0)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> find_grid_optimal_factor(T product_limit,const std::vector<T> & coord_limit,const std::vector<T> & grid)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