P5282 【模板】快速阶乘算法 题解

· · 题解

改了一下第一篇题解的做法，不知道能不能在投一次题解。

发现首先根据 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\%，应该是并行计算了。

放下代码：

#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;
}