[NOI2016]循环之美

StaroForgin

已于 2022-01-24 15:41:34 修改

阅读量667

点赞数 1

分类专栏： ------数学------ # 反演文章标签：算法

于 2022-01-24 11:50:08 首次发布

本文链接：https://blog.csdn.net/Tan_tan_tann/article/details/122663451

版权

------数学------ 同时被 2 个专栏收录

76 篇文章 0 订阅

订阅专栏

反演

12 篇文章 0 订阅

订阅专栏

循环之美

题解

我们记 $n, m, K$ 对应的答案为 $f (n, m, K)$ ，显然有，
$f(n,m,K)=\sum_{x=1}^{n}\sum_{y=1}^{m}[(x,y)=1][y|K^k-1]\\ =\sum_{x=1}^{n}\sum_{y=1}^{m}[(x,y)=1][K^k\equiv 1\mod y]$ 由欧拉降幂可知，当且仅当 $(K, y) = 1$ 时，存在 $[K^k\equiv 1\mod y]$ ，所以该条件可以等价于 $[(K, y) = 1]$ ，有
$f(n,m,k)=\sum_{x=1}^n\sum_{y=1}^m[(x,y)=1][(K,y)=1]\\ =\sum_{x=1}^{n}\sum_{y=1}^{m}[(x,y)=1]\sum_{d|(y,K)}\mu(d)\\ =\sum_{d=1}^{K}\mu(d)\sum_{x=1}^n\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}[(x,dy)=1]\\ =\sum_{d|K}\mu(d)\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{x=1}^n[(y,x)=1][(d,x)=1]\\ =\sum_{d|K}\mu(d)f(\lfloor\frac{m}{d}\rfloor,n,d)$
显然是可以递归求解的，我们需要算一下的就是边界的情况。
首先， $f (0, m, k)$ 与 $f (n, 0, k)$ 都应该是 $0$ ，这是显然的，也就是前两维到边界的情况，避免其 $k$ 一直不变。
我们主要减少的其实还是 $k$ ，也就是说我们大部分边界情况都是 $f (n, m, 1)$ 的模样，考虑这怎么计算。
$f(n,m,1)=\sum_{x=1}^{n}\sum_{y=1}^{m}[(x,y)=1]\\ =\sum_{d=1}\mu(d)\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$ 显然可以通过整除分块加杜教筛快速求出，杜教筛是用来算 $\mu$ 的前缀和的。
由于我们这里要对两个数整除分块，杜教筛记忆化的标号可以建议使用 $\lfloor\frac{n}{d}\rfloor+\lfloor\frac{m}{d}\rfloor$ 来代替，因为这肯定是一个增函数，且我们整除分块每次正好会让其中一个改变。
我们记 $g(n)=\sum_{i=1}^n\mu(i)$ ，显然 $g(n)=1-\sum_{i=2}^{n}g(\lfloor\frac{n}{i}\rfloor)$ ，显然也可以整除分块。

时间复杂度我算不太来，当然杜教筛部分肯定是 $O\left(n^{\frac{2}{3}}\right)$ 的。

源码

#include<bits/stdc++.h>
using namespace std;
#define MAXN 1000005
#define MAXM 1000005
#define lowbit(x) (x&-x)
#define reg register
#define pb push_back
#define mkpr make_pair
#define fir first
#define sec second
#define lson (rt<<1)
#define rson (rt<<1|1)
typedef long long LL;
typedef unsigned long long uLL; 
typedef long double ld;
typedef pair<int,int> pii;
const int INF=0x3f3f3f3f;
const int mo=998244353;
const int mod=1e5+3;
const int inv2=5e8+4;
const int jzm=2333;
const int zero=2000;
const int n1=1000;
const int lim=1000000;
const int orG=3,ivG=332748118;
const long double Pi=acos(-1.0);
const double eps=1e-12;
template<typename _T>
_T Fabs(_T x){return x<0?-x:x;}
template<typename _T>
void read(_T &x){
	_T f=1;x=0;char s=getchar();
	while(s>'9'||s<'0'){if(s=='-')f=-1;s=getchar();}
	while('0'<=s&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=getchar();}
	x*=f;
}
template<typename _T>
void print(_T x){if(x<0){x=(~x)+1;putchar('-');}if(x>9)print(x/10);putchar(x%10+'0');}
int gcd(int a,int b){return !b?a:gcd(b,a%b);}
int add(int x,int y,int p){return x+y<p?x+y:x+y-p;}
void Add(int &x,int y,int p){x=add(x,y,p);}
int qkpow(int a,int s,int p){int t=1;while(s){if(s&1)t=1ll*t*a%p;a=1ll*a*a%p;s>>=1;}return t;}
int N,M,K,mobius[MAXM],pre[MAXM],prime[MAXM],cntp,g[MAXM*2];LL dp[MAXN];
bool oula[MAXM];
struct ming{
	int x,y,z;ming(){x=y=z=0;}
	ming(int X,int Y,int Z){x=X;y=Y;z=Z;}
	bool operator == (const ming &rhs)const{return x==rhs.x&&y==rhs.y&&z==rhs.z;}
	bool operator != (const ming &rhs)const{return x!=rhs.x||y!=rhs.y||z==rhs.z;}
	int Hash(){return (1ll*jzm*(1ll*x*jzm%mod+y)%mod+z)%mod;}
};
struct HashMap{
	ming val[MAXN];int head[MAXM],nxt[MAXN],tot;
	int query(ming x){
		int pos=x.Hash(),now=head[pos];
		while(now&&val[now]!=x)now=nxt[now];return now;
	}
	int insert(ming x){
		int pos=x.Hash(),now=++tot;val[now]=x;
		nxt[now]=head[pos];head[pos]=now;return now;
	}
}mp;
int work(int x){
	if(x<=lim)return pre[x];int y=N/x+M/x;
	if(~g[y])return g[y];g[y]=1;
	for(int l=2,r;l<=x;l=r+1)
		r=x/(x/l),g[y]-=1ll*(r-l+1)*work(x/l);
	return g[y];
}
LL sakura(int n,int m,int k){
	if(!n||!m)return 0;int id=mp.query(ming(n,m,k));
	if(id)return dp[id];id=mp.insert(ming(n,m,k));
	if(k==1){
		if(n>m)swap(n,m);
		for(int l=1,r;l<=n;l=r+1)
			r=min(n/(n/l),m/(m/l)),
			dp[id]+=1ll*(work(r)-work(l-1))*(n/r)*(m/r);;
		return dp[id];
	}
	for(int i=1;i*i<=k;i++)if(k%i==0){
		if(mobius[i])dp[id]+=1ll*mobius[i]*sakura(m/i,n,i);
		if(i!=k/i&&mobius[k/i])dp[id]+=1ll*mobius[k/i]*sakura(m/(k/i),n,k/i);
	}
	//printf("sakura %d %d %d\n",n,m,k,dp[id]);
	return dp[id];
}
void init(){
	mobius[1]=1;
	for(int i=2;i<=lim;i++){
		if(!oula[i])prime[++cntp]=i,mobius[i]=-1;
		for(int j=1;j<=cntp&&1ll*prime[j]*i<=lim;j++){
			oula[i*prime[j]]=1;
			if(i%prime[j]==0)continue;
			mobius[i*prime[j]]=-mobius[i];
		}
	}
	for(int i=1;i<=lim;i++)pre[i]=pre[i-1]+mobius[i];
}
signed main(){
	read(N);read(M);read(K);init();
	memset(g,-1,sizeof(g));
	printf("%lld\n",sakura(N,M,K));
	return 0;
}