BZOJ3309 DZY Loves Math 【莫比乌斯反演】

标签：--   code   nlogn   break   因子   等于   bug   推出   namespace   

题目

对于正整数n，定义f(n)为n所含质因子的最大幂指数。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。
给定正整数a,b，求sigma(sigma(f(gcd(i,j)))) (i=1..a, j=1..b)。

输入格式

第一行一个数T，表示询问数。
接下来T行，每行两个数a,b，表示一个询问。

输出格式

对于每一个询问，输出一行一个非负整数作为回答。

输入样例

4

7558588 9653114

6514903 4451211

7425644 1189442

6335198 4957

输出样例

35793453939901

14225956593420

4332838845846

15400094813

提示

【数据规模】

T<=10000

1<=a,b<=10^7

题解

前面的推导很套路：
\[ans = \sum\limits_{i = 1}^{n} \sum\limits_{j = 1}^{m} f[gcd(i,j)]\]
\[=\sum\limits_{d = 1}^{n} f[d] * \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum\limits_{j = 1}^{\lfloor \frac{m}{d} \rfloor} [gcd(i,j) == 1]\]
\[=\sum\limits_{d = 1}^{n} f[d] * \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(i) * \lfloor \frac{n}{id} \rfloor\lfloor \frac{m}{id} \rfloor\]
\[=\sum\limits_{T = 1}^{n} \lfloor \frac{n}{T} \rfloor\lfloor \frac{m}{T} \rfloor \sum\limits_{d|T} f[d] * \mu(\frac{T}{d})\]
后面那玩意\(g(T) = \sum\limits_{d|T} f[d] * \mu(\frac{T}{d})\)如果能预处理出来，就能\(O(T\sqrt{n})\)计算了
然后我只会\(O(nlogn)\)，，，，
去膜题解

要利用\(\mu(i)\)的性质
显然\(i\)有平方项就不用考虑了
所以\(T = \prod\limits_{i = 1}^{k} p_i^{a_i}\)中每个\(a_i\)最多被取掉\(1\)
设最大为\(r\)
所以\(f(T) = r\)或\(r - 1\)
我们设有\(x\)个这样的指数为\(r\)，那么剩余的\(k - x\)个质因子的指数就可以任选，有\(2^{k - x}\)中选法
如果\(k \ne x\)，\(2^{k - x}\)为偶数，对应\(\mu\)的正负数量相等，最后和为\(0\)
所以只有\(k = x\)时，\(g(T) \ne 0\)
否则我们假使所有的\(f(d)\)都等于\(r\)，那么和依旧为\(0\)，但是实际上当\(k\)个数都被选的时候\(f(d) = r - 1\)，多了一个\(-1\)，根据奇偶性，最后会产生\((-1)^{k + 1}\)的贡献
所以此时\(g(T) = (-1)^{k + 1}\)

具体可以先筛出\(\mu(i)\)，再由\(\mu(i) \ne 0\)的\(i\)推出所有的\(f(i^x)\)，这样做每个数只会被推一次，所以是\(O(n)\)的

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<‘ ‘; puts("");
using namespace std;
const int maxn = 10000005,maxm = 100005,N = 1e7,INF = 1000000000;
inline int read(){
    int out = 0,flag = 1; char c = getchar();
    while (c < 48 || c > 57){if (c == ‘-‘) flag = -1; c = getchar();}
    while (c >= 48 && c <= 57){out = (out << 3) + (out << 1) + c - 48; c = getchar();}
    return out * flag;
}
int p[maxn],pi,isn[maxn],mu[maxn];
LL g[maxn];
void init(){
    mu[1] = 1;
    for (int i = 2; i <= N; i++){
        if (!isn[i]) p[++pi] = i,mu[i] = -1;
        for (int j = 1; j <= pi && i * p[j] <= N; j++){
            isn[i * p[j]] = true;
            if (i % p[j] == 0){
                mu[i * p[j]] = 0;
                break;
            }
            mu[i * p[j]] = -mu[i];
        }
    }
    for (LL i = 2; i <= N; i++)
        if (mu[i] != 0){
            for (LL j = i,t = -mu[i]; j <= N; j *= i)
                g[j] = t;
        }
    for (int i = 1; i <= N; i++) g[i] += g[i - 1];
}
int main(){
    init();
    int T = read(),n,m;
    LL ans;
    while (T--){
        n = read(); m = read(); ans = 0;
        if (n > m) swap(n,m);
        for (int i = 1,nxt; i <= n; i = nxt + 1){
            nxt = min(n / (n / i),m / (m / i));
            ans += 1ll * (n / i) * (m / i) * (g[nxt] - g[i - 1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}

BZOJ3309 DZY Loves Math 【莫比乌斯反演】

标签：--   code   nlogn   break   因子   等于   bug   推出   namespace   

原文地址：https://www.cnblogs.com/Mychael/p/8974560.html