bzoj 2693: jzptab （反演）

题目描述

传送门

题目大意：

∑i=1n∑j=1mlcm(i,j) mod 100000009

题解

哎，虽然式子很恶心，但是还是要化的 ,设n<=m

∑i=1n∑j=1mlcm(i,j)

∑i=1n∑j=1mi∗jgcd(i,j)

∑k=1n∑i=1n∑j=1mi∗jk[gcd(i,j)=k]

∑k=1n1k∑i=1⌊nk⌋∑j=1⌊mk⌋i∗k∗j∗k[gcd(i,j)=1]

∑k=1nk∑i=1⌊nk⌋∑j=1⌊mk⌋i∗j∑d|(i,j)μ(d)

∑k=1nk∑nd=1μ(d)∑i=1⌊nk⌋[d|i]∗i∑j=1⌊mk⌋[d|j]∗j

∑k=1nk∑nd=1μ(d)∑i=1⌊nk∗d⌋d∗i∑j=1⌊mk∗d⌋d∗j

设sum(x,y)=∑i=1xi∑j=1yj 这个式子直接用等差数列求和公式计算即可。

那么上面的式子就是

∑k=1nk∑d=1nμ(d)∗d2sum(⌊nk∗d⌋,⌊mk∗d⌋)

然后设T=k∗d ,转变枚举顺序

∑T=1nsum(⌊nT⌋,⌊mT⌋)∗∑d|nμ(Td)∗Td2∗d

那么问题的关键就转换成了如何快速的求解h(n)=∑d|nμ(Td)∗Td2∗d

这个函数是一个积性函数，可以线性筛啊。a,b互质的时候，直接相乘。

a,b不互质的时候，直接乘上质数。h(pixi)=μ(pi)∗pi2∗pixi−1+μ(1)∗pixi

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define LL long long
#define p 100000009
#define N 10000000
using namespace std;
int n,m,T,pd[N+3],prime[N+3],g[N+3];
LL f[N+3],inv;
void init()
{
f[1]=1;
for (int i=2;i<=N;i++) {
if (!pd[i]) {
prime[++prime[0]]=i;
f[i]=(LL)(i-(LL)i*i%p)%p;
}
for (int j=1;prime[j]*i<=N;j++) {
pd[i*prime[j]]=1;
if (i%prime[j]==0) {
f[i*prime[j]]=(LL)f[i]*prime[j]%p;
break;
}
f[i*prime[j]]=(LL)f[i]*f[prime[j]]%p;
}
}
for (int i=1;i<=N;i++) f[i]=(f[i-1]+f[i])%p;
}
LL calc(LL x,LL y)
{
LL t1=((x+1)*x>>1)%p; LL t2=((y+1)*y>>1)%p;
return t1*t2%p;
}
int main()
{
freopen("a.in","r",stdin);
freopen("my.out","w",stdout);
init(); scanf("%d",&T);
inv=g[4];
while (T--) {
scanf("%d%d",&n,&m);
if (n>m) swap(n,m);
int j=0; LL ans=0;
for (int i=1;i<=n;i=j+1) {
j=min(n/(n/i),m/(m/i));
ans=(ans+calc(n/i,m/i)*(f[j]-f[i-1])%p)%p;
}
printf("%lld\n",(ans%p+p)%p);
}
}