膜拜cdc……他的推导详细到我这种蒟蒻都能看得懂!
所以我附一下代码就好了。
#include<bits/stdc++.h> #define N 10000005 #define yql 20101009 using namespace std; typedef long long ll; int mu[N],prime[N],cnt,s[N],vis[N]; ll n,m,ans,maxn; inline ll sum(ll x,ll y){ return ((x*(x+1)/2)%yql)*((y*(y+1)/2)%yql)%yql; } inline void calcmu(){ memset(vis,1,sizeof(vis));cnt=0;mu[1]=1; for(int i=2;i<=maxn;i++){ if(vis[i]){prime[++cnt]=i;mu[i]=-1;} for(int j=1;j<=cnt;j++){ int t=i*prime[j];if(t>maxn)break; vis[t]=0; if(i%prime[j]==0){mu[t]=0;break;} mu[t]=-mu[i]; } } for(ll i=1;i<=maxn;i++)s[i]=(s[i-1]+(i*i*mu[i])%yql)%yql; } inline ll F(ll x,ll y){ ll ans=0;ll n=min(x,y); for(ll i=1,j=1;i<=n;i=j+1){ j=min(x/(x/i),y/(y/i)); ans=(ans+(s[j]-s[i-1])*sum(x/i,y/i)%yql)%yql; } return ans; } int main(){ cin>>n>>m;maxn=min(n,m); calcmu(); for(ll i=1,j=1;i<=maxn;i=j+1){ j=min(n/(n/i),m/(m/i)); ans=(ans+(i+j)*(j-i+1)/2%yql*F(n/i,m/i)%yql)%yql; } printf("%lld ",(ans+yql)%yql); }