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
28namespace 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