X-ray（NTT）

最新推荐文章于 2025-12-18 16:28:22 发布

原创最新推荐文章于 2025-12-18 16:28:22 发布 · 1.4k 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#题解 #c++

题解专栏收录该内容

348 篇文章

订阅专栏

文章讲述了如何解决一道数论题目，该题目涉及计算数组中每个元素欧拉函数的倒数之积的期望值。在宇宙射线导致数组元素变化后，需要求解新结果的期望。解题策略包括使用快速数论变换（NTT）来处理涉及欧拉函数的复杂计算，以及考虑时间复杂度为O(nlogn)的方法。代码示例展示了具体的实现过程。

题目描述

小 $D$ 最近在复习数论，他遇到一道题：

给你一个长度为 $n$ 的数组 $a$ ，要你求出每个 $a_i$ 的欧拉函数的倒数的积，即求 $\prod\limits_{i=1}^n\dfrac{1}{\phi(a_i)}$
其中 $1\leq n,a_i\leq 10^6$

小 $D$ 刚打完，正准备交，突然，一道宇宙射线射向他的电脑，改变了一部分 $a_i$ 。具体来说，宇宙射线随机删除了两个不同的数 $a_i,a_j$ ，然后又往 $a$ 数组中添加了 $a_i+a_j$ 。小 $D$ 需要你求出最后程序输出结果的期望。输出答案模 $998244353$ 后的值。

输入格式

第一行一个正整数 $n$ 。

第二行 $n$ 个正整数 $a_i$ 。

输出格式

一行一个数，表示最后程序输出结果的期望对 $998244353$ 取模后的值。

样例输入

4
2 3 4 5

样例输出

490456861

数据范围

对于 $20\%$ 的数据， $1\leq n\leq 10^4$
对于 $100\%$ 的数据， $1\leq n,a_i\leq 10^6$

题解

前置知识：快速数论变换（ $NTT$ ）

首先，令 $t=\prod\limits_{i=1}^n\dfrac{1}{\phi(a_i)}$ ，则 $t$ 可以 $O (n)$ 求出。

那么，我们要求的就是

$t\times[(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^n\dfrac{\phi(a_i)\cdot \phi(a_j)}{\phi(a_i+a_j)})-\sum\limits_{i=1}^n\dfrac{\phi(a_i)\cdot \phi(a_i)}{\phi(2a_i)}]$

中括号内的后一个部分也可以 $O (n)$ 求出，现在要求的就是前一个部分了。

我们可以将 $\sum\limits_{i=1}^{n}\sum\limits_{j=1}^n\dfrac{\phi(a_i)\cdot \phi(a_j)}{\phi(a_i+a_j)}$ 变为

$\sum\limits_{k=1}^{mx}\dfrac{1}{\phi(k)}\sum_{a_i+a_j=k}\phi(a_i)\cdot \phi(a_j)$

那么后面这部分就可以用 $NTT$ 来求了。

时间复杂度为 $O(n\log n)$ 。

code

#include<bits/stdc++.h>
using namespace std;
const int N=2200000;
int n,mx=0,z[N+5],p[N+5],ph[N+5];
long long w,wn,ans=1,tmp,jc[N+5],ny[N+5],f[N+5];
const long long G=3,mod=998244353;
long long mi(long long t,long long v){
	if(!v) return 1;
	long long re=mi(t,v/2);
	re=re*re%mod;
	if(v&1) re=re*t%mod;
	return re;
}
void ch(long long *a,int l){
	for(int i=1,j=l/2;i<l-1;i++){
		if(i<j) swap(a[i],a[j]);
		int k=l/2;
		while(j>=k){
			j-=k;k>>=1;
		}
		j+=k;
	}
}
void ntt(long long *a,int l,int fl){
	for(int i=2;i<=l;i<<=1){
		if(fl==1) wn=mi(G,(mod-1)/i);
		else wn=mi(G,mod-1-(mod-1)/i);
		for(int j=0;j<l;j+=i){
			w=1;
			for(int k=j;k<j+i/2;k++,w=w*wn%mod){
				long long t=a[k],u=w*a[k+i/2]%mod;
				a[k]=(t+u)%mod;
				a[k+i/2]=(t-u+mod)%mod;
			}
		}
	}
	if(fl==-1){
		long long ny=mi(l,mod-2);
		for(int i=0;i<l;i++) a[i]=a[i]*ny%mod;
	}
}
void init(){
	ph[1]=1;
	for(int i=2;i<=N;i++){
		if(!z[i]){
			p[++p[0]]=i;
			ph[i]=i-1;
		}
		for(int j=1;j<=p[0]&&i*p[j]<=N;j++){
			z[i*p[j]]=1;
			if(i%p[j]==0){
				ph[i*p[j]]=ph[i]*p[j];
				break;
			}
			ph[i*p[j]]=ph[i]*(p[j]-1);
		}
	}
	jc[0]=1;
	for(int i=1;i<=N;i++) jc[i]=jc[i-1]*i%mod;
	ny[N]=mi(jc[N],mod-2);
	for(int i=N-1;i>=0;i--) ny[i]=ny[i+1]*(i+1)%mod;
}
int main()
{
	init();
	scanf("%d",&n);
	for(int i=1,x;i<=n;i++){
		scanf("%d",&x);
		mx=max(mx,x);
		ans=ans*ny[ph[x]]%mod*jc[ph[x]-1]%mod;
		f[x]=(f[x]+ph[x])%mod;
		tmp=(tmp-1ll*ph[x]*ph[x]%mod*ny[ph[2*x]]%mod*jc[ph[2*x]-1]%mod+mod)%mod;
	}
	ans=ans*mi(1ll*n*(n-1)%mod,mod-2)%mod;
	int len=1;
	while(len<2*(mx+1)) len<<=1;
	ch(f,len);
	ntt(f,len,1);
	for(int i=0;i<len;i++) f[i]=f[i]*f[i]%mod;
	ch(f,len);
	ntt(f,len,-1);
	for(int i=1;i<len;i++){
		tmp=(tmp+f[i]*ny[ph[i]]%mod*jc[ph[i]-1]%mod)%mod;
	}
	ans=ans*tmp%mod;
	printf("%lld",ans);
	return 0;
}