BZOJ2154: Crash的数字表格

2154: Crash的数字表格

Time Limit: 20 Sec  Memory Limit: 259 MB
Submit: 4300  Solved: 1563
[Submit][Status][Discuss]

Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Sample Input

4 5

Sample Output

122
【数据规模和约定】
100%的数据满足N, M ≤ 10^7。

思路{

　　lcm(i,j)=i*j/gcd(i,j)

　　Ans=∑(1<=i<=n,1<=j<=m)i*j/gcd(i,j) 枚举GCD

　　Ans=(min(n,m))∑(d=1){d^2*(F((n/d),(m/d),1))/d}

　　=(min(n,m))∑(d=1){d*(F((n/d),(m/d),1))}

　　定义F(x,y,d)=∑(1<=i<=x,1<=j<=y,gcd(i,j)==d)i*j

　　G(x,y,d)=∑(1<=i<=x,1<=j<=y,gcd(i,j)|d)i*j

　　S(x,y)=∑(1<=i<=x,1<=j<=y)i*j

　　S(x,y)=(x)∑(i=1)i*(y)∑(j=1)j=((x+1)*x/2)*((y+1)*y/2)

　　那么G(x,y,d)=S((x/d),(y/d))*d^2;

　　反演理论得:F(x,y,1)=(min(n,m))∑(i=1)μ(i)*S((x/i),(y/i))*i^2;

　　那么数论分块套数论分块.还有要使用逆元.

}

#include<bits/stdc++.h>
#define RG register
#define il inline
#define db double
#define LL long long
#define N 10000010
#define mod 20101009
using namespace std;
int mu[N],p[N];LL G[N];bool vis[N];
LL S(LL x){return (((x%mod)*((x+1)%mod)%mod)*10050505)%mod;}
void pre(){
  mu[1]=1;
  for(int i=2;i<N;++i){
    if(!vis[i])p[++p[0]]=i,mu[i]=-1;
    for(int j=1;j<=p[0]&&p[j]*i<N;++j){
      vis[p[j]*i]=true;
      if(i%p[j])mu[i*p[j]]=-mu[i];
      else {
    mu[i*p[j]]=0;
    break;
      }
    }
  }
  for(LL i=1;i<N;++i)G[i]=G[i-1]+mu[i]*i*i,G[i]=(G[i]+mod)%mod;
}
LL F(int n,int m){
  if(n>m)swap(n,m);LL Ans(0);
  for(LL l=1,r;l<=n;l=r+1){
    r=min(n/(n/l),m/(m/l));
    Ans+=(((G[r]-G[l-1]+mod)%mod)*(S(n/l)*S(m/l)%mod))%mod;
    Ans=(Ans+mod)%mod;
  }return Ans;
}
int main(){
  int n,m;
  scanf("%d%d",&n,&m);pre();
  if(n>m)swap(m,n);LL ans(0);
  for(LL l=1,r;l<=n;l=r+1){
    r=min(n/(n/l),m/(m/l));
    LL tmp=(r-l+1)*(l+r)/2;
    tmp%=mod;
    ans+=(tmp*F((n/l),(m/l)))%mod;
    if(ans>=mod)ans-=mod;
  }cout<<ans;
  return 0;
}