莫比乌斯反演 | OceanEye's Blog

@OceanEye7年前

04/13

21:39

BZOJ4816 莫比乌斯反演

OceanEye 撰写 BZOJ 莫比乌斯反演无回复

nlogn狄利克雷卷积打表

和普通的差不多，就是和换成积，贡献算一算就可以了

然后死活想不起来写 O(n根号n) 的TLE到死

安利金策大爷的教程

贴完代码跑路qwq

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define ll long long
#define p(x)  ('0'<=x&&x<='9')
char cc;
void read(int &x)
{
	cc=getchar(); x=0;
	while(!p(cc)) cc=getchar();
	while(p(cc)) { x=x*10+cc-48; cc=getchar(); }
}

using namespace std;
#define M 1000000
#define mod 1000000007

void ext__(int a,int b,int &x,int &y)
{
	if(b==0)
	{
		x=1; y=0;
		return;
	}
	ext__(b,a%b,y,x);
	y-=(a/b)*x;
}

int pow_mod(int a,ll b)
{
	int ret=1;
	while(b)
	{
		if(b&1) ret=(ll)ret*a%mod;
		a=(ll)a*a%mod;
		b>>=1;
	}
	return ret;
}

int get_inv(int k)
{
	if(k==1) return 1;
	int x,y;
	ext__(mod,k,x,y);
	return (mod+(y%mod))%mod;
}

int prime[M+10],miu[M+10],_[M+10],pre[M+10],fib[M+10],tot;
bool isprime[M+10];
void init()
{
	fib[0]=0; fib[1]=1;
	for(int i=2;i<=M;++i)
	{
		fib[i]=fib[i-1]+fib[i-2];
		if(fib[i]>mod) fib[i]-=mod;
	}
	miu[1]=1; _[1]=1;
	for(int i=2;i<=M;i++)
	{
		if(!isprime[i])
		{
			prime[tot++]=i;
			miu[i]=-1;
			_[i]=fib[i];
		}
		for(int j=0;j<tot;j++)
		{
			if(i*prime[j]>M) break;
			isprime[i*prime[j]]=true;
			if(i%prime[j]==0)
			{
				miu[i*prime[j]]=0;
				break;
			}
			miu[i*prime[j]]=-miu[i];
		}
	}
	for(int i=0;i<=M;i++) pre[i]=1;
	for(int i=1;i*i<=M;i++)
	{
		for(int j=i;j*i<=M;j++)
			if(i==j)
			{
				if(miu[i]==-1)
				pre[i*j]=(ll)pre[i*j]*get_inv(fib[i])%mod;
				if(miu[i]==1)
				pre[i*j]=(ll)pre[i*j]*fib[i]%mod;
			}
			else
			{
				if(miu[i]==-1)
					pre[i*j]=(ll)pre[i*j]*get_inv(fib[j])%mod;
				if(miu[i]==1)
					pre[i*j]=(ll)pre[i*j]*fib[j]%mod;
				if(miu[j]==-1)
					pre[i*j]=(ll)pre[i*j]*get_inv(fib[i])%mod;
				if(miu[j]==1)
					pre[i*j]=(ll)pre[i*j]*fib[i]%mod;
			}
	}
	for(int i=2;i<=M;i++) pre[i]=(ll)pre[i]*pre[i-1]%mod;
}
int T,A,B,ans;
int main()
{
	init();
	for(read(T);T;T--)
	{
		ans=1;
		read(A); read(B);
		if(A>B) swap(A,B);
		for(int node=1,m;node<=A;node=m+1)
		{
			m=min(A/(A/node),B/(B/node));
			ans=(ll)ans*pow_mod((ll)pre[m]*get_inv(pre[node-1])%mod,(ll)(A/node)*(B/node))%mod;
		}
		printf("%d\n",ans);
	}
	return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

#include <cstdio>

#include <algorithm>

#include <cstring>

#include <cmath>

#define ll long long

#define p(x) ('0'<=x&&x<='9')

char cc;

void read(int &x)

{

cc=getchar(); x=0;

while(!p(cc)) cc=getchar();

while(p(cc)) { x=x*10+cc-48; cc=getchar(); }

}

using namespace std;

#define M 1000000

#define mod 1000000007

void ext__(int a,int b,int &x,int &y)

{

if(b==0)

{

x=1; y=0;

return;

}

ext__(b,a%b,y,x);

y-=(a/b)*x;

}

int pow_mod(int a,ll b)

{

int ret=1;

while(b)

{

if(b&1) ret=(ll)ret*a%mod;

a=(ll)a*a%mod;

b>>=1;

}

return ret;

}

int get_inv(int k)

{

if(k==1) return 1;

int x,y;

ext__(mod,k,x,y);

return (mod+(y%mod))%mod;

}

int prime[M+10],miu[M+10],_[M+10],pre[M+10],fib[M+10],tot;

bool isprime[M+10];

void init()

{

fib[0]=0; fib[1]=1;

for(int i=2;i<=M;++i)

{

fib[i]=fib[i-1]+fib[i-2];

if(fib[i]>mod) fib[i]-=mod;

}

miu[1]=1; _[1]=1;

for(int i=2;i<=M;i++)

{

if(!isprime[i])

{

prime[tot++]=i;

miu[i]=-1;

_[i]=fib[i];

}

for(int j=0;j<tot;j++)

{

if(i*prime[j]>M) break;

isprime[i*prime[j]]=true;

if(i%prime[j]==0)

{

miu[i*prime[j]]=0;

break;

}

miu[i*prime[j]]=-miu[i];

}

for(int i=0;i<=M;i++) pre[i]=1;

for(int i=1;i*i<=M;i++)

{

for(int j=i;j*i<=M;j++)

if(i==j)

{

if(miu[i]==-1)

pre[i*j]=(ll)pre[i*j]*get_inv(fib[i])%mod;

if(miu[i]==1)

pre[i*j]=(ll)pre[i*j]*fib[i]%mod;

}

else

{

if(miu[i]==-1)

pre[i*j]=(ll)pre[i*j]*get_inv(fib[j])%mod;

if(miu[i]==1)

pre[i*j]=(ll)pre[i*j]*fib[j]%mod;

if(miu[j]==-1)

pre[i*j]=(ll)pre[i*j]*get_inv(fib[i])%mod;

if(miu[j]==1)

pre[i*j]=(ll)pre[i*j]*fib[i]%mod;

}

for(int i=2;i<=M;i++) pre[i]=(ll)pre[i]*pre[i-1]%mod;

}

int T,A,B,ans;

int main()

{

init();

for(read(T);T;T--)

{

ans=1;

read(A); read(B);

if(A>B) swap(A,B);

for(int node=1,m;node<=A;node=m+1)

{

m=min(A/(A/node),B/(B/node));

ans=(ll)ans*pow_mod((ll)pre[m]*get_inv(pre[node-1])%mod,(ll)(A/node)*(B/node))%mod;

}

printf("%d\n",ans);

}

return 0;

}

BZOJ4816 莫比乌斯反演

@OceanEye7年前

04/11

23:03

BZOJ2820

OceanEye 撰写莫比乌斯反演无回复

莫比乌斯反演

普通的莫比乌斯反演就是搞出一个分式然后利用[n/i]的连续性制造根号复杂度的东西啊= =

这里的推导也是一样的

先找到simple的f(x) & g(x)

f(x)=sigma(i)sigma(j) gcd(i,j)==x

g(x)=sigma(i)sigma(j) x|gcd(i,j)

上面这个东西随便反演吧= =然后把所有的质数尝试带进x里面

sigma(prime -> p) f(x) == sigma(prime -> p) sigma(dx<=p) g(dx) * μ(d)

把dx换成C，会发现一个C对应着几个 μ(d)，这里是一个质因子对应一个 μ(d)。

这里把上面那行东西的值想成一个奇奇gaygay的数论函数h(C)就可以了。

【最神奇的是这个函数居然可以线性求= =太神了不仔细看都看不出来】

这里就可以对于1-N中的每个数来计算贡献了= =

即是 sigma(C = 1 to N) [N/C]*[M/C] * h(C)

然后套一套根号前缀和就可以了。

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define p(x) ('0'<=x&&x<='9')
#define ll long long
char cc;
void read(int &x)
{
	cc=getchar(); x=0;
	while(!p(cc)) cc=getchar();
	while(p(cc)) { x=x*10+cc-48; cc=getchar(); }
}

using namespace std;
#define M 10000000
int prime[M+10],miu[M+10],mn[M+10],tot;
int g[M];//pre[];
bool isprime[M+10];

void pre()
{
	miu[1]=1;
	for(int i=2;i<=M;i++)
	{
		if(!isprime[i])
		{
			miu[i]=-1;g[i]=1;
			prime[tot++]=i;
		}
		for(int j=0,k;j<tot;j++)
		{
			k=i*prime[j];
			if(k>M) break;
			isprime[k]=true;
			if(i%prime[j]==0)
			{
				miu[k]=0;
				g[k]=miu[i];
				break;
			}
			miu[k]=-miu[i];
			g[k]=miu[i]-g[i];
		}
	}
	for(int i=1;i<=M;i++) g[i]+=g[i-1];
}
int T,A,B;
ll ans;
int main()
{
	pre();
	for(read(T);T--;)
	{
		ans=0;
		read(A); read(B);
		if(A>B) swap(A,B);
		for(int node=1,m;node<=A;node=m+1)
		{
			m=min(A/(A/node),B/(B/node));
			ans+=(ll)(g[m]-g[node-1])*(A/m)*(B/m);
		}
		printf("%lld\n",ans);
	}
	return 0;
}

#include <cstdio>

#include <algorithm>

#include <cstring>

#include <cmath>

#define p(x) ('0'<=x&&x<='9')

#define ll long long

char cc;

void read(int &x)

{

cc=getchar(); x=0;

while(!p(cc)) cc=getchar();

while(p(cc)) { x=x*10+cc-48; cc=getchar(); }

}

using namespace std;

#define M 10000000

int prime[M+10],miu[M+10],mn[M+10],tot;

int g[M];//pre[];

bool isprime[M+10];

void pre()

{

miu[1]=1;

for(int i=2;i<=M;i++)

{

if(!isprime[i])

{

miu[i]=-1;g[i]=1;

prime[tot++]=i;

}

for(int j=0,k;j<tot;j++)

{

k=i*prime[j];

if(k>M) break;

isprime[k]=true;

if(i%prime[j]==0)

{

miu[k]=0;

g[k]=miu[i];

break;

}

miu[k]=-miu[i];

g[k]=miu[i]-g[i];

}

for(int i=1;i<=M;i++) g[i]+=g[i-1];

}

int T,A,B;

ll ans;

int main()

{

pre();

for(read(T);T--;)

{

ans=0;

read(A); read(B);

if(A>B) swap(A,B);

for(int node=1,m;node<=A;node=m+1)

{

m=min(A/(A/node),B/(B/node));

ans+=(ll)(g[m]-g[node-1])*(A/m)*(B/m);

}

printf("%lld\n",ans);

}

return 0;

}

OceanEye's Blog

很多人即使只见过一面，已经算见过了最后一面。

@OceanEye7年前

BZOJ4816 莫比乌斯反演

BZOJ4816 莫比乌斯反演

@OceanEye7年前

BZOJ2820

BZOJ2820

友情链接

近期文章

标签