Loading... 改了一下第一篇题解的做法，不知道能不能在投一次题解。 发现首先根据 Wilson 定理： $$ \begin{aligned} (p-1) &\equiv -1 \pmod{p}\\ \prod_{i=1}^n i \prod_{i=1}^{p-1-n} (-i)&\equiv -1 \pmod{p}\\ n!(p-n-1)!(-1)^{p-n-1}&\equiv -1 \pmod{p} \end{aligned} $$ 于是可以选择 $n$ 和 $p-n-1$ 中较小者来算，这样已经卡掉一半常数了。 然后第一篇题解实际上跑了三秒多，于是加一次循环展开，用两个向量 AVX2 变量存，一次加上 16。这样时间会减少 $40\%$，应该是并行计算了。 放下代码： ```cpp #include<stdio.h> #include<immintrin.h> #define inl inline __attribute((always_inline)) #pragma GCC target("avx2") //__m256i为AVX2指令集的整型变量 static unsigned mod,r,n2_; static __m256i a0,mod1,R,hi32,ans1,ans2,ml1,ml2,ad; long long exgcd(long long a,long long b,long long &x,long long &y) { long long d=a; if(b==0) x=1,y=0; else d=exgcd(b,a%b,y,x),y-=a/b*x; return d; } inl unsigned mul(unsigned x,unsigned y)//蒙哥马利数域内的乘法及取模 { unsigned long long z=(unsigned long long)x*y;//计算乘积 return (z+(unsigned long long)(unsigned(z)*r)*mod)>>32;//对p取模 } inl __m256i add(__m256i _num1,__m256i _num2)//将8个32位整数相加并取模 { __m256i apb=_mm256_add_epi32(_num1,_num2),//将num1,num2内的8个整数对应相加 ret=_mm256_sub_epi32(apb,mod1);//将答案减去mod __m256i cmp=_mm256_cmpgt_epi32(a0,ret),//得到答案内小于0的数 add=_mm256_and_si256(cmp,mod1); return _mm256_add_epi32(add,ret);//将小于0的数加mod } inl __m256i mul(__m256i _num1,__m256i _num2)//将8个32位整数相乘并取模 { __m256i _num3=_num1,_num4,_num5=_num2; _num2=_mm256_mul_epu32(_num1,_num2);//将偶数位的整数相乘（得到4个64位整数并放回） _num1=_mm256_mul_epu32(_mm256_mul_epu32(_num2,R),mod1); _num4=_mm256_srli_epi64(_mm256_add_epi64(_num1,_num2),32);//每个64位整数均左移32位，空出高32位（奇数位） _num1=_mm256_srli_si256(_num3,4);_num2=_mm256_srli_si256(_num5,4);//进行4字节位移，即将奇数位的整数移动到偶数位 _num2=_mm256_mul_epu32(_num1,_num2);//将原奇数位的整数相乘 _num1=_mm256_mul_epu32(_mm256_mul_epu32(_num2,R),mod1); _num1=_mm256_and_si256(_mm256_add_epi64(_num1,_num2),hi32);//通过跟hi32做与取出低32位（偶数位） return _mm256_or_si256(_num1,_num4);//将偶数位乘积和奇数位乘积合并位一个AVX2变量返回 } inl unsigned mon_in(unsigned x){return mul(x,n2_);}//进入蒙哥马利数域 inl unsigned mon_out(unsigned x)//离开蒙哥马利数域 {unsigned ret=((x+(unsigned long long)(unsigned(x)*r)*mod)>>32);return ret<mod?ret:ret-mod;} inl int solve(int n,int p) { unsigned i=1; long long x,y; mod=p; n2_=-(unsigned long long)mod%mod;exgcd(mod,1ll<<32,x,y);r=-unsigned(x);//初始化一些蒙哥马利约减需要使用的变量 a0=_mm256_setzero_si256(),//预处理全0的AVX2变量（做加法有用） mod1=_mm256_set1_epi32(mod),R=_mm256_set1_epi32(r);//预处理全mod1,全r的AVX2变量 hi32=_mm256_set_epi32(-1,0,-1,0,-1,0,-1,0);//预处理奇数位全1的AVX2变量（做乘法有用） ans1=_mm256_set1_epi32(mon_in(1));//初始化ans ans2=_mm256_set1_epi32(mon_in(1));//初始化ans ad=_mm256_set1_epi32(mon_in(16));//初始化全16的AVX2变量 ml1=_mm256_set_epi32(mon_in(8),mon_in(7),mon_in(6),mon_in(5),mon_in(4),mon_in(3),mon_in(2),mon_in(1));//初始化因数 ml2=_mm256_set_epi32(mon_in(16),mon_in(15),mon_in(14),mon_in(13),mon_in(12),mon_in(11),mon_in(10),mon_in(9));//初始化因数 for(;i+16<=n;i+=16) { ans1=mul(ans1,ml1);//将ans的16位与因数的16位依次相乘 ans2=mul(ans2,ml2); ml1=add(ml1,ad);//将因数16位全部加16 ml2=add(ml2,ad); } unsigned *fl1=(unsigned*)&ans1; unsigned as1=mul(mul(mul(fl1[0],fl1[1]),mul(fl1[2],fl1[3])),mul(mul(fl1[4],fl1[5]),mul(fl1[6],fl1[7])));//取出ans的8位并相乘 unsigned *fl2=(unsigned*)&ans2; unsigned as2=mul(mul(mul(fl2[0],fl2[1]),mul(fl2[2],fl2[3])),mul(mul(fl2[4],fl2[5]),mul(fl2[6],fl2[7]))); unsigned as=mon_out(mul(as1,as2));//离开蒙哥马利数域 for(;i<=n;i++)//处理整段未做掉的末几位 as=1ull*as*i%mod; return as; } long long inv(long long n,long long p){ long long x,y; exgcd(n,p,x,y); x%=p; return x>=0?x:x+p; } int main() { int T; scanf("%d",&T); while(T--) { int n,p,res; scanf("%d%d",&n,&p); if(n<=p-1-n) res=solve(n,p); else{ res=inv(solve(p-1-n,p),p); if((p-n)&1) res=p-res; } printf("%d\n",res); } return 0; } ``` 最后修改：2024 年 12 月 02 日 © 允许规范转载 打赏 赞赏作者 赞 如果觉得我的文章对你有用，请随意赞赏