题目

给定一个质数ppp以及一个数列aia_iai，求：∑i=1n∑j=1nf(ai,aj)f(aj,ai)modp\sum_{i=1}^n\sum_{j=1}^nf(a_i,a_j)f(a_j,a_i) \mod p∑i=1n∑j=1nf(ai,aj)f(aj,ai)modp

其中f(x,y)f(x,y)f(x,y)为最小的iii满足∃j,xi=yj(modp)\exist j,x^i=y^j\pmod p∃j,xi=yj(modp)

n≤105n\le 10^5n≤105

p≤1018p\le 10^{18}p≤1018

正解

WC的时候见过类似的……然而这次还是没有做出来……其实主要的问题是不会求阶

看到这个东西首先要想到将每个aia_iai的阶求出来。

求之前先将写个Porllad-Rho将p−1p-1p−1给分解了。想到阶的定义：aaa的阶（记作ord(a)ord(a)ord(a)为最小的kkk满足ak≡1(modp)a^k\equiv 1 \pmod pak≡1(modp)，那么把ord(a)ord(a)ord(a)除以任意一个质因子，这个东西都不成立。

我们考虑分别求ord(a)ord(a)ord(a)的每个质因子的指数。设现在求质因子qqq，它在p−1p-1p−1中的指数为uiu_iui，于是我们要求出最小的ttt满足ap−1qiuqt≡1(modp)a^{\frac{p-1}{q^u_i}q^t}\equiv 1 \pmod paqiup−1qt≡1(modp)。

设x=ap−1qiux=a^{\frac{p-1}{q^u_i}}x=aqiup−1，暴力枚举ttt，每次x→xqx\to x^qx→xq（暴力快速幂），直到它第一次变成111为止。

算下时间复杂度：∣P∣lg⁡p+∑uilg⁡q=∣P∣lg⁡p+lg⁡p|P|\lg p+\sum u_i\lg q=|P|\lg p + \lg p∣P∣lgp+∑uilgq=∣P∣lgp+lgp。其中PPP为p−1p-1p−1的质因子集合。可以看到后面这一部分根本不是瓶颈。前面算xxx的那一部分就相对有点慢。（不过前面的那部分可以用个简单的分治来求，但是意义不大，所以就不说了）

接下来就是推式子时间：考虑求f(x,y)f(x,y)f(x,y)的值。设原根为ggg，则存在a,ba,ba,b满足ga≡x(modp)g^a\equiv x \pmod pga≡x(modp)，gb≡y(modp)g^b\equiv y \pmod pgb≡y(modp)。显然ord(x)=p−1gcd⁡(a,p−1)ord(x)=\frac{p-1}{\gcd(a,p-1)}ord(x)=gcd(a,p−1)p−1。

我们要找到最小的iii满足∃j,xi≡yj(modp)\exist j,x^i\equiv y^j \pmod p∃j,xi≡yj(modp)，即ai≡bj(modp−1)ai\equiv bj \pmod {p-1}ai≡bj(modp−1)。

由于jjj是任取的，用裴蜀定理可以得到gcd⁡(b,p−1)∣ai\gcd(b,p-1) | aigcd(b,p−1)∣ai，进而得到gcd⁡(b,p−1)gcd⁡(a,b,p−1)∣i\frac{\gcd(b,p-1)}{\gcd(a,b,p-1)}|igcd(a,b,p−1)gcd(b,p−1)∣i，由于要求最小的iii，所以取等号。

f(x,y)f(y,x)=gcd⁡(b,p−1)gcd⁡(a,p−1)gcd(a,b,p−1)2f(x,y)f(y,x)=\frac{\gcd(b,p-1)\gcd(a,p-1)}{gcd(a,b,p-1)^2}f(x,y)f(y,x)=gcd(a,b,p−1)2gcd(b,p−1)gcd(a,p−1)

=p−1ord(x)p−1ord(y)gcd⁡(p−1ord(x),p−1ord(y))2=\frac{\frac{p-1}{ord(x)}\frac{p-1}{ord(y)}}{\gcd(\frac{p-1}{ord(x)},\frac{p-1}{ord(y)})^2}=gcd(ord(x)p−1,ord(y)p−1)2ord(x)p−1ord(y)p−1

=lcm(ord(x),ord(y))2ord(x)ord(y)=\frac{lcm(ord(x),ord(y))^2}{ord(x)ord(y)}=ord(x)ord(y)lcm(ord(x),ord(y))2

=ord(x)ord(y)gcd⁡(ord(x),ord(y))2=\frac{ord(x)ord(y)}{\gcd(ord(x),ord(y))^2}=gcd(ord(x),ord(y))2ord(x)ord(y)

也就是求∑d∑i∑j[gcd⁡(ai,aj)=d]ord(ai)ord(aj)d2\sum_d \sum_i \sum_j [\gcd(a_i,a_j)=d]\frac{ord(a_i)ord(a_j)}{d^2}∑d∑i∑j[gcd(ai,aj)=d]d2ord(ai)ord(aj)

这个东西可以看成个高维的minminmin卷积，于是类似FWT、IFWT做高维前缀和，高维差分即可。

代码

using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <climits>
#include <cmath>
#include <cassert>
#include <cstdlib>
#include <ctime>
#define N 100010
#define ll long long
int n;
ll p;
ll a[N],r[N];
ll gcd(ll a,ll b){while (b){ll k=a%b;a=b,b=k;}return a;
}ll multi(ll x,ll y,ll mo){//x%=mo,y%=mo;ll z=(long double)x*y/mo;z=x*y-z*mo;if (z>=mo) z-=mo;else if (z<0) z+=mo;return z;
}
ll qpow(ll x,ll y=p-2,ll mo=p){x%=mo;ll r=1;for (;y;y>>=1,x=multi(x,x,mo))if (y&1)r=multi(r,x,mo);return r;
}
bool mr(ll v){static int p[10]={2,3,5,7,11,13,17,19,23,29};for (int i=0;i<10;++i){if (v==p[i]) return 1;if (v%p[i]==0) return 0;}ll x=v-1,y=0;for (;!(x&1);x>>=1,++y);for (int i=0;i<10;++i){ll t=qpow(p[i],x,v);if (t==1 || t==v-1) continue;for (int j=1;j<y;++j){t=multi(t,t,v);if (t==v-1) break;}if (t!=v-1) return 0;}return 1;
}
ll pr(ll n){while (1){ll a=rand()%n,b=a,c=rand()%n,s=1;for (int i=1,k=2;i;++i){a=(multi(a,a,n)+c)%n;if (a==b) break;ll t=multi(s,abs(a-b),n);if (t==0) return gcd(s,n);s=t;if (i==k || !(i&127)){ll g=gcd(s,n);if (g!=1) return g;if (i==k)k<<=1,b=a;}}}
}
ll q[70],nq,s[70],pro[70];
ll que[70],tail=0;
void divide(ll n){if (mr(n)){que[++tail]=n;return;}ll d=pr(n);divide(d);divide(n/d);
}
void initp(){divide(p-1);sort(que+1,que+tail+1);for (int i=1;i<=tail;++i){if (que[i]!=q[nq])q[++nq]=que[i];s[nq]++;}pro[0]=1;for (int i=1;i<=nq;++i)pro[i]=pro[i-1]*(s[i]+1);
}
ll ord[N],id[N];
void calco(ll x,ll &ord,ll &id){ord=1;for (int i=1;i<=nq;++i){ll d=qpow(x,(p-1)/qpow(q[i],s[i]));if (d==1) continue;ll qs=d;for (int t=1;t<=s[i];++t){qs=qpow(qs,q[i]);ord*=q[i];id+=pro[i-1];if (qs==1) break;}}
}
ll g[200000];
void fwt(){for (int i=1;i<=nq;++i)for (int j=0;j<pro[nq];j+=pro[i])for (int t=s[i];t>=1;--t)for (int k=0;k<pro[i-1];++k)(g[j+(t-1)*pro[i-1]+k]+=g[j+t*pro[i-1]+k])%=p;
}
void ifwt(){for (int i=1;i<=nq;++i)for (int j=0;j<pro[nq];j+=pro[i])for (int t=1;t<=s[i];++t)for (int k=0;k<pro[i-1];++k)(g[j+(t-1)*pro[i-1]+k]+=p-g[j+t*pro[i-1]+k])%=p;
}
ll ans;
void dfs(int x,ll prod,int num){if (x>nq){ans+=multi(g[num],qpow(prod,2*(p-2)),p);return;}for (int i=0;i<=s[x];++i,prod*=q[x],num+=pro[x-1])dfs(x+1,prod,num);
}
int main(){//freopen("in.txt","r",stdin);freopen("wlwl.in","r",stdin);freopen("wlwl.out","w",stdout);scanf("%d%lld",&n,&p);for (int i=1;i<=n;++i)scanf("%lld",&a[i]);srand(time(0));initp();for (int i=1;i<=n;++i){calco(a[i],ord[i],id[i]);(g[id[i]]+=ord[i])%=p; }fwt();for (int i=0;i<pro[nq];++i)g[i]=multi(g[i],g[i],p);ifwt();dfs(1,1,0);ans%=p;printf("%lld\n",ans);return 0;
}

6782. 2020.08.06【NOI2020】模拟T3 乌拉乌拉相关推荐

2020.08.06狂人日记：Python项目转C#项目问题
2020.08.06狂人日记:Python项目转C#项目问题 C#学习笔记问题及解决 C#学习笔记下拉选框中,在界面加入的元素集合和代码中写入的元素集合不会覆盖,即便有相同的元素也不会覆盖,代码中 ...
2020.08.08【NOIP提高组】模拟：奶牛的图片总结
2020.08.08[NOIP提高组]模拟:奶牛的图片总结 Description Farmer John希望给他的 N ( 1 ≤ N ≤ 100 , 000 ) N(1\leq N\leq100 ...
PYTHON学习笔记之（一）2020.08
PYTHON学习笔记之(一)2020.08 Python基础数据类型常见的列表.字典,以及元组.集合. 1 列表 list 1.1 列表转换字符串 stu = ['王一', '李二', '张三'] ...
（十三：2020.08.28）CVPR 2015 追踪之论文纲要（译）
CVPR 2020 追踪之论文纲要(修正于2020.08.27) 讲在前面论文目录讲在前面论坛很多博客都对论文做了总结和分类,但就医学领域而言,对这些论文的筛选信息显然需要更加精细的把控,所以自 ...
（十一：2020.08.28）CVPR 2017 追踪之论文纲要（译）
CVPR 2017 追踪之论文纲要(修正于2020.08.28) 讲在前面论文目录讲在前面论坛很多博客都对论文做了总结和分类,但就医学领域而言,对这些论文的筛选信息显然需要更加精细的把控,所以自 ...
2020年防爆电气模拟考试及防爆电气实操考试视频
题库来源:安全生产模拟考试一点通公众号小程序 2020年防爆电气模拟考试及防爆电气实操考试视频,包含防爆电气模拟考试答案和解析及防爆电气实操考试视频练习.由安全生产模拟考试一点通公众号结合国家防爆电气 ...
（十四：2020.08.28）CVPR 2014 追踪之论文纲要（译）
CVPR 2020 追踪之论文纲要(修正于2020.08.28) 讲在前面论文目录讲在前面论坛很多博客都对论文做了总结和分类,但就医学领域而言,对这些论文的筛选信息显然需要更加精细的把控,所以自 ...
大话西游2服务器维护,大话西游2经典版：2020年06月04日维护公告
??????? 各位大话西游玩家们大家好,又到了每周维护公告的介绍了.为了保证服务器的稳定和服务质量,<大话西游2经典版>将于2020年06月04日(本周四)早上8:00停机,进行每周例行 ...
网络安全与黑客技术学习资源汇总---2020.08更新
网络安全与黑客攻防学习资源汇总,截至 2020.08.15 可正常访问 (目前不可访问的网站未列出). 文章目录漏洞测试学习平台安全资讯平台安全行业协会/组织 POC(验证性测试)提交学习黑 ...

6782. 2020.08.06【NOI2020】模拟T3 乌拉乌拉

题目

正解

代码

6782. 2020.08.06【NOI2020】模拟T3 乌拉乌拉相关推荐

最新文章

热门文章