BZOJ3309 DZY Loves Math

原题链接：https://www.lydsy.com/JudgeOnline/problem.php?id=3309

DZY Loves Math

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
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957

Sample Output

35793453939901
14225956593420
4332838845846
15400094813

HINT

【数据规模】

T<=10000
1<=a,b<=10^7

题解

$\mathcal{BZOJ}$的题面是真的丑，看公式容易瞎。。。

式子的化简还是挺简单的：

$$\begin{align*} &\sum_{i=1}^a\sum_{j=1}^bf(gcd(i,j))\\ &=\sum_{d=1}^{min(a,b)}f(d)\sum_{i=1}^a\sum_{j=1}^b[d=gcd(i,j)]\\ &=\sum_{d=1}^{min(a,b)}f(d)\sum_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{b}{d}\rfloor}\sum_{d'|gcd(i,j)}\mu(d')\\ &=\sum_{d=1}^{min(a,b)}f(d)\sum_{d'=1}^{\lfloor\frac{min(a,b)}{d}\rfloor}\mu(d')\sum_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{b}{d}\rfloor}[d'|gcd(i,j)]\\ &=\sum_{d=1}^{min(a,b)}f(d)\sum_{d'=1}^{\lfloor\frac{min(a,b)}{d}\rfloor}\mu(d')\lfloor\frac{a}{dd'}\rfloor\lfloor\frac{b}{dd'}\rfloor\\ &=\sum_{T=1}^{min(a,b}\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor\sum_{d|T}f(d)\mu(\frac{T}{d}) \end{align*}$$

设$F(T)=\sum_{d|T}f(d)\mu(\frac{T}{d})$，显然这玩意儿不是积性的（废话$F(1)=0$你让我怎么积性）。

但是我们已经没法化简了啊，没办法，硬着头皮讨论一波：

设$T=\prod_{i=1}^kp_i^{a_i}$，对于$\frac{T}{d}$，$p_i$的指数只能是$0$或$1$，否则$\mu(\frac{T}{d})$的值为$0$，没有意义。这样，对于一个确定的$f(d)$，指数和为奇或为偶的方案数相同，直接抵消成$0$。

但是，对于$T=\prod_{i=1}^kp_i^a$，即质数的指数都一样的时候，当$\frac{T}{d}=\prod_{i=1}^kp_i$时，$f(d)=a-1$，而其他方案的贡献都是$a$，无法抵消，所以需要额外减掉$1$，再考虑上$\mu(\frac{T}{d})=(-1)^k$的贡献，此时$F(T)=(-1)^{k+1}$。

综上，虽然$F(T)$不是积性的，但是我们仍然可以通过讨论强行套进线筛板子里，只需要记录一下最小质因子的质数以及最小质因子的幂就可以判断一下是否存在质数指数都相等的情况。

推式子的精髓不在于化简，在于你推不下去的时候能把你搞出来辣鸡函数的筛法讨论出来。

代码

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int M=1e7+5;
ll f[M],prod[M],mi[M];
int p[M/3],T;
bool isp[M];
void pre()
{
    for(int i=2,t;i<M;++i)
    {
        if(!isp[i])p[++p[0]]=i,f[i]=mi[i]=1,prod[i]=i;
        for(int j=1;j<=p[0];++j)
        {
            t=i*p[j];if(t>=M)break;isp[t]=1;
            if(i%p[j]==0)
            {
                mi[t]=mi[i]+1,prod[t]=prod[i]*p[j];
                if(i==prod[i])f[t]=1;
                else f[t]=(mi[i/prod[i]]==mi[t]?-f[i/prod[i]]:0);
                break;
            }
            prod[t]=p[j],mi[t]=1,f[t]=(mi[i]==1?-f[i]:0);
        }
    }
    for(int i=1;i<M;++i)f[i]+=f[i-1];
}
ll calc(int a,int b)
{
    if(a>b)swap(a,b);ll ans=0;
    for(int l=1,r;l<=a;l=r+1)r=min(a/(a/l),b/(b/l)),ans+=1ll*(a/l)*(b/l)*(f[r]-f[l-1]);
    return ans;
}
void in(){scanf("%d",&T);}
void ac(){pre();int a,b;while(T--)scanf("%d%d",&a,&b),printf("%lld\n",calc(a,b));}
int main(){in();ac();}