[BZOJ]2154: Crash的数字表格&&2693: jzptab

题目描述

Description



Input

一个正整数T表示数据组数

接下来T行 每行两个正整数 表示N、M

Output

T行 每行一个整数 表示第i组数据的结果

Sample Input

1

4 5

Sample Output

122

HINT

T <= 10000

N, M<=10000000

解题思路

首先我们先直观地写出题目要我们求什么

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

这个东西就等于

∑i=1n∑j=1ni∗j/gcd(i,j)

令

f(n)=1n,f(n)=∑d|nF(d),F(n)=∑d|nμ(d)f(nd)

原式可以转为

∑i=1n∑j=1ni∗j∗f(gcd(i,j))

=∑i=1n∑j=1ni∗j∑d|i,jF(d)

=∑d=1min(n,m)F(d)∑i=1n∑j=1mi∗j [d|i & d|j]

=∑d=1min(n,m)F(d)∗d∗d∑i=1⌊nd⌋∑j=1⌊md⌋i∗j

令

S=∑i=1⌊nd⌋∑j=1⌊md⌋i∗j

S=(1+⌊n⌋d)⌊n⌋d2∗(1+⌊m⌋d)⌊m⌋d2

之后我们用f函数来代替F函数

∑d=1min(n,m)F(d)∗d∗d∗S

=∑d=1min(n,m)∑k|dμ(k)f(dk)∗d∗d∗S

=∑d=1min(n,m)∑k|dμ(k)kd∗d∗d∗S

=∑d=1min(n,m)∑k|dμ(k)∗k∗d∗S

令

h(d)=∑k|dμ(k)∗k∗d

当d为素数时，h(d)即为d-d^2，h函数为积性函数，可以直接筛出来

所以答案就是

∑d=1min(n,m)h(d)∗S

代码

#include <cstdio>
using namespace std;
#define ll long long
const int mod=100000009;
inline char tc(void)
{
static char fl[1000000],*A=fl,*B=fl;
return A==B&&(B=(A=fl)+fread(fl,1,1000000,stdin),A==B)?EOF:*A++;
}
inline int read(void)
{
int a=0;static char c;
while((c=tc())<'0'||c>'9');
while(c>='0'&&c<='9')
a=a*10+c-'0',c=tc();
return a;
}
inline void swap(int &a,int &b)
{
int t=a;a=b,b=t;
return ;
}
inline int min(int a,int b)
{
return a<b?a:b;
}
const int MAXN=10000001;
int T,tot,prim[3000001],h[MAXN+1];
char p[MAXN+1];
inline void prepare()
{
h[1]=1;
for (int i=2;i<=MAXN;++i)
{
if(!p[i])
prim[++tot]=i,h[i]=(i-((ll)i*i)%mod)%mod;
for (int j=1;j<=tot&&i*prim[j]<=MAXN;++j)
{
p[i*prim[j]]=1;
if(i%prim[j]==0)
{
h[i*prim[j]]=(ll)h[i]*prim[j]%mod;
break;
}
else h[i*prim[j]]=(ll)h[i]*h[prim[j]]%mod;
}
}
for (int i=1;i<=MAXN;++i)
h[i]=(h[i]+h[i-1])%mod;
return ;
}
inline int t(int x,int y)
{
return (((ll)x*(x+1)>>1)%mod)*(((ll)y*(y+1)>>1)%mod)%mod;
}
inline int cal(int x,int y)
{
ll sum=0;
for (int i=1,pos;i<=x;i=pos+1)
pos=min(x/(x/i),y/(y/i)),sum+=(((ll)h[pos]-h[i-1])%mod)*t(x/i,y/i)%mod,sum%=mod;
return (sum+mod)%mod;
}
int main(void)
{
register int i,x,y;
T=read(),prepare();
while(T--)
{
x=read(),y=read();
if(x>y)
swap(x,y);
printf("%d\n",cal(x,y));
}
return 0;
}