数论模板

快速幂&取余运算

$O(\log_2{b})$ 求 $a^b \operatorname{mod} p$。

ll qp(ll a,ll b,ll p)
{
    ll ans=1;
    while(b)
    {
        if(b&1)ans=(ans*a)%p;
        a=(a*a)%p;
        b>>=1;
    }
    return ans;
}

乘法逆元

快速幂

模数是质数版：qp(a,p-2,p)

模数仅满足 $\gcd(a,p)=1$ 版：qp(a,phi(p)-1,p)

exGCD版

如下：

ll inv(ll a,ll p)
{
    ll tx,ty;exgcd(a,p,tx,ty);
    return (tx%p+p)%p;
}

线性预处理

比较懒，直接贺了 oiwiki 的代码

inv[1] = 1;
for (int i = 2; i <= n; ++i) {
  inv[i] = (long long)(p - p / i) * inv[p % i] % p;
}

康托展开

排列$a_i$的排名$num$为

$$num=1+\Sigma^n_{i=1}s_i(n-i)!$$

其中$s_i=\Sigma^n_{j=i}[a_j<a_i]$

int n,nl[maxn],tl[maxn];ll fc=1,ans=1;
void upd(int dc){for(;dc<=n;dc+=dc&(-dc))tl[dc]++;}
int ask(int dc){int ans=0;for(;dc;dc-=dc&(-dc))ans+=tl[dc];return ans;}
int main()
{
    n=in();
    fru(i,1,n)
    {
        nl[i]=in();
    }
    frd(i,n,1)
    {
        upd(nl[i]);
        (ans+=fc*ask(nl[i]-1))%=mod;
        (fc*=n-i+1)%=mod;
    }
    printf("%lld",ans);
    return 0;
}

筛法

埃氏筛

模板谁都会写，不放了。

埃氏筛时间复杂度为 $O(n\log\log{n})$，线性筛为 $O(n)$。

通常不用，但是可以从中得出一个推论：枚举 $i,j(i,j\in\Z^+)$ 使得 $i*{prime}_j\leqslant n$ 的复杂度为 $O(n\log\log n)$。

线性筛（欧拉筛）

注：只是筛一个数的质因子只要从小到大 $O(\sqrt n)$ 得筛就行了，因为合数 $nm$ 肯定会先被 $n$ 和 $m$ 筛掉，一点不剩，所以最后出来的都是质数。

bool p[/*值域*/];
int prime[/*值域*/>>1],tot,n;
inline void getprime(void)
{
    p[1]=true;
    for(int i=2;i<=n;i++)
    {
        if(!p[i])prime[++tot]=i;
        for(int j=1;j<=tot&&prime[j]*i<=n;j++)
        {
            p[prime[j]*i]=true;
            if(i%prime[j]==0)break;
        }
    }
}

线筛计算欧拉函数

欧拉函数：$\varphi(x)$ 表示 $\leqslant{x}$ 且与 $x$ 互质的正整数个数，如：$\varphi(1)=1$，$\varphi(p)=p-1$。

欧拉定理：当 $\gcd(a,m)=1$ 时有 $a^{\varphi(m)}\equiv 1$。

拓展欧拉定理：

$$a^b \equiv \begin{cases} a^{b \bmod \varphi(m)}, &\gcd(a,m) = 1, \\ a^b, &\gcd(a,m)\ne 1, b < \varphi(m), \\ a^{(b \bmod \varphi(m)) + \varphi(m)}, &\gcd(a,m)\ne 1, b \ge \varphi(m). \end{cases} \pmod m$$

bool bl[/*值域*/];
int pme[/*值域*/>>1],tot;
ll phi[/*值域*/],md;
void getphi(void)
{
    bl[1]=true;phi[1]=1;
    fru(i,2,/*值域*/)
    {
        if(!bl[i])
        {
            pme[++tot]=i;phi[i]=i-1;
        }
        for(register int j=1;i*pme[j]<=/*值域*/&&j<=tot;j++)
        {
            bl[i*pme[j]]=true;
            if(i%pme[j])
            {
                phi[i*pme[j]]=phi[i]*phi[pme[j]];
            }
            else
            {
                phi[i*pme[j]]=phi[i]*pme[j];
                break;
            }
        }
    }
}

Pollard rho质因数分解

依赖于米勒罗宾质数判断，复杂度 $O(n^{\frac{1}{4}})$。

mt19937 rnd(time(0));
ll al[]={2,3,5,7,11,13,17,19,23,29,31,37};
inline ll poww(ll a,ll b,ll p)
{
    ll ans=1;
    while(b)
    {
        if(b&1)ans=((__int128)ans*a)%p;
        a=((__int128)a*a)%p;
        b>>=1;
    }
    return ans;
}
bool mr(ll x)
{
    if(x<3||x%2==0)return x==2;
    ll t1=x-1,t2=0,tmp;register int j;
    while(!(t1&1)){t1>>=1;t2++;}
    fru(i,0,11)
    if(x%al[i])
    {
        tmp=poww(al[i],t1,x);
        if(tmp==1)continue;
        for(j=0;j<t2;j++)
        {
            if(tmp==x-1)break;
            tmp=((__int128)tmp*tmp)%x;
        }
        if(j==t2)return false;
    }
    return true;
}
ll pr(ll x)
{
    ll c=rnd()%(x-2)+1,s=0,t=0,val,tmp;
    for(register int lin=1;;lin<<=1,s=t,val=1)
    {
        fru(j,1,lin)
        {
            t=((__int128)t*t+c)%x;
            val=(__int128)val*abs(s-t)%x;
            if(lin%127==0)
            {
                tmp=__gcd(val,x);
                if(tmp>1)return tmp;
            }
        }
        tmp=__gcd(val,x);
        if(tmp>1)return tmp;
    }
}
void fac(ll x,ll& ans)
{
    if(x==1)return ;
    if(mr(x))
    {
        ans=max(ans,x);
        return ;
    }
    ll tmp=pr(x);
    while(!(x%tmp))
    x/=tmp;
    fac(tmp,ans);fac(x,ans);
}

二次剩余

Cipolla

mt19937 rnd((unsigned)time(NULL));
ll qp(ll a,ll b,ll p)
{
    ll ans=1;
    while(b)
    {
        if(b&1)(ans*=a)%=p;
        (a*=a)%=p;
        b>>=1;
    }
    return ans;
}
ll pwi,mod;
struct cpl
{
    ll a,b;
    cpl operator*(cpl y)
    {
        return (cpl){(a*y.a+pwi*(b*y.b%mod))%mod,(a*y.b+b*y.a)%mod};
    }
    bool operator==(cpl y)
    {
        return  a==y.a&&b==y.b;
    }
};
cpl qp(cpl a,ll b)
{
    cpl ans=(cpl){1,0};
    while(b)
    {
        if(b&1)ans=ans*a;
        a=a*a;
        b>>=1;
    }
    return ans;
}
ll cip(ll a,ll p)
{
    if(qp(a,(p-1)>>1,p)==p-1)return -1;
    mod=p;
    ll r;
    while(1)
    {
        r=rnd()%p;
        pwi=((r*r-a)%p+p)%p;
//        cout<<r<<endl;
        if(qp(pwi,(p-1)>>1,p)==p-1)break;
    }
    return qp((cpl){r,1},(p+1)>>1).a;
}

最大公约数（GCD）

递归版：

int gcd(int a,int b)
{
    return b?gcd(b,a%b):a;
}

非递归版：

int gcd(int a,int b)
{
    int tm;
    while(b)
    {
        tm=a%b;
        a=b;
        b=tm;
    }
    return a;
}

拓展欧几里得（exGCD）

根据贝祖定理，不定方程 $ax+by=c$ 当且仅当 $\gcd(a,b)|c$ 的时候有解。

因此只需对 $ax+by=\gcd(a,b)$ 求解即可。

以下代码在 $O(\log n)$ 的时间求出 $\gcd(a,b)$ 并且给出了一组特解 $x_0$，$y_0$。

之后即可使用

$$\begin{cases} x=x_0+\cfrac {kb}{\gcd(a,b)} \\ y=y_0-\cfrac {ka}{\gcd(a,b)} \end{cases}(k\in\Z)$$

得到其他所有解。

int exgcd(int a,int b,int& x,int& y)
{
    if(!b)
    {
        x=1;y=0;
        return a;
    }
    int td=exgcd(b,a%b,x,y),tx;tx=x;
    x=y;
    y=(tx-a/b*y);
    return td;
}

中国剩余定理（CRT）

太简单了，所以……

拓展中国剩余定理（exCRT）

ll exgcd(ll a,ll b,ll& x,ll& y)
{ 
    if(!b)
    {
        x=1;y=0;
        return a;
    }
    ll tp=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return tp;
}
ll excrt(int n,ll* al,ll* bl)
{
    ll md,tp,ans,tx,ty;
    ans=al[1];md=bl[1];
    fru(i,2,n)
    {
        tp=exgcd(md,bl[i],tx,ty);
        if((al[i]-ans)%tp!=0)return -1;
        tx=((__int128)((al[i]-ans)/tp)*tx)%(md*(bl[i]/tp));
        (ans+=(__int128)tx*md%(md*=(bl[i]/tp)))%=md;
    }
    return (ans%md+md)%md;
}

离散对数/大步小步算法（BSGS）

以下代码拥有 $O(\sqrt{p}\log{\lceil\sqrt{p}\rceil})$ 的时间复杂度。

如果题目卡最后的 $\log$，就换成unordered_map，运气好可以 $O(1)$。

不过考虑到unordered_map使用哈希表实现，运气不好会出事，平时不卡就用map吧。

ll bsgs(ll a,ll b,ll p)//仅在a,p互质时有用，无解返回-1
{
    map<int,int>mp;
    ll t=sqrt(p-1)+1,dc=1,tp=a;
    if(b%p==1)return 0;
    fru(i,0,t-1)
    {
        mp[b*dc%p]=i;
        (dc*=tp)%=p;
    }
    dc=1;tp=qp(a,t,p);
    fru(i,1,t)
    {
        (dc*=tp)%=p;
        auto rs=mp.find(dc);
        if(rs!=mp.end())return i*t-(*rs).second;
    }
    return -1;
}

拓展 BSGS/exBSGS

细节真的超级多😢 。

ll exgcd(ll a,ll b,ll& x,ll& y)
{
    if(!b)
    {
        x=1;y=0;
        return a;
    }
    ll tp=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return tp;
}
ll inv(ll a,ll p)
{
    ll tx,ty;
    exgcd(a,p,tx,ty);
    return (tx%p+p)%p;
}
ll qp(ll a,ll b,ll p)
{
    ll ans=1;
    while(b)
    {
        if(b&1)(ans*=a)%=p;
        (a*=a)%=p;
        b>>=1;
    }
    return ans;
}
ll bsgs(ll a,ll b,ll p)
{
    // cout<<a<<' '<<b<<' '<<p<<endl;
    map<int,int>mp;
    ll t=sqrt(p-1)+1,dc=1,tp=a;
    if(b%p==1||p==1)return 0;//注意为0有n^0==1和p==1两个情况
    fru(i,0,t-1)
    {
        mp[b*dc%p]=i;
        (dc*=tp)%=p;
    }
    dc=1;tp=qp(a,t,p);
    fru(i,1,t)
    {
        (dc*=tp)%=p;
        auto rs=mp.find(dc);
        if(rs!=mp.end())return i*t-(*rs).second;
    }
    return -1;
}
ll exbsgs(ll a,ll b,ll p)
{
    ll ad=1,k=0,mod=p,tb=b,tp,dc=1;//防止b的变化，新开tb
    while((tp=__gcd(mod,a))!=1)
    {
        // cout<<dc<<" "<<b<<' '<<(dc-b)%p<<endl;
        if((dc-b)%p==0)return k;//在<k就有解时可以避免向下迭代，从而有解。
        if(tb%tp)return -1;
        mod/=tp;tb/=tp;
        ad=ad*(a/tp)%mod;
        // cout<<ad<<endl;
        (dc*=a)%=p;
        k++;
    }
    // cout<<b<<' '<<ad<<' '<<inv(ad,mod)<<endl;
    tp=bsgs(a,tb*inv(ad,mod)%mod,mod);
    return tp==-1?-1:(tp+k);
}

FFT/虚数

吐槽：最后的代码真的只有FFT算子本身，思维过程的体现一点没有😭 ，难怪我看了几小时代码愣是啥都没看出来。

计算矩阵（这里设 $0\leqslant{i,j}\leqslant{n-1}$）：

DFT：$A_{i,j}=W_n^{ij}$

IDFT：$A_{i,j}=\cfrac{-W_n^{ij}}n$

以上两矩阵之积正好为一个单位矩阵，

注意运算前一定要把 $n$ 补成 $2^k$ 的形式。

//虚数类
struct cpl
{
    double a,b;
    cpl operator+(cpl y)
    {
        return (cpl){a+y.a,b+y.b};
    }
    cpl operator-(cpl y)
    {
        return (cpl){a-y.a,b-y.b};
    }
    cpl operator*(cpl y)
    {
        return (cpl){a*y.a-b*y.b,b*y.a+a*y.b};
    }
    bool operator==(cpl y)
    {
        return a==y.a&&b==y.b;
    }
};
//FFT函数
int v[maxn];
void fft(cpl cl[],int len,bool ck)
{
    int lt=(1<<len);
    fru(i,0,lt-1)
    {
        v[i]=v[i>>1]>>1;
        if(i&1)
        {
            v[i]|=(lt>>1);
        }
        if(v[i]>i)swap(cl[i],cl[v[i]]);
    }
    cpl w,w1,T1,T2;
    for(register int i=1;i<lt;i<<=1)
    {
        w1=ck?(cpl){cos(-2.0*pi/(i<<1)),sin(-2.0*pi/(i<<1))}:(cpl){cos(2.0*pi/(i<<1)),sin(2.0*pi/(i<<1))};
        for(register int j=0;j<lt;j+=(i<<1))
        {
            w=(cpl){1,0};
            for(register int k=j;k<i+j;k++)
            {
                // cout<<k<<' '<<k+i<<endl;
                T1=cl[k]+(w*cl[k+i]);
                T2=cl[k]-(w*cl[k+i]);
                cl[k]=T1;cl[k+i]=T2;
                w=w*w1;
            }
        }
    }
    if(ck)fru(i,0,lt-1)cl[i].a/=lt;
}