有两个注意的地方：
1.预处理到$max(n,m)$的范围就行了，要不在大视野上超时
2. (x+1)*x/2*(y+1)*y/2 在中间会炸long long，记得在中间mod
自己的n<=1000000做法：
主要的限制是我的因为要除一个东西，没法线性筛，只能log筛，n>2000000就会超时
$$ ans=sum_{i=1}^nsum_{j=1}^m lcm(i,j) $$
$$ ans=sum_{i=1}^nsum_{j=1}^mfrac{i*j}{gcd(i,j)} $$
这是题意
我的是 设f(i) 为 gcd(x,y)==i的x*y的和 F(i) 为i|gcd(x,y) 的x*y的和
反演后
$$ f(i)=sum_{i|d}mu(frac{d}{i})F(d)，其中另x=lfloorfrac{n}{i} floor，y=lfloorfrac{m}{i} floor 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; inline int read() { 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}^nsum_{j=1}^mfrac{i*j}{gcd(i,j)} $$
$$ ans=sum_{d=1}^{min(n,m)}frac{d^2*f(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor)}{d} $$
f(x,y) 为gcd(x,y)==1的x*y的和，反演后
$$ f(x,y)=sum_{d=1}^{min(n,m)}d^2mu(d)Sum(lfloorfrac{x}{d} floor,lfloorfrac{y}{d} floor)$$
$$ Sum(x,y)=(x+1)*x/2*(y+1)*y/2 $$
然后就是$O(sqrt{n}sqrt{n})$
用到的思想就是把求$gcd(x,y)==d$转化成了求
#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; inline int read() { 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(); }