链接

[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);
}

分类: 文章

3 条评论

litble · 2018年2月3日 11:49 上午

STO cai Orz

konnyakuxzy · 2018年2月3日 11:30 上午

STO cai Orz

发表评论

电子邮件地址不会被公开。 必填项已用*标注

你是机器人吗? =。= *