test/src/number/enumerate_quotient/divide_both.test.cpp

View this file on GitHub
Last update: 2023-05-11 13:32:09+09:00
Problem: https://atcoder.jp/contests/abc206/tasks/abc206_e

Depends on

商列挙 (library/number/enumerate_quotient.hpp)

Code

#define PROBLEM "https://atcoder.jp/contests/abc206/tasks/abc206_e"

#include <iostream>
#include <unordered_map>

#include "library/number/enumerate_quotient.hpp"

using namespace suisen;

constexpr int MAX = 1000010;

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    int l, r;
    std::cin >> l >> r;
    l = std::max(l - 1, 1);
    std::unordered_map<long long, long long> memo;
    auto f = [&memo](auto self, int l, int r) -> long long {
        const long long key = (long long) l * MAX + r;
        {
            auto it = memo.find(key);
            if (it != memo.end()) return it->second;
        }
        long long res = (long long) (r - l) * (r - l);
        for (auto [lg, rg, qs] : multiple_quotients_ranges(std::array<int, 2>{l, r})) {
            if (lg == 1 and rg == 1) continue;
            auto [nl, nr] = qs;
            res -= (rg - lg + 1) * self(self, nl, nr);
        }
        return memo[key] = res;
    };
    long long whole_pairs = (long long) (r - l) * (r - l);
    long long divisor_pairs = 0;
    for (auto [lg, rg, q] : quotient_ranges(r)) {
        divisor_pairs += (long long) std::max(0, (rg - std::max(lg - 1, l))) * q;
    }
    divisor_pairs += divisor_pairs - (r - l);
    std::cout << whole_pairs - f(f, l, r) - divisor_pairs << std::endl;
    return 0;
}

#line 1 "test/src/number/enumerate_quotient/divide_both.test.cpp"
#define PROBLEM "https://atcoder.jp/contests/abc206/tasks/abc206_e"

#include <iostream>
#include <unordered_map>

#line 1 "library/number/enumerate_quotient.hpp"



#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <tuple>
#include <vector>

namespace suisen {
    /**
     * @return #{ d>0 | exists an integer x s.t. floor(n/x)=d }
    */
    template <typename T, std::enable_if_t<std::is_integral_v<T>, std::nullptr_t> = nullptr>
    T quotient_num(const T& n) {
        const T h = sqrtl(n);
        return 2 * h - (h * (h + 1) > n);
    }

    /**
     * @brief Runs f(q) for each q>0 s.t. (there exists an integer x s.t. floor(n/x)=q) in ascending order.
     * @tparam T integer type
     */
    template <typename T, typename Fun,
        std::enable_if_t<std::conjunction_v<std::is_integral<T>, std::is_invocable<Fun, T>>, std::nullptr_t> = nullptr>
    void foreach_quotients(const T &n, Fun &&f) {
        assert(0 <= n);
        if (n == 0) return;
        const int h = sqrtl(n), m = 2 * h - (h == n / h);
        for (int l = 1; l <= h; ++l) f(l);
        if (n < 1LL << 53) { // the size of the fraction part of double is 52=53-1.
            for (int l = h - (m & 1); l >= 1; --l) f(static_cast<T>(static_cast<double>(n) / l));
        } else {
            for (int l = h - (m & 1); l >= 1; --l) f(n / l);
        }
    }
    /**
     * @brief Calculates an ascending vector of S = { d>0 | exists an integer x s.t. floor(n / x) = d } in O(√n) time.
     * @tparam T integer type
     * @return an ascending vector S = { d>0 | exists an integer x s.t. floor(n / x) = d }.
     */
    template <typename T, std::enable_if_t<std::is_integral_v<T>, std::nullptr_t> = nullptr>
    std::vector<T> enumerate_quotients(T n) {
        std::vector<T> q(quotient_num(n));
        auto it = q.begin();
        foreach_quotients(n, [&it](const T& v) { *it++ = v; });
        return q;
    }

    /**
     * @brief Calculates vector of { l = min S, r = max S, q }, where S = { d | floor(n / d) = q } in O(√n) time.
     * @tparam T integer type
     * @return vector of { l = min S, r = max S, q }, where S = { d | floor(n / d) = q }.
     * It is guaranteed that l, r is ascending and q is descending.
     */
    template <typename T, std::enable_if_t<std::is_integral_v<T>, std::nullptr_t> = nullptr>
    std::vector<std::tuple<T, T, T>> quotient_ranges(T n) {
        assert(0 <= n);
        const std::vector<T> qs = enumerate_quotients(n);
        const int m = qs.size();
        std::vector<std::tuple<T, T, T>> res(qs.size());
        for (int i = 0; i < m; ++i) {
            T l = i ? qs[i - 1] + 1 : 1, r = qs[i], q = qs[m - i - 1];
            res[i] = std::make_tuple(l, r, q);
        }
        return res;
    }

    /**
     * @brief Calculates vector of { l = min S, r = max S, (q[0], ..., q[k-1]) }, where S = { d | floor(vs[i] / d) = q[i] for i = 0,...,k-1 } in O(|vs| * Σ√v) time.
     * @tparam Container vector or array
     * @return vector of { l = min S, r = max S, (q[0], ..., q[k-1]) }, where S is { d | floor(vs[i] / d) = q[i] for i = 0,...,k-1 }.
     * It is guaranteed that l, r is ascending and q[i] is descending. 
     */
    template <typename Container, std::enable_if_t<std::is_integral_v<typename Container::value_type>, std::nullptr_t> = nullptr>
    std::vector<std::tuple<typename Container::value_type, typename Container::value_type, Container>> multiple_quotients_ranges(const Container& vs) {
        using T = typename Container::value_type;
        int n = vs.size();
        T max_val = *std::max_element(vs.begin(), vs.end());
        assert(*std::min_element(vs.begin(), vs.end()) >= 0);
        std::vector<std::tuple<T, T, Container>> res;
        for (T l = 1, r = 1; l <= max_val; l = r + 1) {
            Container qs{};
            if constexpr (std::is_same_v<Container, std::vector<T>>) qs.resize(n);
            r = std::numeric_limits<T>::max();
            for (int i = 0; i < n; ++i) {
                qs[i] = vs[i] / l;
                r = std::min(r, qs[i] == 0 ? std::numeric_limits<T>::max() : vs[i] / qs[i]);
            }
            res.emplace_back(l, r, std::move(qs));
        }
        return res;
    }
} // namespace suisen



#line 7 "test/src/number/enumerate_quotient/divide_both.test.cpp"

using namespace suisen;

constexpr int MAX = 1000010;

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    int l, r;
    std::cin >> l >> r;
    l = std::max(l - 1, 1);
    std::unordered_map<long long, long long> memo;
    auto f = [&memo](auto self, int l, int r) -> long long {
        const long long key = (long long) l * MAX + r;
        {
            auto it = memo.find(key);
            if (it != memo.end()) return it->second;
        }
        long long res = (long long) (r - l) * (r - l);
        for (auto [lg, rg, qs] : multiple_quotients_ranges(std::array<int, 2>{l, r})) {
            if (lg == 1 and rg == 1) continue;
            auto [nl, nr] = qs;
            res -= (rg - lg + 1) * self(self, nl, nr);
        }
        return memo[key] = res;
    };
    long long whole_pairs = (long long) (r - l) * (r - l);
    long long divisor_pairs = 0;
    for (auto [lg, rg, q] : quotient_ranges(r)) {
        divisor_pairs += (long long) std::max(0, (rg - std::max(lg - 1, l))) * q;
    }
    divisor_pairs += divisor_pairs - (r - l);
    std::cout << whole_pairs - f(f, l, r) - divisor_pairs << std::endl;
    return 0;
}