Bài này search ra thấy có đề cập ở đây: https://stackoverflow.com/questions/26463250/find-the-sum-of-least-common-multiples-of-all-subsets-of-a-given-set Searhch thêm nữa thì ra Project Euler Problem 858 cũng khá giống (khác 1 chỗ LCM của tập rỗng = 1). Nó là dãy OEIS A226037
Duyệt trâu thì mất O(n2^n) với n là kích thước của mảng input, n \leq 100 nên duyệt trâu có thể lên tới 100 * 2^{100} là bất khả thi.
Cái gợi ý ở trang Stack Overflow kia là đại khái như thế này: vì giá trị của a_i \leq 500 nên nếu bỏ các số khi mà phân tích thừa số nguyên tố (tsnt) ra có số nguyên tố “lớn” hơn 19 thì bỏ qua 1 bên tính sau. Các số còn lại có tsnt nhỏ (2, 3, 5, 7, 11, 13, 17, 19) thì tuy có thể lên tới 2^m tập hợp con với m \leq n nhưng giá trị BCNN (LCM - least common multiple) của các tập hợp này chỉ có thể có tối đa tầm 70000 giá trị, vì khi phân tích tsnt ra rồi lấy LCM thì LCM của x số có số mũ là max của các số mũ của tsnt, mà 1 số có tsnt 2 thì tối đa chỉ có số mũ là 2^8, 3 là 3^5, , 5 là 5^3, 3 là 7^3, 11, 13, 17, 19 có số mũ tối đa là 2, vậy giá trị LCM có thể là 2^{0..8}3^{0..5}5^{0..3}7^{0..3}11^{0..2}13^{0..2}17^{0..2}19^{0..2} hay chỉ có thể có 9*6*4*4*3*3*3*3 = 69984 giá trị. Vì vậy làm 1 cái mảng dp có 69984 phần tử, dp[i] là số lượng số tập hợp có LCM = i.
Ví dụ mảng arr = {2, 4, 5} thì các tập hợp con của arr là \{\}, \{2\}, \{4\}, \{5\}, \{2,4\}, \{2,5\}, \{4,5\}, \{2,4,5\} có LCM lần lượt là 1, 2, 4, 5, 4, 10, 20, 20 thì dp có thể là 1 hash map (cho dễ hình dung) là {1: 1, 2: 1, 4: 2, 5: 1, 10: 1, 20: 2}. Lưu ý ở đây tính LCM của tập rỗng là 1.
Để thêm 1 phần tử mới có giá trị tạm gọi là new_value vào hash map dp đã tính rồi thì có thể tính đơn giản duyệt trâu, tối đa chỉ lặp 70000 lần:
new_dp = copy(dp)
for k, v in dp:
new_set_lcm = LCM(k, new_value)
new_dp[new_set_lcm] += v
swap(dp, new_dp)
Ví dụ thêm 10 và arr = {2, 4, 5} thì nó chạy lần lượt là:
new_dp = copy(dp)
k = 1, v = 1: LCM(1, 10) = 10, vậy new_dp[10] += 1
k = 2, v = 1: LCM(2, 10) = 10, vậy new_dp[10] += 1
k = 4, v = 2: LCM(4, 10) = 20, vậy new_dp[20] += 2
k = 5, v = 1: LCM(5, 10) = 10, vậy new_dp[10] += 1
k = 10, v = 1: LCM(10, 10) = 10, vậy new_dp[10] += 1
k = 20, v = 2: LCM(20, 10) = 20, vậy new_dp[20] += 2
--> new_dp = {1: 1, 2: 1, 4: 2, 5: 1, 10: 5, 20: 6}
Vì thuật toán duyệt trâu thế này nên cho LCM tập rỗng là 1 để tiện đường tính LCM(1, new_value).
Tổng số values() của hash map dp chính là số tập con của mảng n phần tử và bằng 2^n. Ở đây ta thấy khi thêm 1 phần tử vào mảng thành mảng n + 1 phần tử thì thay vì phải tính LCM 2^{n-1} lần ta chỉ tính tối đa 70000 lần, làm 100 lần thì cũng chỉ tốn 7 * 10^6 lần tính LCM, khả thi.
Nhưng mảng arr còn có thể chứa phần tử có tsnt “lớn” như 23, 29, ..., 499 (số mũ của các tsnt lớn này luôn là 1) thì trong cái solution sketch ở cái competition tút tận năm 2012 kia ko đề cập, chắc là “quá dễ, để độc giả tự tìm cách” 
Tôi bí 2 ngày cuối tuần rốt cuộc cũng mò ra: đầu tiên là nhóm các số có tsnt lớn lại theo tsnt đó, và chia số đó cho tsnt đó ra. Ví dụ 46, 23, 276, 93, 310 thì nhóm lại thành {23: [2, 1, 12], 31: [3, 10]} rồi lần lượt “đếm” và “gom” các nhóm này vào dp. Khi “gom” vào xong thì values() trong dp ko còn là 2^n nữa mà là 1 cái tổng kì dị nào đó
Cách “đếm” ví dụ cho tsnt lớn là
dp_lg = {} // chỉ khởi tạo khi thêm số có tsnt = 23 đầu tiên (lg là viết tắt của large)
new_dp_lg = copy(dp_lg)
for k in <dp.keys() giao dp_lg.keys()>:
new_set_lcm = LCM(k, new_value_without_large_prime)
new_dp_lg[new_set_lcm] += dp[k] + dp_lg[k]
swap(dp_lg, new_dp_lg)
Ở đây ta tạo 1 mảng mới dp_lg rồi thêm vào tương tự như số có tsnt nhỏ, nhưng thay vì chỉ duyệt keys là k trong dp ta còn phải duyệt keys trong dp_lg nữa, vì các số mới có thể có tsnt mới so với các tsnt trong dp. Ngoài ra ta còn phải lấy LCM của new_value đã chia cho tsnt lớn, ví dụ 46 thì lấy LCM của 2 thôi chứ ko phải LCM của 46, và khi cộng vào thì ta cộng thêm dp_lg[k] vào nữa, ko chỉ cộng có mỗi dp[k].
Cộng hết các phần tử mới có tsnt lớn là p rồi thì ta “gom” dp_lg lại vào dp bằng cách:
for k in <dp.keys() giao dp_lg.keys()>:
dp[k] += v * p
Rồi tiếp tục tạo dp_lg = {} mới cho nhóm có tsnt lớn mới.
Toy viết ra 2 maps dp và dp_lg vậy để dễ hình dung thôi chứ khi implement toy chỉ xài có 1 map (có value type là pair<Int, Int>) để tránh cái chỗ tìm tập hợp keys giao của 2 maps. Mà viết ra tới đây thì nghĩ ra cách tránh tìm giao mới bằng lập 2 vòng mỗi vòng cho 1 map và cộng dp[k] hoặc dp_lg[k] chứ ko cộng cả 2
Thôi dù sao cũng đã code rồi, code của toy đây: https://godbolt.org/z/cT3fMKcab
#include <iostream>
#include <algorithm>
#include <vector>
#include <array>
#include <cmath>
#include <numeric>
#include <unordered_map>
#include <fmt/format.h>
#include <fmt/ranges.h>
#include <chrono>
namespace cron = std::chrono;
struct IntMod {
static constexpr int kMod = 1'000'000'007;
int val;
explicit IntMod(int val = 0) : val{val} {}
explicit IntMod(int64_t val) : val{static_cast<int>(val % kMod)} {}
IntMod& operator+=(const IntMod& other) {
val += other.val;
if (val >= kMod) val -= kMod;
return *this;
}
IntMod operator+(const IntMod& other) const {
IntMod res = *this;
res += other;
return res;
}
IntMod& operator*=(const IntMod& other) {
val = (static_cast<int64_t>(val) * other.val) % kMod;
return *this;
}
IntMod operator*(const IntMod& other) const {
IntMod res = *this;
res *= other;
return res;
}
};
consteval bool isSmallPrime(int n) {
for (int i = 2; i * i <= n; ++i)
if (n % i == 0) return false;
return true;
}
consteval unsigned countSmallPrimesLeqSqrt(int n) {
unsigned res{};
for (int i = 2; i * i <= n; ++i)
if (isSmallPrime(i)) ++res;
return res;
}
template <unsigned N>
consteval auto getSmallPrimesLeqSqrt() {
std::array<int, countSmallPrimesLeqSqrt(N)> res{};
auto out = begin(res);
for (int i = 2; i * i <= N; ++i)
if (isSmallPrime(i)) *out++ = i;
return res;
}
consteval int logBase(int n, int base) {
// Find largest e such that base^e <= n
int res{};
int m = base;
for (; m <= n; m *= base) ++res;
return res - static_cast<int>(m == n);
}
template <unsigned N>
consteval auto getSmallPrimeMaxExponents() {
std::array<int, countSmallPrimesLeqSqrt(N)> res{};
auto out = begin(res);
for (int p : getSmallPrimesLeqSqrt<N>()) *out++ = logBase(N, p);
return res;
}
template <unsigned NMax>
struct SmallPrimeFactorizer {
// Optimization: user must provide arrays of primes, maxExponents
static constexpr auto kPrimes = getSmallPrimesLeqSqrt<NMax>();
static constexpr auto kMaxExponents = getSmallPrimeMaxExponents<NMax>();
using ExponentType = std::array<int, kPrimes.size()>;
static std::pair<ExponentType, int> factorize(int n) {
ExponentType ex{};
size_t i{};
for (int p : kPrimes) {
for (; n % p == 0; n /= p) ex[i]++;
++i;
}
return {ex, n};
}
static constexpr int getExponentsMaxArraySize() {
return std::accumulate(begin(kMaxExponents), end(kMaxExponents), 1, [](int s, int n) { return s * (n + 1); });
}
static int toInt(const ExponentType& e) {
int res{};
int baseMult{1};
for (size_t i = 0; i < kMaxExponents.size(); ++i) {
res += e[i] * baseMult;
baseMult *= kMaxExponents[i] + 1;
}
return res;
}
static ExponentType fromInt(int id) {
ExponentType res{};
int baseMult{1};
for (size_t i = 0; i < kMaxExponents.size(); ++i) {
baseMult = kMaxExponents[i] + 1;
res[i] = id % baseMult;
id /= baseMult;
}
return res;
}
static IntMod toValue(const ExponentType& e) {
IntMod res{1};
for (size_t i = 0; i < kPrimes.size(); ++i) {
const int v = static_cast<int>(std::pow(kPrimes[i], e[i]));
res *= IntMod{v};
}
return res;
}
static ExponentType lcm(const ExponentType& lhs, const ExponentType& rhs) {
ExponentType res{};
for (size_t i = 0; i < res.size(); ++i) res[i] = std::max(lhs[i], rhs[i]);
return res;
}
static int sumLcmSubsets(const std::vector<int>& arr) {
// Separate numbers that have large prime factors (and factorize all numbers)
std::vector<ExponentType> arrSm;
std::unordered_map<int, std::vector<ExponentType>> mapLg;
for (auto n : arr) {
auto [exponents, prime] = factorize(n);
(prime == 1 ? arrSm : mapLg[prime]).push_back(exponents);
}
// Solve for number with small prime factors
constexpr int kSz = getExponentsMaxArraySize();
std::vector<IntMod> dp(kSz);
dp[0] = IntMod{1};
for (const auto& fn : arrSm) {
auto dp2 = dp;
for (int k = 0; k < kSz; ++k) {
if (dp[k].val == 0) continue;
const auto fx = lcm(fromInt(k), fn);
const auto ix = toInt(fx);
dp2[ix] += dp[k];
}
dp.swap(dp2);
}
// Convert dp<KeyType, int> to dp<KeyType, <int, int>>
std::vector<std::pair<IntMod, IntMod>> dpLg(kSz);
for (int i = 0; i < kSz; ++i) dpLg[i] = std::make_pair(dp[i], IntMod{});
// Solve for number with 1 large prime factor
for (auto& [p, arrLg] : mapLg) {
for (const auto& fs : arrLg) {
auto dp2 = dpLg;
for (int k = 0; k < kSz; ++k) {
const auto fx = lcm(fromInt(k), fs);
dp2[toInt(fx)].second += dpLg[k].first + dpLg[k].second;
}
dpLg.swap(dp2);
}
for (int k = 0; k < kSz; ++k) {
auto& v = dpLg[k];
v.first += v.second * IntMod{p};
v.second = IntMod{};
}
}
return std::accumulate(
begin(dpLg), end(dpLg), IntMod{},
[i = 0](auto sum, const auto& v) mutable { return sum + v.first * toValue(fromInt(i++)); })
.val -
1;
}
};
int main() {
const std::array N{4, 20, 50, 100, 300, 500, /*800*/};
//(+1 for ProjectEuler 858)
const std::array A{87, 108589718, 907004818, 591116486, 962860665, 618948632, 973077198};
for (int i = 0; int n : N) {
std::vector<int> arr(n);
std::iota(begin(arr), end(arr), 1);
fmt::print("arr = {{1..{}}}: ", n);
const auto start = cron::steady_clock::now();
int res = n <= 100 ? SmallPrimeFactorizer<100>::sumLcmSubsets(arr)
: n <= 200 ? SmallPrimeFactorizer<200>::sumLcmSubsets(arr)
: n <= 300 ? SmallPrimeFactorizer<300>::sumLcmSubsets(arr)
: n <= 400 ? SmallPrimeFactorizer<400>::sumLcmSubsets(arr)
: n <= 500 ? SmallPrimeFactorizer<500>::sumLcmSubsets(arr)
: SmallPrimeFactorizer<800>::sumLcmSubsets(arr);
cron::duration<double> elapsed = cron::steady_clock::now() - start;
if (res != A[i++])
fmt::println("X incorrect");
else
fmt::println("{} in {}s", res, elapsed.count());
}
}
Code này để giải Project Euler 858 nên mấy cái test ko có giá trị trùng lặp, ai rảnh thì test chắc chắn cũng đúng thôi Và cái kMod là 1e9 + 7 chứ ko phải 1e5 + 7 như bài này yêu cầu.
Sửa kMod thành 1e5 + 7 và gọi SmallPrimeFactorizer<500>::sumLcmSubsets(arr) là solve đúng hết cho bài chủ thớt hỏi. Để chạy lẹ hơn có thể tìm max(arr) rồi if else so sánh với mấy NMax = 100, 200, 300, 400 . Ví dụ NMax = 100 thì cái mảng primes chỉ có 2, 3, 5, 7 và maxExponents chỉ là 6, 4, 2, 2 hay dp cho số có tsnt nhỏ chỉ có tối đa 7*5*3*3 = 315 phần tử thôi, lẹ hơn tính 70000 phần tử nhiều. Để tính chính xác cho mọi giá trị max(arr) thì có thể xài std::vector thay cho std::array nhưng mà tạo nhiều std::vector quá nó cấp phát động nhiều quá làm chậm chương trình đi rất nhiều, ví dụ tìm đáp số của Project Euler 858 xài std::array chỉ mất 6 giây trên máy toy, còn xài std::vector thì mất 88 giây Khủng hơn thì có thể if else 500 lần thì cũng lẹ nhất có thể mà thôi toy ko bị khùng
83% thành viên diễn đàn không hỏi bài tập, còn bạn thì sao?