BZOJ 2154 Crash的数字表格 && BZOJ 2693 jzptab

发表于 | 更新于

Link And This

求

∑i<=N∑j<=Mlcm(i,j)\sum_{i<=N}\sum_{j<=M}lcm(i,j)

Solution

Ans=∑i<=N∑j<=Mlcm(i,j)=∑i<=N∑j<=Mijgcd(i,j)Ans=\sum_{i<=N}\sum_{j<=M}lcm(i,j)=\sum_{i<=N}\sum_{j<=M}\dfrac {ij}{\gcd(i,j)}

改成枚举d=gcd(i,j)d=\gcd(i,j)

Ans=∑d∑i<=N∑j<=Mijd[gcd(i,j)==d]=∑d∑i<=N∑j<=Mijd[gcd([id],[jd])==1]Ans=\sum_d\sum_{i<=N}\sum_{j<=M}\dfrac {ij}d[\gcd(i,j)==d]=\sum_d\sum_{i<=N}\sum_{j<=M}\dfrac {ij}d[\gcd([\dfrac id],[\dfrac jd])==1]

令x=[id],y=[jd]x=[\dfrac id],y=[\dfrac jd]，则

Ans=∑d∑x<=[Nd]∑j<=[Md]xyd[gcd(x,y)==1]=∑dd{∑x<=[Nd]∑j<=[Md]xy[gcd(x,y)==1]}Ans=\sum_d\sum_{x<=[\frac Nd]}\sum_{j<=[\frac Md]}xyd[\gcd(x,y)==1]=\sum_dd\{\sum_{x<=[\frac Nd]}\sum_{j<=[\frac Md]}xy[\gcd(x,y)==1]\}

好了，现在的问题是：求

∑i<=a∑j<=bij[gcd(i,j)==1]\sum_{i<=a}\sum_{j<=b}ij[\gcd(i,j)==1]

定义F(n)=∑i<=a∑j<=bij[n∣gcd(i,j)],f(n)=∑i<=a∑j<=bij[gcd(i,j)==n]F(n)=\sum_{i<=a}\sum_{j<=b}ij[n|\gcd(i,j)],f(n)=\sum_{i<=a}\sum_{j<=b}ij[\gcd(i,j)==n]，则

F(n)=∑n∣df(d)F(n)=\sum_{n|d}f(d)f(n)=∑n∣dμ(dn)F(d)f(n)=\sum_{n|d}\mu(\dfrac dn)F(d)

要求的是

f(1)=∑dμ(d)F(d)f(1)=\sum_{d}\mu(d)F(d)

考虑求F(n)F(n)，

TeX parse error: Undefined control sequence \(

好了，

TeX parse error: Undefined control sequence \(

好了,

Ans=∑pp{∑x<=[Np]∑j<=[Mp]xy[gcd(x,y)==1]}Ans=\sum_pp\{\sum_{x<=[\frac Np]}\sum_{j<=[\frac Mp]}xy[\gcd(x,y)==1]\}

TeX parse error: Undefined control sequence \(

令T=pdT=pd，则

TeX parse error: Undefined control sequence \(

为什么我想推O(n)O(n)的却直接推出了O(n)O(\sqrt n)的 不过好歹算是自己推出来了，，虽然用了半个小时。。主要是变量重名把自己绕晕了，就重新推。。推反演式子的时候多用几个字母。。重名晕了就GG了。。 前面分块，后面前缀和。。 而后面的前缀和是一个积性函数。

积性函数

一个数论函数ff，如果在gcd(x,y)==1\gcd(x,y)==1的时候有f(xy)==f(x)f(y)f(xy)==f(x)f(y)，则称为积性函数。

  1. 积性函数的积是积性函数，这个很显然。
  2. 积性函数的约数和是积性函数 对于积性函数ff，定义g(n)=∑d∣nf(d)g(n)=\sum_{d|n} f(d) 那么g(x)g(y)=∑d1∣xf(d1)∑d2∣yf(d2)g(x)g(y)=\sum_{d_1|x}f(d_1)\sum_{d_2|y}f(d_2) 因为gcd(x,y)==1gcd(x,y)==1，所以∑d1∣xd1∑d2∣xd2\sum_{d_1|x}d_1\sum_{d_2|x}d_2依次相乘展开就是xyxy的所有约数，所以∑d1∣xf(d1)∑d2∣xf(d2)\sum_{d_1|x}f(d_1)\sum_{d_2|x}f(d_2)依次相乘展开的和就是∑d∣xyf(d)=g(xy)\sum_{d|xy}f(d)=g(xy)

有了第二个结论之后，h(T)=T∑d∣Tμ(d)dh(T)=T\sum_{d|T}\mu(d)d就是个积性函数了。因为μ\mu是个积性函数，而idid也是个积性函数(id(n)=nid(n)=n)。 积性函数就可以用线性筛筛了 

1
2
3
4
5
6
7
8
...
if(!notp[i])  prime[++pcnt]=i;
for(int j=1;j<=pcnt && prime[j]*i<=N;j++)
{
    notp[prime[j]*i]=1;
    if(i%prime[j]==0) break;
}
...

如果 ，说明h(i×primej)==h(i)h(primej)h(i\times prime_j)==h(i)h(prime_j) 否则，对于这个函数，i×primeji\times prime_j比ii多出来的因数一定都是primej2×xprime_j^2\times x的形式，它们的μ\mu都是00，所以不会对答案产生贡献。所以对答案的影响只是前面的系数ii变成了i×primeji\times prime_j，即h(i×primej)=i×primej∑d∣iμ(d)dh(i\times prime_j)=i\times prime_j\sum_{d|i}\mu(d)d

Code

BZOJ 2154

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
//Code by Lucida
#include<bits/stdc++.h>
#define red(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int P=20101009,N=1e7,MAXN=N+10;
using std::swap;
using std::min;

typedef long long LL;
int prime[MAXN],pcnt,mu[MAXN];
LL Smud2[MAXN];
void Shake(int N)
{
	static bool notp[MAXN];
	mu[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!notp[i])
		{
			prime[++pcnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=pcnt && (LL)prime[j]*i<=N;j++)
		{
			notp[prime[j]*i]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else
				mu[i*prime[j]]=-mu[i];
		}
	}
	//Smud2[1]=mu[1]*1*1;
	Smud2[0]=0;
	for(int i=1;i<=N;i++)
		Smud2[i]=(Smud2[i-1]+(LL)mu[i]*i*i%P)%P;
}
inline LL S(LL x,LL y){return ( (x*(x+1)/2%P) * (y*(y+1)/2%P) )%P;}
inline LL SMUD2(int l,int r)
{
	LL res=Smud2[r]-Smud2[l-1];
	((res%=P)+=P)%=P;
	return res;
}
LL F(int n,int m)
{
	if(n>m) swap(n,m);
	LL res=0;
	for(int i=1,last;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		(res+=(SMUD2(i,last)*S(n/i,m/i))%P)%=P;
	}
	return res;
}
inline LL Sum(LL l,LL r){return ((l+r)*(r-l+1)/2)%P;}
int main()
{
	int n,m;red(n),red(m);
	if(n>m) swap(n,m);
	Shake(n);
	LL ANS=0;
	for(int i=1,last;i<=n;i=last+1)//last=i+1..
	{
		last=min(n/(n/i),m/(m/i));
		(ANS+=(Sum(i,last)*F(n/i,m/i))%P)%=P;
	}
	((ANS%=P)+=P)%=P;
	printf("%lld\n",ANS);
	return 0;
}

BZOJ 2693

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
//Code by Lucida
#include<bits/stdc++.h>
#define red(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int N=10000000,MAXN=N+10,P=1e8+9,MAXT=10000+10;
using std::min;
using std::swap;
using std::max;
typedef long long LL;
int prime[MAXN],pcnt,mu[MAXN];
LL H[MAXN],SH[MAXN];
void Shake(int N)
{
	static bool notp[MAXN];
	mu[1]=H[1]=1;
	for(int i=2;i<=N;i++)
	{
		if(!notp[i])
		{
			prime[++pcnt]=i;
			mu[i]=-1;
			H[i]=(((LL)i*(1*mu[1]+i*mu[i])%P)+P)%P;
		}
		for(int j=1;j<=pcnt && (LL)prime[j]*i<=N;j++)
		{
			notp[prime[j]*i]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				H[i*prime[j]]=H[i]*prime[j]%P;
				break;
			}
			else
			{
				mu[i*prime[j]]=-mu[i];
				H[i*prime[j]]=H[i]*H[prime[j]]%P;
			}
		}
	}
	/*for(int i=1;i<=N;i++)
		for(int j=i;j<=N;j+=i)
			H[j]+=((LL)j*mu[i]*i)%P,(H[j]+=P)%=P;
	*/
	SH[0]=0;
	for(int i=1;i<=N;i++)
		SH[i]=(SH[i-1]+H[i])%P;
}
LL S(LL x,LL y)
{
	return ( ((x*(x+1)/2)%P)*((y*(y+1)/2)%P) )%P;
}
LL HSum(LL l,LL r)
{
	LL res=SH[r]-SH[l-1];
	((res%=P)+=P)%=P;
	return res;
}
LL Calc(LL n,LL m)
{
	LL res=0;
	if(n>m) swap(n,m);
	for(int i=1,last;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		(res+=HSum(i,last)*S(n/i,m/i))%=P;
	}
	return res;
}
int main()
{
	//freopen("input","r",stdin);
	static int Q[MAXT][2];
	int T;red(T);
	int n=0;
	for(int i=1;i<=T;i++)
	{
		red(Q[i][0]),red(Q[i][1]);
		chkmx(n,max(Q[i][0],Q[i][1]));
	}
	Shake(n);
	for(int i=1;i<=T;i++)
		printf("%lld\n",Calc(Q[i][0],Q[i][1]));
	return 0;
}