跳转至

数论

BinaryGCD 快速最大公约数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
namespace BinaryGCD
{
    int gcd(int a,int b)
    {
        int az=__builtin_ctz(a),bz=__builtin_ctz(b);
        int z=min(az,bz),tmp; b>>=bz;
        while(a) a>>=az,tmp=a-b,az=__builtin_ctz(tmp),b=min(a,b),a=abs(tmp);
        return b<<z;
    }
}

线性同余方程

求解 \(ax + by = \gcd(a,b)\)

1
2
3
4
5
6
7
8
9
namespace exEuclidean
{
    int exgcd(int a,int b,int &x,int &y)
    {
        if(!b) return x=1,y=0,a;
        int d=exgcd(b,a%b,y,x);
        return y-=x*(a/b),d;
    }
}

中国剩余定理

求解 \(x \equiv b_i (\mod a_i)\)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
namespace ChineseRemainderTheorem
{
    int N,a[MAX],b[MAX],m[MAX];
    bool flag;
    inline int CRT()
    {
        int M=1,ans=0,x;
        for(int i=1;i<=N;++i) M*=b[i];
        for(int i=1;i<=N;++i) m[i]=M/b[i],exgcd(m[i],b[i],x,*new int),(ans+=a[i]*m[i]*x+M)%=M;
        return (ans+M)%M;
    }
    inline pair<bool,int> EXCRT()
    {
        int M=b[1],ans=a[1];
        for(int i=2;i<=N;++i)
        {
            int y,d=exgcd(M,b[i],*new int,y),c=(ans-a[i])/d;
            flag|=((ans-a[i])%d!=0),M=M/__gcd(M,b[i])*b[i];
            ans=(a[i]+(__int128_t)b[i]*c%M*y)%M;
        }
        return mp(flag,(ans+M)%M);
    }
}

离散对数同余方程

求解 \(a^x \equiv b (\mod p)\)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
namespace BabyStepGiantStep
{
    inline int Inv(int x,int p)
    {
        auto exgcd=[&](auto exgcd,int a,int b,int &x,int &y)->int{if(!b) return x=1,y=0,a; int d=exgcd(exgcd,b,a%b,y,x); return y-=x*(a/b),d;};
        return exgcd(exgcd,x,p,x,*new int),(x%p+p)%p;
    }
    inline int BSGS(int a,int b,int p)
    {
        int A=1,B=sqrt(p)+1; a%=p,b%=p;
        unordered_map<int,int> mp;
        for(int i=1;i<=B;++i) mp[(A=A*a%p)*b%p]=i;
        for(int i=1,AA=A;i<=B;++i,AA=AA*A%p) if(mp.find(AA)!=mp.end()) return i*B-mp[AA];
        return -1; 
    }
    inline int exBSGS(int a,int b,int p)
    {
        int d=__gcd(a,p),cnt=0; a%=p,b%=p;
        while(d>1)
        {
            if(b==1) return cnt;
            if(b%d)  return -1;
            p/=d,b=b/d*Inv(a/d,p)%p,d=__gcd(p,a%=p),++cnt;
        }
        int ans=BSGS(a,b,p);
        return ~ans?ans+cnt:-1;
    }
}

二次剩余

求解 \(x^2 \equiv n (\mod p)\)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
namespace QuadraticResidue
{
    int P,I;
    struct Complex
    {
        int x,y;
        Complex(int X=0,int Y=0):x(X),y(Y) {}
        inline friend Complex operator + (Complex a,Complex b) {return Complex((a.x+b.x)%P,(a.y+b.y)%P);}
        inline friend Complex operator * (Complex a,Complex b) {return Complex((a.x*b.x%P+a.y*b.y%P*I%P)%P,(a.x*b.y%P+a.y*b.x%P)%P);}
        inline friend bool operator == (Complex a,Complex b) {return a.x==b.x&&a.y==b.y;}
    };
    inline int qpow(int a,int b) {int res=1; while(b) {if(b&1) (res*=a)%=P; (a*=a)%=P,b>>=1;} return res;}
    inline Complex qpow(Complex a,int b) {Complex res(1,0); while(b) {if(b&1) res=res*a; a=a*a,b>>=1;} return res;}
    inline int randoom(int l,int r)
    {
        static mt19937_64 sd(20070707^time(NULL));
        static uniform_int_distribution<int> range(l,r);
        return range(sd);
    }
    inline pii Cipolla(int n,int p)
    {
        P=p; int x,x1,x2;
        if(!n) return mp(0,0);
        if(qpow(n,(P-1)>>1)!=1) return mp(-1,-1);
        do x=randoom(1,P-1); while(qpow((x*x+P-n)%P,(P-1)>>1)==1);
        I=(x*x+P-n)%P,x1=qpow(Complex(x,1),(P+1)>>1).x,x2=P-x1;
        return x1>x2?mp(x2,x1):mp(x1,x2);
    }
}

Lucas 定理

求解 \(\binom{n}{m} \mod p\)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
namespace LUCAS
{
    int q,qk,tp;
    int ans,M,fac[MAX];
    inline int Inv(int x,int p)
    {
        auto exgcd=[&](auto exgcd,int a,int b,int &x,int &y)->int{if(!b) return x=1,y=0,a; int d=exgcd(exgcd,b,a%b,y,x); return y-=x*(a/b),d;};
        return exgcd(exgcd,x,p,x,*new int),(x%p+p)%p;
    }
    inline int qpow(int a,int b,int p) {int res=1; while(b) {if(b&1) (res*=a)%=p; (a*=a)%=p,b>>=1;} return res;}
    inline int num(int n) {int s=0; while(n) s+=n/=q; return s;}
    inline int calc(int n) {return n?calc(n/q)*(tp&&(n/qk&1)?qk-1:1)%qk*fac[n%qk]%qk:1;}
    inline int solve(int n,int m)
    {
        tp=(qk&1)||(qk<=4);
        fac[0]=1; for(int i=1;i<qk;++i) fac[i]=fac[i-1]*(i%q?i:1)%qk;
        int pw=num(n)-num(m)-num(n-m),res=calc(n)*Inv(calc(m)*calc(n-m)%qk,qk)%qk;
        while(pw--) (res*=q)%=qk;
        return res;
    }
    inline int exLucas(int n,int m,int p)
    {
        ans=0,M=1;
        for(int i=2;i<=p;++i)
        {
            if(i*i>p) i=p;
            if(p%i==0)
            {
                q=i,qk=1; while(p%i==0) p/=i,qk*=i;
                ans+=(solve(n,m)-ans%qk+qk)*Inv(M,qk)%qk*M,ans%=(M*=qk);
            }
        }
        return ans;
    }
}

Pollard_Rho 分解质因数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
namespace PrimeFactorization
{
    int mix,seed,step;
    static const int P[]={2,3,5,7,11,13,17,19,23,29,31,37};
    static const int Pcnt=12;
    inline int qmul(int a,int b,int mod) {int ans=a*b-(int)((long double)a*b/mod+0.5)*mod; return ans<0?ans+mod:ans;}
    inline int qpow(int a,int b,int mod) {int res=1; while(b) {if(b&1) res=qmul(res,a,mod); a=qmul(a,a,mod); b>>=1;} return res;}
    inline int f(int x,int mod) {return (qmul(x,x,mod)+seed)%mod;}
    inline int randoom(int l,int r)
    {
        static mt19937_64 sd(20070707^time(NULL));
        static uniform_int_distribution<int> range(l,r);
        return range(sd);
    }
    inline bool Miller_Robin(int x)
    {
        if(x<=2) return x==2;
        int y_=x-1; while(!(y_&1)) y_>>=1;
        for(int i=0;i<Pcnt;++i)
        {
            if(x==P[i]) return 1;
            int flag=0,y=y_,z=qpow(P[i],y,x);
            if(z==1) flag=1;
            else while(y<x-1)
            {
                if(z==x-1) {flag=1; break;}
                y<<=1,z=qmul(z,z,x);
            }
            if(!flag) return 0;
        }
        return 1;
    }
    inline int floyd(int x)
    {
        seed=randoom(0,x-1);
        int fast,slow,res=1; fast=slow=randoom(0,x-1);
        fast=f(fast,x);
        for(int i=0;slow!=fast;++i)
        {
            res=qmul(res,(fast-slow)%x+x,x);
            if(!res) res=(fast-slow)%x+x;
            if(i%step==0) {int g=gcd(res,x); if(g!=1) return g; res=1;}
            slow=f(slow,x),fast=f(f(fast,x),x);
        }
        return gcd(res,x);
    }
    inline void Pollard_Rho(int x)
    {
        if(x==1) return;
        if(Miller_Robin(x)) return cmax(mix,x);
        int k=1;
        step=((int)log(x))<<1|1;
        while(k==1) k=floyd(x);
        Pollard_Rho(k),Pollard_Rho(x/k);
    }
    inline int solve(int n)
    {
        mix=0; Pollard_Rho(n);
        if(mix==n) return -1; else return mix;
    }
}

找原根

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
namespace PrimordialRoot
{
    constexpr int N=100000;
    int vis[N+10],h[N+10],phi[N+10],P[MAX],Pcnt;
    inline void PreWork()
    {
        h[1]=h[2]=h[4]=1,phi[1]=1;
        for(int i=2;i<=N;++i)
        {
            if(!vis[i]) P[++Pcnt]=i,h[i]=1,phi[i]=i-1;
            for(int j=1;j<=Pcnt&&i*P[j]<=N;++j)
            {
                vis[i*P[j]]=1;
                if(i%P[j]) phi[i*P[j]]=phi[i]*phi[P[j]];
                else      {phi[i*P[j]]=phi[i]*P[j]; break;}
            }
        }
        for(int i=2;i<=Pcnt;++i)
        {
            for(int j=1;j*P[i]<=N;h[j*=P[i]]=1);
            for(int j=2;j*P[i]<=N;h[j*=P[i]]=1); 
        }
    }
    inline vector<int> solve(int n)
    {
        auto qpow=[](int a,int b,int p)->int{int res=1; while(b) {if(b&1) (res*=a)%=p; (a*=a)%=p,b>>=1;} return res;};
        vector<int> ans,dp,tmp;
        if(!h[n]) {ans.eb(0); return ans;}
        int p=phi[n],g=0;
        for(int i=1;i<=Pcnt&&P[i]<=p;++i) if(p%P[i]==0) dp.eb(P[i]);
        for(int i=1,flag=1;i<n;++i,flag=1)
        {
            if(__gcd(i,n)!=1) continue;
            for(auto d:dp) if(qpow(i,p/d,n)==1) {flag=0; break;}
            if(flag) {g=i; break;}
        }
        for(int s=1,Pow=g;((int)tmp.size())<phi[phi[n]];++s,(Pow*=g)%=n) if(__gcd(p,s)==1) tmp.eb(Pow);
        sort(tmp.begin(),tmp.end()),ans.eb(phi[phi[n]]); for(auto it:tmp) ans.eb(it);
        return ans;
    }
}