【noi.ac #2036】 乌拉乌拉

题目

http://noi.ac/contest/358/problem/2036

思路

显然 f(a, b) 必不可能等于 0，因为 i = p 1 时一定有 ai ≡ 1 ≡ bp 1(%p)。
定义小于质数 p 的正整数 a 关于 p 的阶为使得 ak ≡ 1(%p) 的最小的正整数 a，记
作 ord(a)。考虑如何求一个数 a 的阶。显然 %p 意义下的阶一定是 p 1 的约数，所以
我们可以对 p−1 质因数分解，然后求阶关于 p−1 的每个质因子的幂次。具体就是枚举
p−1 的每个质因子 qi，不妨设其对应的指数为 ui，然后从小到大枚举 qi 的指数 ti，如
果到了某个 ti 满足 a p 1 qui ti i ≡ 1(%p)，则表明 a 的阶在 qi 下的指数就是 ti。 记 c(p 1) 为 p 1 的不同质因子个数，那么求一次阶要枚举 c(p−1) 个质因子，对
于每个质因子 qi，都要用 O(log p) 的时间求出 a p 1 qui i 的值，而每次枚举指数可以 O(log qi)
转移，所以转移一个质因子的复杂度是 O(log qui i + log p) 的，这样暴力求一次阶的总
复杂度是 O(c(p−1)log p)。当然预处理 a p 1 qui i 时可以用分治来优化，复杂度可以降到
O(log c(p−1)log p)。具体分治方法如下：每次分治会求一个区间 [l, r] 内每个 i 对应
的 a p 1 qui i ，分治进这个区间前已经求得了 now = a∏i∈/[l,r] qui i , 分治进左右区间时左边代
入 now = now
∏i∈[mid+1,r] qui i , 右边代入 now = now
∏i∈[l,mid] qui i ，这样分治过程中每层都是
O(log P) 的，一共 O(log c(p 1)) 层。
令 r 为 p 的任一原根，a ≡ rx(%p), 1 < X < p 1, 则有 ord(a) = p 1
gcd(x,p 1)。令
a ≡ rx(%p), b ≡ ry(%p), 1 < x, y < p 1，则由裴蜀定理可得 ∃j > 0 满足 ai ≡ bj (%p)
当且仅当 gcd(y, p 1)|xi。那么就有
f(a, b) = gcd(lcm(x, y), p p 1)
gcd(x, p 1)
= p 1
gcd(x,p 1)
p 1
lcm(gcd(x,p 1),gcd(y,p 1))
= p 1
gcd(x,p 1)
gcd( p 1
gcd(x,p 1) , p 1
gcd(y,p 1) ) =
ord(a)
gcd(ord(a), ord(b))
那么对于一对 a, b, f(a, b) × f(b, a) = ord(a)ord(b)
gcd2(ord(a), ord(b))
，总答案就是
∑d ∑a ∑b [gcd(ord(a), ord(b)) = d]
ord(a)ord(b) d2
注意到这个过程本质上就是对每个质因子次数取 min 的卷积，因此可以通过类似 FWT
和高位前缀和的方法来实现。
使用 Pollard_Rho 算法分解质因数，总复杂度为 O(p 14 + n log c(p 1)log p + d(p
1)c(p 1))，其中 c(p 1) 为 p 1 质因子数量，d(p 1) 为 p 1 因子个数，可以验证
在 1018 范围内它们分别不超过 15 和 105000。

代码

#include<bits/stdc++.h>
using namespace std;
#define I inline int
#define V inline void
#define ld long double
#define ll long long int
#define FOR(i,a,b) for(int i=a;i<=b;i++)
#define ROF(i,a,b) for(int i=a;i>=b;i--)
const int N=2e5+1;
const int pri[9]={2,3,5,7,11,13,17,19,23};
int n,m,tot;
map<ll,ll>dp;
ll a[N],d[N],t[20],mod,ans;
ll lowbit(ll x){return x&-x;}
V check(ll&x){x-=mod,x+=x>>63&mod;}
ll mul(ll x,ll y,ll p){return(x*y-(ll)((ld)x/p*y)*p+p)%p;}
I test(int x,ll n){
	ll s=1,t=x,d=lowbit(n-1),m=n/d;
	for(;m;m>>=1,t=mul(t,t,n))if(m&1)s=mul(s,t,n);
	if(s==1||s==n-1)return 1;
	while(d>>=1){
		t=s,s=mul(s,s,n);
		if(s==1)return t==n-1;
	}
	return s==1;
}
I Mr(ll n){
	if(n==1)return 0;
	FOR(i,0,8)if(n==pri[i])return 1;
	FOR(i,0,8)if(!test(pri[i],n))return 0;
	return 1;
}
ll gcd(ll x,ll y){return!y?x:gcd(y,x%y);}
ll f(ll x,ll c,ll n){return(mul(x,x,n)+c)%n;}
V ins(ll x){FOR(i,1,m)if(t[i]==x)return;t[++m]=x;}
ll o(){
	static mt19937 seed(time(0));
	return abs((ll)seed());
}
ll once(ll n){
	FOR(i,0,8)if(n%pri[i]==0)return pri[i];
	ll x=0,y=0,c=o()%(n-1)+1,t;
	do x=f(x,c,n),y=f(f(y,c,n),c,n),t=gcd(abs(x-y),n);
	while(x!=y&&t==1);
	return t;
}
V Pr(ll n,ll p=0){
	if(n==1)return;
	if(Mr(n))return ins(n);
	for(p=n;p==n;p=once(n));
	while(n%p==0)n/=p;
	Pr(n),Pr(p);
}
V input(){
	scanf("%d%lld",&n,&mod);
	FOR(i,1,n)scanf("%lld",a+i);
}
ll Pow(ll t,ll x,ll s=1){
	for(;x;x>>=1,t=mul(t,t,mod))if(x&1)s=mul(s,t,mod);
	return s;
}
V dfs(ll x,int p){
	if(p>m)return void(d[++tot]=x);
	if(dfs(x,p+1),x%t[p]==0)dfs(x/t[p],p);
}
ll calc(ll x,ll s=mod-1){
	FOR(i,1,m)while(s%t[i]==0&&Pow(x,s/t[i])==1)s/=t[i];
	return s;
}
ll pw(ll x){return mul(x,x,mod);}
V init(){
	Pr(mod-1);
	dfs(mod-1,1),sort(d+1,d+1+tot);
	FOR(i,1,n)a[i]=calc(a[i]),check(dp[a[i]]+=a[i]);
}
V work(){
	FOR(i,1,m)ROF(j,tot,1)if(d[j]%t[i]==0)check(dp[d[j]/t[i]]+=dp[d[j]]);
	FOR(i,1,tot)dp[d[i]]=pw(dp[d[i]]);
	FOR(i,1,m)FOR(j,1,tot)if(d[j]%t[i]==0)check(dp[d[j]/t[i]]+=mod-dp[d[j]]);
	FOR(i,1,tot)check(ans+=Pow(pw(d[i]),mod-2,dp[d[i]]));
	cout<<ans;
}
int main(){
	input();
	init();
	work();
	return 0;
}