莫比乌斯3题小结 hdu4746,1695 && poj2154

我们约定,[x]表示x取下整

hdu4746:

   n   m

求∑  ∑  f(gcd(a,b))<=p       f为因子个数

   i=1j=1

                            min(n,m) min(n/d,m/d)

我们可以看成是:∑          ∑   u(t)* [n/dt]*[m/dt]      其中d的因子个数<=p

                             d=1       t=1

令T=dt

  min(n,m)                                                                                              min(n,m)

有 ∑       ∑  u(T/d)[n/T][m/T]     ------->[n/T][m/T]移到前面去,也就是   ∑    [n/T][m/T]  ∑ u(T/d)  &&d的因子数<=p

   T=1   d|T                                                                                              T=1                  d|T


如果能先处理出G(T)=∑ u(T/d) d的因子数<=p就可以sqrt(n)分块了。

                                   d|T

一开始我想直接每次跑一下,发现时间太长,然后YY了一个G[T][p]表示T的所有因子数<=p的d的u(T/d)的和

然后不知所措了。。。。              p一定不会超过20(logn约等于18.8),这个一定可以预处理。

然后发现大神都是先处理G[T][p]中d==p的u(T/d)之和,然后之前的定义看成当前数组的前缀和。搞下前缀和就ok

(orz前缀和)

然后对于每个p搞一下前缀,就可以分块做了。

(吐槽:离成功就差那么一点点,卧槽没想到前缀和啊)

思路在代码里面:

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<cmath>
#define ll long long
using namespace std;

/*
        min(n,m)
推出来: ∑    [n/T][m/T]∑  u(T/d)
         T=1             d|T&&f(d)<=p
*/
// f[T][p]表示T的所有d|T且f(d)<=p的  是g[T][p](d|T&&f(d)==p)的前缀和
const ll maxn=500000+5;
ll num[maxn]; 
ll prime[maxn];
bool flag[maxn];
ll tot;
ll u[maxn];
ll f[maxn][25];
void init()
{
	u[1]=1;
	num[1]=0;
	tot=0;
	for(ll i=2;i<=maxn;i++)
	{
		if(!flag[i])
		{
			prime[++tot]=i;
			u[i]=-1;
			num[i]=1;
		}
		for(ll j=1;j<=tot&&i*prime[j]<=maxn;j++)
		{
			flag[i*prime[j]]=1;
			num[i*prime[j]]=num[i]+1;
			if(i%prime[j]==0)
			{
				u[i*prime[j]]=0;
				break;
			}
			else u[i*prime[j]]=-u[i];
		}
	}
	memset(f,0,sizeof(f));
	for(ll i=1;i<=maxn;i++)
	{
		for(ll j=1;i*j<=maxn;j++)
		{
			f[i*j][num[i]]+=(ll)u[j];
		}
	}
	for(ll i=1;i<=maxn;i++)
	{
		for(ll j=1;j<=20;j++)
		{
			f[i][j]+=f[i][j-1];
		}
	}
	for(ll i=1;i<=maxn;i++)
	{
		for(ll j=0;j<=20;j++)f[i][j]+=f[i-1][j];
	}
} 
int main()
{
	init();
	ll T;
	scanf("%I64d",&T);
	while(T--)
	{
		ll n,m,p;
		scanf("%I64d%I64d%I64d",&n,&m,&p);
		ll h=min(n,m);
		ll pos;
		ll ans=0;
		if(p>=19)ans=(ll)n*m;
		else for(ll i=1;i<=h;i=pos+1)
		{
			pos=min(n/(n/i),m/(m/i));
			ans+=(f[pos][p]-f[i-1][p])*(ll)(n/i)*(ll)(m/i);
		}
		printf("%I64d\n",ans);
	}
	return 0;
}



hdu1695:

    b   d

求∑   ∑  [gcd(i,j)==k]*(f[j][i]=f[i][j]^=1)            (吐槽一下:我之前没看到约定a=c=1,不停的想特殊情况分类讨论,狂汗!!!)

   i=1 j=1

f[i][j]^=1是神马,不要问我,我只是装一下B,这样保证i,j这一对数只被访问一次。

                          b       d

通俗点,就是求∑     ∑ [gcd(i,j)==k] ,但要保证相同的(i,j)和(j,i)只能算一次。

                          i=1  j=1                                                                                     min(b/k,d/k)

主要是去重。 首先不管那个条件,就是裸的bzoj1101(或者说是problem b),∑   u(t)*[b/kt]*[d/kt]分块搞一下就好了。

                                                                                                                              t=1

然后我们把算重了的(i,j)去掉一次。

考虑[1,b],[1,d] 我们约定b<=d(可以swao嘛QAQ)

                                      b    b

那么,只有我们在计算∑   ∑ 的时候,会出现重复(会重复的是不重复的2倍,减去1/2就好了),当for j=b+1;j->d的时候不

                                     i=1 j=1                                 

会出现重复的i,j.所以我们只需要再计算i,j从[1,b],减去它的一半就好了。

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<iostream>
#include<queue>
#include<cmath>
#define ll long long
using namespace std;
/*
b   d
∑  ∑  [gcd(i,j)==k] 
i=a j=c 

转为f[b][d]+f[a-1][c-1]-f[a-1][d]-f[b][c-1]

        a   b                     a/k  b/k
f[a][b]=∑  ∑   [gcd(i,j)==k]  = ∑    ∑  [gcd(i,j)==1] 
        i=1 j=1                   i=1  j=1
        
        min(a',b')
f[a][b]=∑      u(t)  [a'/d][b'/d]
        d=1 
*/
const int maxn=100000+5;
int prime[maxn];
bool flag[maxn];
int tot;
int u[maxn];
int sum[maxn];
void init()
{
	u[1]=1;
	tot=0;
	for(int i=2;i<maxn;i++)
	{
		if(!flag[i])
		{
			prime[++tot]=i;
			u[i]=-1;
		} 
		for(int j=1;j<=tot&&i*prime[j]<maxn;j++)
		{
			flag[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				u[i*prime[j]]=0;
				break;
			}
			else u[i*prime[j]]=-u[i];
		}
	}
	for(int i=1;i<maxn;i++)sum[i]=sum[i-1]+u[i];
}
ll f(int a,int b,int k)
{
	a=a/k;
	b=b/k;
	int h=min(a,b);
	ll ans=0;
	int pos;
	for(int i=1;i<=h;i=pos+1)
	{
		pos=min(a/(a/i),b/(b/i));
		ans+=(ll)(sum[pos]-sum[i-1])*(a/i)*(b/i);//先把前面转成ll型防止爆炸 
	}
	return ans;
}
int a,b,c,d,k;
int main()
{
	init();
	int T;
	scanf("%d",&T);
	int cas=1;
	while(T--)
	{
		scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
		if(!k)
		{
			printf("Case %d: 0\n",cas++);
			continue;
		}
		if(b>d)swap(b,d);
		ll ans=0;
		ans=f(b,d,k)-f(b,b,k)/2;
		printf("Case %d: %I64d\n",cas++,ans);
	}
	return 0;
}



然后吐槽一下我犯低级错误的这道题QAQ

是一道polya的题。

思路见代码~~~~~~

最后我们推出ans=∑ n^(d-1)  phi(n/d)      然后我果断想到了线性筛,发现10亿(筛==GG)

                             d|n

然后吐槽一下我居然忘了phi的求法--》           phi(n)=n(1-1/p1)(1-1/p2)(1-1/p3)````(1-1/pk)

然后每次我们回答的时候就是直接求。复杂度X*sqrt(n)*log logn  ^2  (枚举i从1到sqrt(n),还有二分的优化,还有分解质因数也是log n级别)


复杂度这个我是大致估计的,没有严格证明,不过做了这么多类似的题感觉上应该是这个。QAQ

(不然它是怎么跑进2s的,3500多组数据呀!)

代码&&思路:

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<cstdlib>
#include<queue>
using namespace std;
/*f(d)表示最少转d能和自己重合的有多少种
//g(d)表示转d能和自己重合的有多少种
//g(n)=∑f(d)
      d|n
g(i)=n^i

ans=∑f(d) /d =∑  ∑ n^t *u(d/t)*1/d         为什么要/d? 因为对于一段d你把这一段最左边接到最右边 
   d|n        d|n t|d                         在环上是一样的,所以这d种会重复 
   
t|d|n  d=t*d'

 ∑      ∑     n^t * u(d')*1/d
t*d'|n  t|t*d'

 ∑ n^(t-1) phi(n/t)
t|n
*/
int prime[36000];
bool flag[36000];
int tot;
void init()
{
	tot=0;
	for(int i=2;i<=36000;i++)
	{
		if(!flag[i])prime[++tot]=i;
		for(int j=1;j<=tot&&i*prime[j]<=36000;j++)
		{
			flag[i*prime[j]]=1;
			if(i%prime[j]==0)break;
		}
	}
}
int qk(int x,int y,int mod)
{
	int res=1;
	x%=mod;
	while(y)
	{
		if(y&1)res=res*x%mod;
		y>>=1;
		x=x*x%mod;
	}
	return res%mod ;
}
int phi(int n,int p)
{
	int ans=n;
	for(int i=1;prime[i]*prime[i]<=n;i++)
	{
		if(n%prime[i]==0)
		{
			ans=ans-ans/prime[i];
		    while(n%prime[i]==0)n/=prime[i];
		}
		
	}
	if(n>1)ans=ans-ans/n;//话说我是怎么脑残把ans=ans-ans/n打成了  ans=ans=ans/n QAQ白白贡献几个WA 
	return ans%p;
}
void work(int n,int p)
{
	int ans=0;
	for(int i=1;i*i<=n;i++)
	{
		if(i*i==n)
		{
			ans+=qk(n,i-1,p)*phi(n/i,p);
			ans%=p;
		}
		else if(n%i==0)
		{
			ans+=qk(n,i-1,p)*phi(n/i,p)+qk(n,n/i-1,p)*phi(i,p);
			ans%=p;
		}
	}
	printf("%d\n",ans%p);
}
int n,p;
int main()
{
	init();
	/*马丹,又想到线性筛去了,弱爆了,都求到这一步了还不会算~~*/
	int T;
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d%d",&n,&p);
		work(n,p);
	}	
	return 0;
}



评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值