【BZOJ3309】DZY Loves Math

3309: DZY Loves Math
Time Limit: 20 Sec Memory Limit: 512 MB
Submit: 411 Solved: 161
[Submit][Status][Discuss]
Description

对于正整数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)。

Input

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

Output

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

Sample Input
4

10000

7558588 9653114

6514903 4451211

7425644 1189442

6335198 4957

Sample Output
35793453939901

14225956593420

4332838845846

15400094813
HINT

【数据规模】

T ≤ 10000

1 ≤ a,b ≤ 10 7

Source
为了方便代入反演，我们把原题里那个f函数在这里用g代替
令f(i)为 gcd(x,y)=i,1≤x≤a,1≤y≤b 的数对数目。
然后有 f(i)=∑i|dμ(di)⌊nd⌋⌊md⌋
对于每个询问，我们要求的答案就是

当然就会想到对后面那个和式做前缀和，然后搭配枚举除法使用，把ProblemB里那个莫比乌斯函数的前缀和替换成这里这个和式的前缀和。
要怎么求？
当然我们要先排除掉那些 μ 为0的没有意义的项。
把d和i拆分成素数乘积的形式：

我们已经知道i是d的因子
那肯定就有对于i的每一个质因子 pi,Ai≤ai
下面我们来代入线筛试试看怎么求后面的和式
在这里我进行了一些讨论
假设我们有两个数 i∗prime[j],i
①如果 imodprime[j]=0
这时候i和prime[j]就满足上面d和i的关系了。
这时候先来看一种极为特殊的情况： i=prime[j] ，显然这个时候和式的值为1
再说另一种特殊情况：如果d中每个质因子 pi 的指数 ai 都相等呢？这时候每个质因子都有可能对d的函数值 g(d) 有贡献，而满足这种情况也只可能是 i=pa11∗pa22∗...∗pai−1i∗...∗pann 时(ai-1≠0)。那么这时候 ans[d]=−ans[i]
为什么？
再加入这个新的p_i之后我们看一下那个和式的内容，增加了很多项对吧。这些项有一个共同点，数值的绝对值都一样，只是正负号不一样，而符号的决定就与质因子数目的奇偶有关了。这个稍微脑补一下就能发现是对的。。。（我怎么发现的？打表大法好！）
最后就是普遍情况了，完全是一堆没有规律的a摆在那里，满足里面必定存在ai≠aj，和式的值是多少？

什么你觉得我一直打表找规律是在扯淡？
那就证！
我们知道，所有a里有一部分是对 g(d) 有贡献的，另一部分是没有贡献的。
同样数值的正负号与质因子个数有关。
我们把d拆分成i和 di （不含平方因子，否则 μ 为0就没有意义了)两个因子，偶然发现让 di 有偶个质因子和奇个平方因子的方案数恰好是相同的，在这个条件下，每种方案总有数值为其相反数的方案与之对应。所以最后和式的值就成0了。
我怎么可能在做题时候自己证着玩呢233反正打表找到规律了先写出来代码A了再说，这些事都是做完题之后发生的。
第一种情况终于搞完了= =
② imodprime[j]≠0
跟上面差不多了，不同的是如果i为互异素数组成， ans[i∗prime[j]]=−ans[i] ，否则为0.
终于全部讨论完了= =
不知道我扯淡了这么多看这篇题解的人看懂了没有（说不定早就看烦了QAQ）如果还是不懂的话，请您翻别人的题解吧浪费您这么多时间真是抱歉QAQ网上好的题解一大堆呢。。。
直接愉悦的贴代码吧。。。
(这是一份曾经在BZOJ CE了9次还是10次最后AC的代码QUQ)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdlib>
using namespace std;
#define MAXN 10000010
int num[MAXN];
bool not_prime[MAXN];
int prime[MAXN],top;
int mu[MAXN],Last[MAXN];
long long sum[MAXN];
long long ans;
int T;
int a,b;
void in(int &x)
{
    char ch=' ';x=0;
    while (!(ch>='0'&&ch<='9')) ch=getchar();
    while (ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
}
void check_prime()
{
    for (int i=2;i<=MAXN-10;i++)
    {
        if (!not_prime[i])
            prime[++top]=i,mu[i]=-1,num[i]=1,sum[i]=1,Last[i]=i;
        for (int j=1;j<=top&&i*prime[j]<=MAXN-10;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                num[i*prime[j]]=num[i]+1;
                Last[i*prime[j]]=Last[i]*prime[j];
                if (i/Last[i]==1) sum[i*prime[j]]=1;
                else 
                if (num[i/Last[i]]==num[i*prime[j]]) sum[i*prime[j]]=-sum[i/Last[i]];
                else sum[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
            num[i*prime[j]]=1;
            Last[i*prime[j]]=prime[j];
            if (num[i]==1) sum[i*prime[j]]=-sum[i];
            else sum[i*prime[j]]=0;
        }
    }
}
long long get_num(int a,int b)
{
    long long last=0;
    long long ret=0;
    for (long long i=1;i<=min(a,b);i=last+1)
    {
        last=min(a/(a/i),b/(b/i));
        ret+=(long long)(a/i)*(b/i)*(sum[last]-sum[i-1]);
    }
    return ret;
}
int main()
{
    check_prime();
    num[0]=0;num[1]=0;num[2]=1;mu[0]=0;mu[1]=1;
    for (int i=1;i<=MAXN-10;i++) sum[i]+=sum[i-1];
    in(T);
    while (T--)
    {
        in(a);in(b);
        printf("%I64d\n",get_num(a,b));
    }
}