BZOJ2154 Crash的数字表格及 BZOJ2693 JZPTAB

做这两道题做得心力交瘁，太恶心了，好久没做数学题，套路差点都全忘光了，推式子的时候犯了一堆错：

• 没有注意是否仍然满足乘法分配率
• 模数看错

Problem

BZOJ2154

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

Solution

变一下原式

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

$$\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m[\gcd (i,j)=d]\frac{ij}{d}$$

$$\sum _{d=1}^n\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}[\gcd (i,j)=1]ijd$$

$$\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}[\gcd (i,j)=1]ij$$

---

记$x=n/d,y=m/d$，那么我们来考虑下式怎么计算：

$$\sum _{i=1}^x\sum _{j=1}^y[\gcd (i,j)=1]ij$$

设$$F(d)=\sum _{i=1}^x\sum _{j=1}^y[\gcd (i,j)=d]ij$$

$$G(x)=\sum _{x|d}F(d)$$

通过莫反我们可以知道$F(n)={\sum }_{n|d}\mu (\frac{d}{n})G(d)$，为了方便表达，我们记$S(n)={\sum }_{i=1}^{n}i=\frac{n(n+1)}{2}$

$$G(d)=\sum _{i=1}^x\sum _{j=1}^y[d|\gcd (i,j)]ij=d^2\sum _{i=1}^{x/d}\sum _{j=1}^{y/d}ij=d^{2}S(\frac{x}{d})S(\frac{y}{d})$$

我们要求的就是
$$F(1)=\sum _{i=1}^x\mu (i)G(i)=\sum _{i=1}^x\mu (i)i^{2}S(\frac{x}{i})S(\frac{y}{i})$$

数论分块之，一次计算$O(\sqrt{n})$

---

$$\sum _{d=1}^{n}dF(1)=\sum _{d=1}^{n}d\sum _{i=1}^x\mu (i)G(i)$$

$$\sum _{d=1}^{n}d\sum _{i=1}^x\mu (i)i^{2}S(\frac{x}{i})S(\frac{y}{i})$$

对$x$分块，然后调用上一个讨论的计算函数计算后面式子的值

$O(\sqrt{n}\times \sqrt{n})=O(n)$

Code

#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=10000010,mod=20101009;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    if(f) x=-x;
}
int n,m,tot,ans,pri[maxn],vis[maxn],mu[maxn],sqr[maxn];
int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
void init()
{
	mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!vis[i]) mu[i]=mod-1,pri[++tot]=i;
		for(int j=1;j<=tot&&i*pri[j]<=n;j++)
		{
			vis[i*pri[j]]=1;
			mu[i*pri[j]]=dec(0,mu[i]);
			if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
		}
	}
	for(int i=1;i<=n;i++) sqr[i]=(sqr[i-1]+(ll)i*i%mod*mu[i])%mod;
}
int query(int x,int y)
{
	int res=0;
	for(int i=1,j,tmp;i<=x;i=j+1)
	{
		j=min(x/(x/i),y/(y/i));
		tmp=((ll)(x/i+1)*(x/i)/2%mod)*((ll)(y/i+1)*(y/i)/2%mod)%mod;
		res=(res+(ll)tmp*dec(sqr[j],sqr[i-1]))%mod;
	}
	return res;
}
int main()
{
	#ifndef ONLINE_JUDGE
	freopen("in.txt","r",stdin);
	#endif
	read(n);read(m);
	if(n>m) swap(n,m);
	init();
	for(int i=1,j,tmp;i<=n;i=j+1)
	{
		j=min(n/(n/i),m/(m/i));
		tmp=(ll)(i+j)*(j-i+1)/2%mod;
		ans=(ans+(ll)tmp*query(n/i,m/i))%mod;
	}
	printf("%d\n",ans);
	return 0;
}

Problem

BZOJ2693
同上题，有$10^4$组询问。

Solution

直接拿之前那题的结论接着做：
$$\sum _{d=1}^{n}d\sum _{i=1}^{n/d}\mu (i)i^{2}S(\frac{n}{id})S(\frac{m}{id})$$

枚举$id$

$$\sum _{id=1}^{n}S(\frac{n}{id})S(\frac{m}{id})\sum _{i|id}\mu (i)i^2\times \frac{id}{i}$$

$$\sum _{id=1}^{n}S(\frac{n}{id})S(\frac{m}{id})\sum _{i|id}id\mu (i)i$$

后面的是个积性函数，前面的数论分块。

注意$id$不能拆出第二个求和，否则会不满足乘法分配率导致分块做的前缀和有问题。。

时间复杂度$O(n+m\sqrt{n})$

Code

#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=10000010,mod=1e8+9;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    if(f) x=-x;
}
int z,n,m,tot,ans,pri[maxn],vis[maxn],val[maxn];
int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
void init()
{
	val[1]=1;
	for(int i=2;i<=10000000;i++)
	{
		if(!vis[i]) val[i]=dec(i,(ll)i*i%mod),pri[++tot]=i;
		for(int j=1;j<=tot&&i*pri[j]<=10000000;j++)
		{
			vis[i*pri[j]]=1;
			if(i%pri[j]==0){val[i*pri[j]]=(ll)val[i]*pri[j]%mod;break;}
			val[i*pri[j]]=(ll)val[i]*val[pri[j]]%mod;
		}
	}
	for(int i=2;i<=10000000;i++) val[i]=pls(val[i],val[i-1]);
}
int main()
{
	#ifndef ONLINE_JUDGE
	freopen("in.txt","r",stdin);
	#endif
	read(z);
	init();
	while(z--)
	{
		read(n);read(m);ans=0;
		if(n>m) swap(n,m);
		for(int i=1,j,tmp;i<=n;i=j+1)
		{
			j=min(n/(n/i),m/(m/i));
			tmp=((ll)(n/i)*(n/i+1)/2%mod)*((ll)(m/i)*(m/i+1)/2%mod)%mod;
			ans=(ans+(ll)tmp*dec(val[j],val[i-1]))%mod;
		}
		printf("%d\n",ans);
	}
	return 0;
}