# bzoj_2154: Crash的数字表格 (自己想了70分)

1.预处理到$max(n,m)$的范围就行了，要不在大视野上超时

2. (x+1)*x/2*(y+1)*y/2 在中间会炸long long，记得在中间mod

$$ans=\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)$$

$$ans=\sum_{i=1}^n\sum_{j=1}^m\frac{i*j}{gcd(i,j)}$$

$$f(i)=\sum_{i|d}\mu(\frac{d}{i})F(d)，其中另x=\lfloor\frac{n}{i}\rfloor，y=\lfloor\frac{m}{i}\rfloor F(i)=i^2(x+1)*x/2*(y+1)*y$$

$$ans=\sum_{d=1}^{min(n,m)}x*(x+1)*y*(y+1)\sum_{i|d}\frac{\mu(\frac{d}{i})}{4*i}$$

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#define ll long long
#define mem(a,b) memset(a,b,sizeof(a))
using namespace std;
{
char q=getchar();int ans=0;
while(q<'0'||q>'9')q=getchar();
while(q>='0'&&q<='9'){ans=ans*10+q-'0';q=getchar();}
return ans;
}
const int mod=20101009;
const int N=2000006;

int prime[N],cnt;
bool he[N];
int mu[N];
ll ni[N],ji[N];

void chu()
{
ni[1]=1;
for(int i=2;i<N;++i)
ni[i]=(ll)(mod-mod/i)*ni[mod%i]%mod;

mu[1]=1;//ji[1]=1;
for(int i=2;i<N;++i)
{
if(!he[i])
{
prime[++cnt]=i;
mu[i]=-1;
//ji[i]=(ni[i]-1+mod)%mod;
}
for(int j=1;j<=cnt&&prime[j]*i<N;++j)
{
he[i*prime[j]]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
//ji[i*prime[j]]=ji[i]*ni[prime[j]]%mod;
break;
}
//ji[i*prime[j]]=ji[i]*ji[prime[j]]%mod;
mu[i*prime[j]]=-mu[i];
}
}

for(int i=1;i<N;++i)
for(int j=i;j<N;j+=i)
ji[j]=(ji[j]+mu[j/i]*ni[i]%mod*ni[4]%mod*j%mod*j%mod+mod)%mod;
for(int i=1;i<N;++i)
ji[i]=((ji[i]+ji[i-1])%mod+mod)%mod;
}

int n,m;

ll work()
{
if(n>m)
swap(n,m);
ll ans=0,tt;
int nx;
for(int i=1;i<=n;i=nx+1)
{
nx=min( n/(n/i),m/(m/i) );
tt=(ll)(n/i+1)%mod*(n/i)%mod*(m/i)%mod*(m/i+1)%mod;
ans=(ans+tt*(ji[nx]-ji[i-1]+mod)%mod)%mod;
}
return (ans+mod)%mod;
}

int main(){

freopen("in.in","r",stdin);
//freopen("nt2011_table.in","r",stdin);
//freopen("nt2011_table.out","w",stdout);

chu();
scanf("%d%d",&n,&m);
cout<<work();
}


$$ans=\sum_{i=1}^n\sum_{j=1}^m\frac{i*j}{gcd(i,j)}$$

$$ans=\sum_{d=1}^{min(n,m)}\frac{d^2*f(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}{d}$$

f(x,y) 为gcd(x,y)==1的x*y的和，反演后

$$f(x,y)=\sum_{d=1}^{min(n,m)}d^2\mu(d)Sum(\lfloor\frac{x}{d}\rfloor,\lfloor\frac{y}{d}\rfloor)$$

$$Sum(x,y)=(x+1)*x/2*(y+1)*y/2$$

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#define ll long long
#define mem(a,b) memset(a,b,sizeof(a))
using namespace std;
{
char q=getchar();int ans=0;
while(q<'0'||q>'9')q=getchar();
while(q>='0'&&q<='9'){ans=ans*10+q-'0';q=getchar();}
return ans;
}
const int mod=20101009;
const int N=10000006;

int prime[N],cnt;
bool he[N];
int mu[N];
ll pr1[N],pr2[N];
int n,m;

void chu()
{
int q1=max(n,m);

mu[1]=1;
for(int i=2;i<=q1;++i)
{
if(!he[i])
{
prime[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt&&prime[j]*i<=q1;++j)
{
he[i*prime[j]]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}

for(int i=1;i<=q1;++i)
{
pr1[i]=(pr1[i-1]+i)%mod;
pr2[i]=(pr2[i-1]+((ll)mu[i]*i%mod*i%mod)+mod)%mod;
}
}

ll sum(ll x,ll y)
{
return ( ((x+1)*x/2%mod) * ((y+1)*y/2%mod) )%mod;
}

ll f(int x,int y)
{
if(x>y)
swap(x,y);
ll ans=0;
int nx;
for(int i=1;i<=x;i=nx+1)
{
nx=min( x/(x/i),y/(y/i) );
ans=(ans+sum(x/i,y/i)*(pr2[nx]-pr2[i-1]+mod)%mod)%mod;
//nx=min( x/(x/i),y/(y/i) );
//ans=(ans+((x/i+1)*(x/i)/2)*((y/i+1)*(y/i)/2)%mod*(pr2[nx]-pr2[i-1]+mod)%mod)%mod;
}
return (ans+mod)%mod;
}

ll work()
{
if(n>m)
swap(n,m);
ll ans=0;
int nx;
for(int i=1;i<=n;i=nx+1)
{
nx=min( n/(n/i),m/(m/i) );
ans=(ans+(pr1[nx]-pr1[i-1]+mod)%mod*f(n/i,m/i))%mod;
}
return (ans+mod)%mod;
}

int main(){

//freopen("in.in","r",stdin);
//freopen("nt2011_table.in","r",stdin);
//freopen("nt2011_table.out","w",stdout);

scanf("%d%d",&n,&m);
chu();
cout<<work();
}


posted @ 2017-10-16 16:42  A_LEAF  阅读(114)  评论(0编辑  收藏  举报