链接

[ProjectEuler466]

题目大意

定义 $P(N, M)$表示在整数 $x\in[1, N], y\in[1, M]$时，$x\times y$能够表示的数的数量

$P(3, 4) = 8$

求 $P(64, 10^{16})$

题解

设 $A_i$表示集合 ${i\times y\ |\ y\in[1, M]}$

则 $$P(N, M) = \left|\bigcup\limits_{k=1}^NA_k\right| = \sum_{S\subset [1, M]\atop {S\neq \emptyset}}(-1)^{|S|-1}
\left|\bigcap\limits_{j\in S}A_j\right|$$

考虑交集的部分，一个数要是所有 $S$内所有数字的倍数当且仅当它是 $\mathrm{lcm}(S)$的倍数，再考虑 $S$中最小的数，它限制了交集中最大的数也只有 $\min(S)\times N$，考虑 $1$到 $\min(S)\times N$中的数，这些数中是 $\mathrm{lcm}(S)$的倍数的数的也一定是交集中的数，所以 $$\left|\bigcap\limits_{j\in S}A_j\right|=\left\lfloor\frac{\min(S)\times M}{\mathrm{lcm}(S)}\right\rfloor$$

设 $$f_{i, j}=\sum_{{S\subset[1, M]} \atop {\min(S)=i \atop {\mathrm{lcm}(S)=j}}}(-1)^{|S|-1}$$

那么 $$P(N, M) = \sum_{i\in[1, M]\atop j=\mathrm{lcm}}\left(\left\lfloor\frac{i\times M}{j}\right\rfloor\times f_{i, j}\right)$$

所以只需要求出 $f$即可

好像已经没有办法优化了，那么我们就用一个 map 来存下所有的状态，然后从小到大枚举元素，每次尝试在集合中加入一个新元素，这样得到了一个新的状态 $(i, j)$，如果发现 $j>i\times N$，那么就剪枝，因为 $j$是单调不减，$i$是单调不增，如果现在满足这个条件也就是 $\lfloor\frac{i\times N}{j}\rfloor=0$，那么由这个状态拓展出的所有状态也会满足等于 0，所以可以剪枝，最后将 $f$值等于 0 的状态减掉，然后发现只要 13s，差不多 1G 的内存

实际上 $f$中不为 0 的状态只有 8085879 个

$P(64, 10^{16}) = 258381958195474745$

代码

#include <map>
#include <cstdio>
#include <cstring>
#include <algorithm>

#define DEBUG(args...) fprintf(stderr, args)

#define FOR(i, l, r) for(int i = (l), i##_end = (r); i <= i##_end; ++i)
#define REP(i, l, r) for(int i = (l), i##_end = (r); i <  i##_end; ++i)
#define DFR(i, l, r) for(int i = (l), i##_end = (r); i >= i##_end; --i)
#define DRP(i, l, r) for(int i = (l), i##_end = (r); i >  i##_end; --i)

typedef long long LL;

template<class T>T Min(const T &a, const T &b) {return a < b ? a : b;}
template<class T>T Max(const T &a, const T &b) {return a > b ? a : b;}
template<class T>bool Chkmin(T &a, const T &b) {return a > b ? a = b, 1 : 0;}
template<class T>bool Chkmax(T &a, const T &b) {return a < b ? a = b, 1 : 0;}

class FastInput {
private:
    static const int SIZE = 1 << 15 | 1;
    char buf[SIZE], *front, *back;

    void Next(char &c) {
        if(front == back) back = (front = buf) + fread(buf, 1, SIZE, stdin);
        c = front == back ? (char)EOF : *front++;
    }

public :
    template<class T>void operator () (T &x) {
        char c, f = 1;
        for(Next(c); c > '9' || c < '0'; Next(c)) if(c == '-') f = -1;
        for(x = 0; c >= '0' && c <= '9'; Next(c)) x = x * 10 + c - '0';
        x *= f;
    }
    void operator () (char &c, char l = 'a', char r = 'z') {
        for(Next(c); c > r || c < l; Next(c)) ;
    }
}read;

const int SN = 100000 + 47;
const int SX = 47;

LL Calc(int, LL);
LL Gcd(LL, LL);

int main() {

    LL n, m;

    n = 64, m = 10000000000000000LL;
    //n = 3, m = 10000000000000000LL;

    printf("%lld\n", Calc(n, m));

    return 0;

}

LL Calc(int x, LL y) {
    int min;
    LL lcm, t, ans = 0;
    std::map<std::pair<int, LL>, LL> f, g;
    std::map<std::pair<int, LL>, LL>::iterator it;
    g.insert(std::make_pair(std::make_pair(x + 1, 1), -1));
    FOR(i, 1, x) {
        f = g; DEBUG("%d\n", i);
        for(it = f.begin(); it != f.end(); ++it) {
            min = Min(i, (it->first).first);
            lcm = (it->first).second;
            t = i / Gcd(i, lcm);
            if(lcm > 1LL * y * min / t) continue;
            g[std::make_pair(min, lcm * t)] -= it->second;
        }
        for(it = g.begin(); it != g.end(); ++it)
            if(!it->second)
                g.erase(it);
    }
    for(it = g.begin(); it != g.end(); ++it)
        if((it->first).first <= x)
            ans += 1LL * y * (it->first).first /
                (it->first).second * it->second;
    DEBUG("%lu\n", g.size());
    return ans;
}

LL Gcd(LL a, LL b) {
    return !b ? a : Gcd(b, a % b);
}