没有题目地址…被删了…
题目大意
对于 T<=10000T<=10000 组数据 n,m<=20000000n,m<=20000000
求
∑i=1n∑j=1nLCM(i,j)∑i=1n∑j=1nLCM(i,j)
开始沉浸于数学的海洋吧
先改一改题目
∑i=1n∑j=1mi∗jgcd(i,j)∑i=1n∑j=1mi∗jgcd(i,j)
这样不好处理
我们干脆枚举gcd
对于
T=gcd(i,j)T=gcd(i,j)
N=min(n,m)N=min(n,m)
f(n,m,T)=∑i=1n∑j=1m(i∗j)[gcd(i,j)=T]f(n,m,T)=∑i=1n∑j=1m(i∗j)[gcd(i,j)=T]
Ans=∑T=1Nf(n,m,T)TAns=∑T=1Nf(n,m,T)T
然后你会惊讶地发现,若设
F(n,m,T)=∑i=1n∑j=1mi∗j[T|gcd(i,j)]F(n,m,T)=∑i=1n∑j=1mi∗j[T|gcd(i,j)]
sum(n,m)=n(n+1)2∗m(m+1)2sum(n,m)=n(n+1)2∗m(m+1)2
则有
F(n,m,1)=sum(n,m)F(n,m,1)=sum(n,m)
F(n,m,T)=T2∗sum(⌊nT⌋,⌊mT⌋)F(n,m,T)=T2∗sum(⌊nT⌋,⌊mT⌋)
F(n,m,T)=∑T|tNf(n,m,t)F(n,m,T)=∑T|tNf(n,m,t)
f(n,m,T)=∑i=1NTF(nT,mT,iT)∗μ(i)f(n,m,T)=∑i=1NTF(nT,mT,iT)∗μ(i)
于是
Ans=∑T=1NT∗f(⌊nT⌋,⌊mT⌋,1)=∑T=1NT∗∑nt=1TF(⌊nT⌋,⌊mT⌋,nt)∗μ(nt)Ans=∑T=1NT∗f(⌊nT⌋,⌊mT⌋,1)=∑T=1NT∗∑nt=1TF(⌊nT⌋,⌊mT⌋,nt)∗μ(nt)
Ans=∑T=1NT∗∑nt=1Tnt2sum(⌊nntT⌋,⌊mntT⌋)∗μ(nt)Ans=∑T=1NT∗∑nt=1Tnt2sum(⌊nntT⌋,⌊mntT⌋)∗μ(nt)
强行换一下位置
由于 nt2sum(⌊nntT⌋,⌊mntT⌋)nt2sum(⌊nntT⌋,⌊mntT⌋)完全是由μ(nt)μ(nt)决定的
换句话说 ntTntT 是重复出现的 并且 nn 也都是一样的
所以 我们只要知道了 μμ 的前缀和之类的东西 就能快速求出一段的和
所以 交换位置得
Ans=∑ntT=1Nsum(⌊nntT⌋,⌊mntT⌋)∑nt|ntTμ(nt)nt2ntTntAns=∑ntT=1Nsum(⌊nntT⌋,⌊mntT⌋)∑nt|ntTμ(nt)nt2ntTnt
后面那一堆 全都是积性函数
分类讨论一下吧
对于 ∑k|Kμ(k)k2Kk∑k|Kμ(k)k2Kk
KK为质数 毫无疑问 就是
K=LK∗prime(nk)K=LK∗prime(nk) LK的值已知 nk在LK中已经出现过了 那么nk的 μ(nk)=0μ(nk)=0 贡献就为0了
但是由于 KLK=nkKLK=nk 所以K的值要乘一个nk
K=LK∗prime(nk)K=LK∗prime(nk) LK的值已知 nk与LK互质 那不就两个值乘起来
然后就分块处理一下
附上理解代码
#include <iostream>
#include <cstdio>
using namespace std;
const long long mod=1e8+9;
inline long long input()
{
char c=getchar();long long o;
while(c>57||c<48)c=getchar();
for(o=0;c>47&&c<58;c=getchar())o=(o<<1)+(o<<3)+c-48;
return o;
}
long long sum[10001234],h[10001234],prime[20000];
bool mark[10001234];
void init()
{
h[1]=sum[1]=1;
int cnt=0;long long num;
for(long long i=2;i<=10001000;i++)
{
if(!mark[i]){h[i]=(i-i*i)%mod;prime[++cnt]=i;}
for(int j=1;j<=cnt;j++)
{
num=i*prime[j];
if(num>10001000)break;
mark[num]=1;
if(i%prime[j]==0){h[num]=prime[j]*h[i]%mod;break;}
h[num]=h[prime[j]]*h[i]%mod;
}
sum[i]=(sum[i-1]+h[i])%mod;
}
}
inline long long CAL(long long a,long long b){return (a*(1+a)/2%mod)*(b*(1+b)/2%mod)%mod;}
void RES(long long n,long long m)
{
long long pos,res=0;
if(n>m)swap(n,m);
for(long long i=1;i<=n;i=pos+1)
{
pos=min(n/(n/i),m/(m/i));
res=(res+CAL(n/i,m/i)*(sum[pos]-sum[i-1]))%mod;
}
printf("%lld\n",(res+mod)%mod);
}
int main()
{
freopen("jzptab.in","r",stdin);
freopen("jzptab.out","w",stdout);
int T;long long n,m;
scanf("%d",&T);
init();
while(T--)
{
n=input();m=input();
RES(n,m);
}
return 0;
}