bzoj 2693 jzptab（莫比乌斯反演+积性函数线性筛一般套路）

题目：
Description
求 ${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$
Input
一个正整数T表示数据组数，接下来T行每行两个正整数表示N和M
(T<=10000,N,M<=10000000)
Output
T行，每行一个整数表示第i组数据的结果
Sample Input
1
4 5
Sample Output
122

参考博客：积性函数
参考1
参考2

这题先参考上一篇博文得到的最终式子：last

我们会得到这个式子：
${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)=$${\sum }_{d=1}^{n}d{\sum }_{i=1}^{n/d}u(i)\ast i^{2}sum(\frac{n}{d\ast i},\frac{m}{d\ast i})$

然后我们可以化简为：$\frac{令T=d\ast i}{=}{\sum }_{T=1}^{n}sum(\frac{n}{T},\frac{m}{T}){\sum }_{i|T}{i}^{2}u(i)\frac{T}{i}$
这个化简说实话，我是想不到。

我们设 $h(T)$$={\sum }_{i|T}{i}^{2}u(i)\frac{T}{i}$，可以知道函数h为积性函数，因为我们知道正因子之和是积性函数，莫比乌斯函数也是积性函数，两个积性函数相乘也为积性函数，如何判断此函数是否为积性函数，像我这种萌新，更多的只能依靠平时的积累加经验（大部分题目给的，你都能看成积性函数去做，只要你能想到推导公式）。

既然现在我们知道函数h为积性函数，那么我们就可以线性筛把函数h求出来，再弄个前缀和就行了。

线性筛分为3部分：
1.n本身是素数，这个根据积性函数的定义可得，很容易求。
2.i%prime[j]!=0，这个也是根据积性函数的性质可得，即f(a)f(b)=f(ab)。
3.i%prime[j]==0，这个可能需要找规律加灵光一闪。

1，对于素数p，它的函数$h(p)=p-p\ast p$;
2，i%prime[j]!=0，此时$h[i\ast prime[j]]=h[i]\ast$ $h[prime[j]]$


3，i%prime[j]==0，此时$h[i\ast prime[j]]=h[i]\ast$ $prime[j]$，i*prime[j] 的唯一分解可以写成一部分素数乘积 i 与 另一部分在前一部分中出现过的素数的乘积 j，也就是存在质因子的指数大于1，它新增的每一个因子的 μ 值都是0，没有意义，只有统计时D变成了原来的 j 倍

所以 此时 $h(i\ast prime[j])=h(i)\ast prime[j]$。
最后我们处理下函数h的前缀和，再加个分块就可以了。

#include<bits/stdc++.h>

using namespace std;
typedef long long LL;

const LL mod=1e8+9;
const int N=1e7+10;

int prime[N],tot;
bool vis[N];

LL f[N],sum[N];
int n,m;

void init()
{
    f[1]=sum[1]=1LL;
    for(int i=2;i<N;i++)
    {
        if(!vis[i]){
            prime[++tot]=i;
            f[i]=(i-1LL*i*i+mod)%mod;///i为素数
        }
        for(int j=1;j<=tot&&i*prime[j]<N;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){///第三种情况
                f[i*prime[j]]=f[i]*prime[j]%mod;
                break;
            }
            ///第二种情况
            f[i*prime[j]]=f[i]*f[prime[j]]%mod;
        }
        sum[i]=(sum[i-1]+f[i])%mod;

    }
}

LL Sum(LL x,LL y){

    x=((x+1)*x/2)%mod;
    y=((y+1)*y/2)%mod;
    return x*y%mod;
}

LL solve()
{
    int r;
    LL re=0;
    ///整数分块
    for(int l=1;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        re=(re+Sum(n/l,m/l)*(sum[r]-sum[l-1])%mod)%mod;
    }
    return re;
}

int main()
{
    init();

    int ncase;
    scanf("%d",&ncase);

    while(ncase--)
    {
        scanf("%d%d",&n,&m);
        if(n>m) swap(n,m);

        printf("%lld\n",solve());
    }
    return 0;
}

我的标签：做个有情怀的程序员。