小练习

据说洛谷有此题，但我没找到。若找到了，请@我，谢谢。

Decribe:

求：

$$\sum _{i=1}^n\sum _{j=1}^{m}i\times j\times gcd(i,j)$$

Solution:

令 $k=gcd(i,j)$。

$$\sum _{k=1}^{min(n,m)}{k}^3\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }j\sum _{d|gcd(i,j)}\mu (d)$$

$$\sum _{k=1}^{min(n,m)}{k}^3\sum _{d=1}^{min(\lfloor \frac{n}{k}\rfloor ,\lfloor \frac{m}{k}\rfloor )}{d}^2\mu (d)\sum _{i=1}^{\lfloor \frac{n}{kd}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{kd}\rfloor }j$$

令 $kd=T$。

$$\sum _{T=1}^{min(n,m)}{T}^2\sum _{d|T}d\mu (\frac{T}{d})\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{T}\rfloor }j$$

$$\sum _{T=1}^{min(n,m)}{T}^2\varphi (T)\sum _{i=1}^{\lfloor \frac{n}{T}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{T}\rfloor }j$$

Code:

bool _Start;
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
namespace IO
{
	#define TP template<typename T>
	#define TP_ template<typename T,typename ... T_>
	#ifdef DEBUG
	#define gc() (getchar())
	#else
	char buf[1<<20],*p1,*p2;
	#define gc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<20,stdin),p1==p2)?EOF:*p1++)
	#endif
	#ifdef DEBUG
	void pc(const char &c)
	{
		putchar(c);
	}
	#else
	char pbuf[1<<20],*pp=pbuf;
	void pc(const char &c)
	{
		if(pp-pbuf==1<<20)
			fwrite(pbuf,1,1<<20,stdout),pp=pbuf;
		*pp++=c;
	}
	struct IO{~IO(){fwrite(pbuf,1,pp-pbuf,stdout);}}_;
	#endif
	TP void read(T &x)
	{
		x=0;static int f;f=0;static char ch;ch=gc();
		for(;ch<'0'||ch>'9';ch=gc())ch=='-'&&(f=1);
		for(;ch>='0'&&ch<='9';ch=gc())x=(x<<1)+(x<<3)+(ch^48);
		f&&(x=-x);
	}
	TP void write(T x)
	{
		if(x<0)
			pc('-'),x=-x;
		static T sta[35],top;top=0;
		do
			sta[++top]=x%10,x/=10;
		while(x);
		while(top)
			pc(sta[top--]^48);
	}
	TP_ void read(T &x,T_&...y){read(x);read(y...);}
	TP void writeln(const T x){write(x);pc('\n');}
	TP void writesp(const T x){write(x);pc(' ');}
	TP_ void writeln(const T x,const T_ ...y){writesp(x);writeln(y...);}
	TP inline T max(const T &a,const T &b){return a>b?a:b;}
	TP_ inline T max(const T &a,const T_&...b){return max(a,max(b...));}
	TP inline T min(const T &a,const T &b){return a<b?a:b;}
	TP_ inline T min(const T &a,const T_&...b){return min(a,min(b...));}
	TP inline void swap(T &a,T &b){static T t;t=a;a=b;b=t;}
	TP inline T abs(const T &a){return a>0?a:-a;}
	#undef TP
	#undef TP_
}
using namespace IO;
using std::cerr;
using LL=long long;
constexpr int N=1e7+5;
constexpr LL P=20101009;
int prime[N],pr;
LL phi[N],s[N];
bool v[N];
void init(int n=N-1)
{
	phi[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!v[i])
			prime[++pr]=i,phi[i]=i-1;
		for(int j=1;j<=pr&&prime[j]*i<=n;j++)
		{
			int k=prime[j]*i;
			v[k]=1;
			if(!(i%prime[j]))
			{
				phi[k]=phi[i]*prime[j]%P;
				break;
			}
			phi[k]=phi[i]*phi[prime[j]]%P;
		}
	}
	for(int i=1;i<=n;i++)
	{
		phi[i]=(phi[i]*i%P*i%P+phi[i-1])%P;
		s[i]=(i+s[i-1])%P;
	}
}
LL calc(LL n)
{
	return n*(n+1ll)%P*(2ll*n+1ll)%P*(16750841ll)%P;
}
bool _End;
int main()
{
//	fprintf(stderr,"%.2 MBlf\n",(&_End-&_Start)/1048576.0);
	LL n,m;
	read(n,m);
	if(n>m)
		swap(n,m);
	init(m);
	LL ans=0;
	for(LL l=1,r;l<=n;l=r+1)
	{
		r=min(n/(n/l),m/(m/l));
		ans=(ans+((phi[r]-phi[l-1])+P)%P*s[n/l]%P*s[m/l]%P)%P;
	}
	writeln(ans);
	return 0;
}