数论模板

Gavin
  2026-03-11 08:41 2026-03-11 08:41 创建   2026-04-15 11:49 2026-04-15 11:49 更新   6.4k 字

exgcd 求解逆元

1
2
3
// ax%b==1
exgcd(a,b,x,y);
ans = (x+b)%b;

Miller-Rabin 素性测试 & Pollard-Rho

Code
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#include <random>

const int TEST_TIME=16;
random_device rd;
mt19937_64 gen(rd());

ll bmul(ll a,ll b,ll m){
  ull c=(ull)a*(ull)b-(ull)((ld)a/m*b+0.5L)*(ull)m;
  if(c<(ull)m){
      return c;
  }
  return c+m;
}

ll qpow(ll a,ll b,ll P){
    ll res=1;
    while(b>0){
        if(b&1){
            res = bmul(res,a,P);
        }
        a = bmul(a,a,P);
        b >>= 1;
    }
    return res;
}

bool MillerRabin(ll n){
    if(n<2){
        return false;
    }
    if(n==2||n==3){
        return true;
    }

    ll u=n-1,t=0;
    while(u%2==0){
        u /= 2;
        t++;
    }

    uniform_int_distribution<ll> dist(0,n-4);

    for(int i=1;i<=TEST_TIME;i++){
        ll a=dist(gen)+2,v=qpow(a,u,n),s=0;
        if(v==1){
            continue;
        }
        for(s=0;s<t;s++){
            if(v==n-1){
                break;
            }
            v = bmul(v,v,n);
        }
        if(s==t){
            return false;
        }
    }

    return true;
}

ll gcd(ll a,ll b){
    return b==0?a:gcd(b,a%b);
}

ll PollardRho(ll x){
    ll s=0,t=0,c=uniform_int_distribution<>(0,x-2)(gen)+1,val=1;

    for(int gl=1;;gl*=2,s=t,val=1){
        for(int stp=1;stp<=gl;stp++){
            t = (bmul(t,t,x)+c)%x;
            val = bmul(val,abs(t-s),x);
            if((stp%127)==0){
                ll d=gcd(val,x);
                if(d>1){
                    return d;
                }
            }
        }
        ll d=gcd(val,x);
        if(d>1){
            return d;
        }
    }
}

ll mxfac;
void maxFac(ll x,bool fst=true){
    if(fst){
        mxfac = 0;
    }

    if(x<=mxfac||x<2){
        return;
    }

    if(MillerRabin(x)){
        mxfac = max(mxfac,x);
        return;
    }

    ll p=x;
    while(p>=x){
        p = PollardRho(x);
    }
    while((x%p)==0){
        x /= p;
    }

    maxFac(x,false);
    maxFac(p,false);
}

#include <vector>

vector<ll> facs;
void getFacs(ll x,bool fst=true){
    if(fst){
        facs.clear();
    }

    if(x==1){
        return;
    }

    if(MillerRabin(x)){
        facs.push_back(x);
        return;
    }

    ll p=x;
    while(p==x){
        p = PollardRho(x);
    }

    getFacs(p,false);
    getFacs(x/p,false);
}

CRT & exCRT

CRT
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){
        x = 1;
        y = 0;
    }
    else{
        exgcd(b,a%b,y,x);
        y -= a/b*x;
    }
}

ll crt(int k,ll* a,ll* r){
    ll n=1,ans=0;
    for(int i=1;i<=k;i++){
        n = n*a[i];
    }
    for(int i=1;i<=k;i++){
        ll m=n/a[i],b,y;
        exgcd(m,a[i],b,y);
        ans = (ans+r[i]*m*b%n)%n;
    }
    return (ans%n+n)%n;
}
exCRT
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
/******************************
 - @GavinCQTD / 2026-03-13 15:12:07
 - "P4777 【模板】扩展中国剩余定理（EXCRT）" From Luogu
 - # https://www.luogu.com.cn/problem/P4777
 - 1000 ms / 500 MB
******************************/

#include <iostream>
#include <cassert>
using namespace std;
using ll = __int128_t;
using ull = unsigned long long;
using ld = long double;

ll n,a[100005],b[100005],A=1,B;

ll exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){
        x = 1;
        y = 0;
        return a;
    }
    else{
        ll res=exgcd(b,a%b,y,x);
        y -= a/b*x;
        return res;
    }
}

// FastIO...

int main(){
    read(n);
    for(int i=1;i<=n;i++){
        read(a[i],b[i]);
        b[i] %= a[i];

        ll x,y,gres=exgcd(A,a[i],x,y);
        // (B-r[i])%gres!=0: No solution
        x = (B-b[i])/gres*x;
        B = B-A*x;
        A = a[i]/gres*A;
        B = (B%A+A)%A;
    }

    put((B%A+A)%A);

    return 0;
}

Lagrange Interpolation

Code
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
#include <vector>

ll qpow(ll a,ll b,ll P){
    ll res=1;
    while(b>0){
        if(b&1){
            res = res*a%P;
        }
        a = a*a%P;
        b >>= 1;
    }
    return res;
}

ll inv(ll a,ll P){
    return qpow(a,P-2,P);
}

vector<int> lagrangeInterpolation(const vector<int> &x,
                                  const vector<int> &y,
                                  ull P){
    int n=x.size();
    vector<int> M(n+1),px(n,1),f(n);
    M[0] = 1;

    for(int i=0;i<n;i++){
        for(int j=i;j>=0;j--){
            M[j+1] = (M[j]+M[j+1])%P;
            M[j] = (ll)M[j]*(P-x[i])%P;
        }
    }

    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i!=j){
                px[i] = (ll)px[i]*(x[i]-x[j]+P)%P;
            }
        }
    }

    for(int i=0;i<n;i++){
        ll t=(ll)y[i]*inv(px[i],P)%P,k=M[n];
        for(int j=n-1;j>=0;j--){
            f[j] = (f[j]+k*t)%P;
            k = (M[j]+k*x[i])%P;
        }
    }

    return f;
}

Lucas & exLucas

Lucas实现1
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
ll qpow(ll a,ll b,ll P){
    ll res=1;
    while(b>0){
        if(b&1){
            res = res*a%P;
        }
        a = a*a%P;
        b >>= 1;
    }
    return res;
}

ll inv(ll a,ll P){
    return qpow(a,P-2,P);
}

ll fac[100005],ifac[100005];

void init(int P){
    fac[0] = 1;
    for(int i=1;i<P;i++){
        fac[i] = fac[i-1]*i%P;
    }
    ifac[P-1] = inv(fac[P-1],P);
    for(int i=P-1;i>=1;i--){
        ifac[i-1] = ifac[i]*i%P;
    }
}

ll C(int n,int m,int P){
    if(n<m){
        return 0;
    }
    if(n<P&&m<P){
        return fac[n]*ifac[n-m]%P*ifac[m]%P;
    }
    return C(n/P,m/P,P)*C(n%P,m%P,P)%P;
}
Lucas实现2
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
class BinomModPrime{
    int p;
    vector<int> fac,ifac;

    int calc(int n,int k){
        if(n<k){
            return 0;
        }
        ll res=fac[n];
        res = (res*ifac[k])%p;
        res = (res*ifac[n-k])%p;
        return res;
    }

public:
    BinomModPrime(int p):p(p),fac(p),ifac(p){
        fac[0] = 1;
        for(int i=1;i<p;i++){
            fac[i] = (ll)fac[i-1]*i%p;
        }
        ifac[p-1] = p-1;
        for(int i=p-1;i>=1;i--){
            ifac[i-1] = (ll)ifac[i]*i%p;
        }
    }

    int binomial(ll n,ll k){
        ll res=1;
        while(n>0||k>0){
            res = (res*calc(n%p,k%p))%p;
            n /= p;
            k /= p;
        }
        return res;
    }

    int operator()(const ll &n,const ll &k){
        return binomial(n,k);
    }
};
exLucas
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
#include <vector>

void exgcd(int a,int b,int &x,int &y){
    if(b==0){
        x = 1;
        y = 0;
    }
    else{
        exgcd(b,a%b,y,x);
        y -= a/b*x;
    }
}

int inv(int a,int P){
    int x,y;
    exgcd(a,P,x,y);
    return (x%P+P)%P;
}

int crtCoeff(int mi,int m){
    ll mm=m/mi;
    mm *= inv(mm,mi);
    return mm%m;
}

class BinomModPrimePower{
    int p,a,pa;
    vector<int> f;

    ll nu(ll n){
        ll cnt=0;
        do{
            n /= p;
            cnt += n;
        }while(n>0);
        return cnt;
    }

    ll factMod(ll n){
        bool neg=(p!=2)||(pa<=4);
        ll res=1;
        while(n>1){
            if((n/pa)&neg){
                res = pa-res;
            }
            res = res*f[n%pa]%pa;
            n /= p;
        }
        return res;
    }

public:
    BinomModPrimePower(int p,int a,int pa):p(p),a(a),pa(pa),f(pa){
        f[0] = 1;
        for(int i=1;i<pa;i++){
            f[i] = i%p!=0?(ll)f[i-1]*i%pa:f[i-1];
        }
    }

    int binomial(ll n,ll k){
        ll v=nu(n)-nu(n-k)-nu(k);
        if(v>=a){
            return 0;
        }
        auto res=factMod(n-k)*factMod(k)%pa;
        res = factMod(n)*inv(res,pa)%pa;
        for(;v!=0;v--){
            res *= p;
        }
        return res%pa;
    }
};

class BinomMod{
    int m;
    vector<BinomModPrimePower> bp;
    vector<ll> crtM;

public:
    BinomMod(int n):m(n){
        for(int p=2;p*p<=n;p++){
            if(n%p==0){
                int a=0,pa=1;
                for(;n%p==0;n/=p,a++,pa*=p);
                bp.emplace_back(p,a,pa);
                crtM.emplace_back(crtCoeff(pa,m));
            }
        }
        if(n>1){
            bp.emplace_back(n,1,n);
            crtM.emplace_back(crtCoeff(n,m));
        }
    }

    int binomial(ll n,ll k){
        ll res=0;
        for(size_t i=0;i!=bp.size();i++){
            res = (bp[i].binomial(n,k)*crtM[i]+res)%m;
        }
        return res;
    }

    int operator()(const ll &n,const ll &k){
        return binomial(n,k);
    }
};

原根

求解原根
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
/******************************
 - @GavinCQTD / 2026-03-23 09:25:56
 - "P6091 【模板】原根" From Luogu
 - # https://www.luogu.com.cn/problem/P6091
 - 3000 ms / 250 MB
******************************/

#include <iostream>
#include <cassert>
using namespace std;
using ll = long long;
using ull = unsigned long long;
using ld = long double;

#include <bitset>
#include <vector>
#include <random>

random_device rd;
mt19937_64 gen(rd());

ll randint(ll l,ll r){
    return uniform_int_distribution<ll>(l,r)(gen);
}

vector<ll> prms;

ll qpow(ll a,ll b,ll P){
    ll res=1;
    while(b>0){
        if(b&1){
            res = res*a%P;
        }
        a = a*a%P;
        b >>= 1;
    }
    return res;
}

void pri(ll x){
    prms.clear();
    for(ll i=2;i*i<=x;i++){
        if(x%i==0){
            prms.push_back(i);
            while(x%i==0){
                x /= i;
            }
        }
    }
    if(x>1){
        prms.push_back(x);
    }
}

ll phi(ll x){
    pri(x);
    for(ll i:prms){
        x = x/i*(i-1);
    }
    return x;
}

bitset<1000005> u,v;

int T,n,d;

void solve(){
    u.set();
    v.reset();

    cin >> n >> d;

    int m=phi(n),g=0;
    pri(m);

    for(int i=1,r;i<=200;i++){
        r = randint(1,n-1);
        if(qpow(r,m,n)==1){
            bool flg=true;
            for(ll j:prms){
                if(qpow(r,m/j,n)==1){
                    flg = false;
                    break;
                }
            }
            if(flg){
                g = r;
                break;
            }
        }
    }

    if(g==0){
        cout << "0\n\n";
        return;
    }

    for(ll i:prms){
        for(ll j=1;j<i&&i*j<=n;j++){
            u[i*j] = false;
        }
    }

    for(ll i=1;i<=n;i++){
        for(ll j:prms){
            if(j>i||i*j>n){
                break;
            }
            u[i*j] = false;
            if(i%j==0){
                break;
            }
        }
    }

    for(ll i=1,p=g;i<=n;i++,p=p*g%n){
        if(u[i]){
            v[p] = 1;
        }
    }

    vector<int> ans;
    for(ll i=1;i<=n;i++){
        if(v[i]){
            ans.push_back(i);
        }
    }

    cout << ans.size() << "\n";
    for(ll i=1;i<=(int)ans.size()/d;i++){
        cout << ans[i*d-1] << " ";
    }
    cout << "\n";
}

int main(){
    cin >> T;

    for(int i=1;i<=T;i++){
        solve();
    }

    return 0;
}
原根判定
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
/******************************
 - @GavinCQTD / 2026-03-23 11:20:46
 - "Primitive Root" From SPOJ - Classical
 - # https://www.spoj.com/problems/PROOT/en/
 - 1000 ms / 1536 MB
******************************/

#include <iostream>
#include <cassert>
using namespace std;
using ll = long long;
using ull = unsigned long long;
using ld = long double;

#include <bitset>
#include <vector>

bitset<100005> ntp;
vector<int> prms;

void init(){
    for(int i=2;i<=1e5;i++){
        if(!ntp[i]){
            prms.push_back(i);
        }
        for(int pj:prms){
            if(i*pj>1e5){
                break;
            }
            ntp[i*pj] = true;
            if(i%pj==0){
                break;
            }
        }
    }
}

ll qpow(ll a,ll b,ll P){
    ll res=1;
    while(b>0){
        if(b&1){
            res = res*a%P;
        }
        a = a*a%P;
        b >>= 1;
    }
    return res;
}

vector<int> facs;
void getf(int x){
    facs.clear();
    for(int prm:prms){
        if(x%prm==0){
            facs.push_back(prm);
        }
        while(x%prm==0){
            x /= prm;
        }
    }
}

bool check(int p,int x){
    for(int fac:facs){
        int u=(p-1)/fac,tmp=qpow(x,u,p);
        if(tmp==1){
            return false;
        }
    }
    return true;
}

int p,n;

int main(){
    init();

    while(cin>>p>>n){
        if(p==0&&n==0){
            break;
        }

        getf(p-1);
        for(int i=1;i<=n;i++){
            int x;
            cin >> x;
            cout << (check(p,x)?"YES\n":"NO\n");
        }
    }

    return 0;
}

BSGS & exBSGS

Code
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#include <unordered_map>
#include <cmath>

unordered_map<ll,ll> mp;

ll qpow(ll a,ll b,ll P){
    ll res=1;
    while(b>0){
        if(b&1){
            res = res*a%P;
        }
        a = a*a%P;
        b >>= 1;
    }
    return res;
}

ll exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){
        x = 1;
        y = 0;
        return a;
    }

    ll res=exgcd(b,a%b,y,x);
    y -= a/b*x;
    return res;
}

ll inv(ll a,ll P){
    ll x,y;
    exgcd(a,P,x,y);
    return (x%P+P)%P;
}

ll BSGS(ll a,ll b,ll P){
    if(a%P==b%P){
        return 1;
    }
    if(a%P==0&&b!=0){
        return -1;
    }

    mp.clear();

    ll uni=(ll)ceil(sqrt(P)),tmp=qpow(a,uni,P);
    mp.reserve(4*uni);

    for(int i=0;i<=uni;i++){
        mp[b] = i;
        b = b*a%P;
    }

    b = 1;

    for(int i=1;i<=uni;i++){
        b = b*tmp%P;
        if(mp.count(b)){
            return i*uni-mp[b];
        }
    }

    return -1;
}

ll exBSGS(ll a,ll b,ll P){
    a %= P;
    b %= P;

    if(b==1||P==1){
        return 0;
    }

    ll x,y,g=exgcd(a,P,x,y),k=0,tmp=1;

    while(g!=1){
        if(b%g!=0){
            return -1;
        }
        k++;
        b /= g;
        P /= g;
        tmp = tmp*(a/g)%P;
        if(tmp==b){
            return k;
        }
        g = exgcd(a,P,x,y);
    }

    ll ans=BSGS(a,b*inv(tmp,P)%P,P);
    if(ans==-1){
        return -1;
    }
    return ans+k;
}

康托展开&逆展开

Code
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
/***************
 - @GavinCQTD
***************/

#include <iostream>
#include <cassert>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
using ll = long long;
using ull = unsigned long long;
using ld = long double;

class Fenwick{
private:
    int n;
    vector<int> tree;

    int lowbit(int x){
        return x&(-x);
    }

public:
    int query(int x){
        int ans=0;
        while(x!=0){
            ans += tree[x];
            x -= lowbit(x);
        }
        return ans;
    }

    void update(int x,int k){
        while(x<=n){
            tree[x] += k;
            x += lowbit(x);
        }
    }

    void fill(){
        for(int i=1;i<=n;i++){
            tree[i] += lowbit(i);
        }
    }

    int kth(int k){
        int ps=0,x=0;
        for(int i=(int)log2(n);i>=0;i--){
            int nxtx=x+(1<<i);
            if(nxtx<=n&&ps+tree[nxtx]<k){
                x = nxtx;
                ps += tree[x];
            }
        }
        return x+1;
    }

    Fenwick(int n_):n(n_),tree(n_+1,0){}
};

const int MOD=998244353;

ll cantor(int n,const vector<int>& a){
    Fenwick fen(n+1);
    ll fac=1,res=0;
    for(int i=n-1;i>=0;i--){
        res = (res+(ll)fen.query(a[i]-1)*fac)%MOD;
        fen.update(a[i],1);
        fac = (fac*(n-i))%MOD;
    }

    return (res+1)%MOD;
}

vector<int> revCantor(int n,ll k){
    k--;
    vector<int> lem(n);
    for(int i=1;i<=n;i++){
        lem[n-i] = k % i;
        k /= i;
    }
    Fenwick fen(n+1);
    fen.fill();
    vector<int> res(n);
    for(int i=0;i<n;i++){
        res[i] = fen.kth(lem[i]+1);
        fen.update(res[i],-1);
    }
    return res;
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int T;
    if(!(cin >> T)) return 0;

    for(int i=1;i<=T;i++){
        int op, n;
        cin >> op;

        if(op==0){
            cin >> n;
            vector<int> a(n);
            for(int j=0;j<n;j++){
                cin >> a[j];
            }
            cout << solve1(n,a) << "\n";
        }

        else{
            ll rk;
            cin >> n >> rk;
            auto res=solve2(n,rk);
            for(int j=0;j<n;j++){
                cout << res[j] << " ";
            }
            cout << "\n";
        }
    }

    return 0;
}

还有一个“可重元素集的康托展开”，例题是 [HAOI2010] 计数。

高次剩余

二次剩余
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
/******************************
 - @GavinCQTD / 2026-04-14 11:52:45
 - "P5491 【模板】二次剩余" From Luogu
 - # https://www.luogu.com.cn/problem/P5491
 - 1000 ms / 125 MB
******************************/

#include <iostream>
#include <cassert>
using namespace std;
using ll = long long;
using ull = unsigned long long;
using ld = long double;

#include <utility>
#include <random>

ll p,v;

struct Poly{
    ll a,b;
    Poly(ll a=0,ll b=0):a(a),b(b){}

    friend Poly operator*(const Poly &x,const Poly &y){
        return Poly((x.a*y.a+v*(x.b*y.b%p))%p,(x.a*y.b+x.b*y.a)%p);
    }
};

Poly qpow(Poly a,ll b){
    Poly res(1,0);
    while(b>0){
        if(b&1){
            res = res*a;
        }
        a = a*a;
        b >>= 1;
    }
    return res;
}

ll qpow(ll a,ll b){
    a %= p;
    ll res=1;
    while(b>0){
        if(b&1){
            res = res*a%p;
        }
        a = a*a%p;
        b >>= 1;
    }
    return res;
}

mt19937_64 gen(random_device{}());

ll cipolla(ll a,ll _p){
    p = _p;
    if(a==0){
        return 0;
    }
    else if(qpow(a,(p-1)/2)==p-1){
        return -1;
    }
    else{
        ll r;
        for(r=gen()%p;;r=gen()%p){
            if(qpow(((r*r%p-a)%p+p)%p,(p-1)/2)==p-1){
                break;
            }
        }
        v = ((r*r%p-a)%p+p)%p;
        return qpow(Poly(r,1),(p+1)/2).a;
    }
}

pair<ll,ll> getAns(ll a,ll _p){
    ll cans=cipolla(a,_p);
    if(cans<=0){
        return {cans,-1};
    }
    else{
        ll ans2=(_p-cans)%_p;
        if(ans2<cans){
            swap(cans,ans2);
        }
        return {cans,ans2};
    }
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

    int T;
    cin >> T;

    while(T--){
        int n,p;
        cin >> n >> p;

        auto [ans1,ans2]=getAns(n,p);
        if(ans1==-1){
            cout << "Hola!\n";
        }
        else if(ans1==0){
            cout << "0\n";
        }
        else{
            cout << ans1 << " " << ans2 << "\n";
        }
    }

    return 0;
}
任意模数二次剩余
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
/***************
 - @GavinCQTD
***************/

#include <iostream>
#include <cassert>
using namespace std;
using ll = long long;
using ull = unsigned long long;
using ld = long double;

#include <algorithm>
#include <bitset>
#include <vector>
#include <random>

using i128 = __int128_t;
const int LINEAR_LIMIT = 300000;

bitset<LINEAR_LIMIT+5> ntp;
vector<int> prms;

void init(){
    for(int i=2;i<=LINEAR_LIMIT;i++){
        if(!ntp[i]){
            prms.push_back(i);
        }
        for(int pj:prms){
            if(i*pj>LINEAR_LIMIT){
                break;
            }
            ntp[i*pj] = true;
            if(i%pj==0){
                break;
            }
        }
    }
}

ll bmul(ll a,ll b,ll P){
    ull c=(ull)a*(ull)b-(ull)((ld)a/P*b+0.5L)*(ull)P;
    if(c<(ull)P){
        return c;
    }
    return c+P;
}

ll qpow(ll a,ll b,ll P){
    a %= P;
    ll res=1;
    while(b>0){
        if(b&1){
            res = bmul(res,a,P);
        }
        a = bmul(a,a,P);
        b >>= 1;
    }
    return res;
}

ll exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){
        x = 1;
        y = 0;
        return a;
    }
    ll res=exgcd(b,a%b,y,x);
    y -= a/b*x;
    return res;
}

ll inv(ll a,ll P){
    ll x,y;
    exgcd(a,P,x,y);
    return (x%P+P)%P;
}

ll v,MOD;

struct Poly{
    ll a,b;
    Poly(ll a=0,ll b=0):a(a),b(b){}

    friend Poly operator*(const Poly &x,const Poly &y){
        return Poly((bmul(x.a,y.a,MOD)+bmul(v,bmul(x.b,y.b,MOD),MOD))%MOD
                    ,(bmul(x.a,y.b,MOD)+bmul(x.b,y.a,MOD))%MOD);
    }
};

Poly qpow(Poly a,ll b){
    Poly res(1,0);
    while(b>0){
        if(b&1){
            res = res*a;
        }
        a = a*a;
        b >>= 1;
    }
    return res;
}

const int TEST_TIME=16;
mt19937_64 gen(random_device{}());

bool MillerRabin(ll n){
    if(n<2){
        return false;
    }
    if(n==2||n==3){
        return true;
    }

    ll u=n-1,t=0;
    while(u%2==0){
        u /= 2;
        t++;
    }

    uniform_int_distribution<ll> dist(0,n-4);

    for(int i=1;i<=TEST_TIME;i++){
        ll a=dist(gen)+2,v=qpow(a,u,n),s=0;
        if(v==1){
            continue;
        }
        for(s=0;s<t;s++){
            if(v==n-1){
                break;
            }
            v = bmul(v,v,n);
        }
        if(s==t){
            return false;
        }
    }

    return true;
}

ll PollardRho(ll x){
    ll s=0,t=0,c=uniform_int_distribution<>(0,x-2)(gen)+1,val=1;

    for(int gl=1;;gl*=2,s=t,val=1){
        for(int stp=1;stp<=gl;stp++){
            t = (bmul(t,t,x)+c)%x;
            val = bmul(val,abs(t-s),x);
            if((stp%127)==0){
                ll d=gcd(val,x);
                if(d>1){
                    return d;
                }
            }
        }
        ll d=gcd(val,x);
        if(d>1){
            return d;
        }
    }
}

vector<ll> facs;
void getFacs(ll x,bool fst=true){
    if(fst){
        facs.clear();
    }

    if(x==1){
        return;
    }

    if(MillerRabin(x)){
        facs.push_back(x);
        return;
    }

    ll p=x;
    while(p==x){
        p = PollardRho(x);
    }

    getFacs(p,false);
    getFacs(x/p,false);
}

ll solve2k(ll a,ll k){
    if(a%(1ll<<k)==0){
        return 0;
    }
    a %= (1ll<<k);

    ll t=0,res=1;
    while(a%2==0){
        a /= 2;
        t++;
    }

    if((a&7)^1){
        return -1;
    }
    if(t&1){
        return -1;
    }
    k -= t;
    for(int i=4;i<=k;i++){
        res = (res+(a%(1ll<<i)-res*res)/2)%(1ll<<k);
    }
    res %= (1ll<<k);
    if(res<0){
        res += (1ll<<k);
    }
    return res<<(t>>1);
}

ll solveP(ll a,ll P){
    a %= P;
    if(qpow(a,(P-1)/2,P)==P-1){
        return -1;
    }
    ll b;
    MOD = P;
    uniform_int_distribution<ll> dist(0,P-1);
    while(true){
        b = dist(gen);
        v = (bmul(b,b,P)-a+P)%P;
        if(qpow(v,(P-1)/2,P)==P-1){
            break;
        }
    }
    ll ans=qpow(Poly(b,1),(P+1)/2).a;
    return min(ans,P-ans);
}

ll solvePK(ll a,ll k,ll p,ll mod){
    if(a%mod==0){
        return 0;
    }
    ll t=0,hmod=1;
    while(a%p==0){
        a /= p;
        t++;
        hmod *= (t&1)?p:1;
    }
    if(t&1){
        return -1;
    }
    k -= t;
    mod /= hmod*hmod;
    ll res=solveP(a,p);
    if(res==-1){
        return -1;
    }
    Poly tmp(res,1);
    v = a;
    MOD = mod;
    tmp = qpow(tmp,k);
    res = bmul(tmp.a,inv(tmp.b,MOD),MOD);
    return res*hmod;
}

ll rm[105],md[105];
ll crt(ll P){
    ll ans=0;
    for(int i=1;i<=(int)facs.size();i++){
        ans = (ans+bmul(bmul(P/md[i],inv(P/md[i],md[i]),P),rm[i],P))%P;
    }
    return ans;
}

ll solve(ll a,ll pmod){
    a %= pmod;
    ll crtp=pmod;
    getFacs(pmod);
    sort(facs.begin(),facs.end());
    facs.erase(unique(facs.begin(),facs.end()),facs.end());
    for(int i=1;i<=(int)facs.size();i++){
        ll nw=0,rmod=1;
        while(pmod%facs[i-1]==0){
            pmod /= facs[i-1];
            nw++;
            rmod *= facs[i-1];
        }
        md[i] = rmod;
        if(facs[i-1]==2){
            rm[i] = solve2k(a,nw);
        }
        else{
            rm[i] = solvePK(a,nw,facs[i-1],rmod);
        }
        if(rm[i]==-1){
            return -1;
        }
    }
    return crt(crtp);
}

int main(){
    init();

    int T;
    cin >> T;

    while(T--){
        ll a,p;
        cin >> a >> p;
        cout << solve(a,p) << "\n";
    }

    return 0;
}
N次剩余
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
/******************************
 - @GavinCQTD / 2026-04-15 10:18:45
 - "P5668 【模板】N 次剩余" From Luogu
 - # https://www.luogu.com.cn/problem/P5668
 - 2000 ms / 128 MB
******************************/

#include <algorithm>
#include <iostream>
#include <cassert>
using namespace std;
using ll = long long;
using ull = unsigned long long;
using ld = long double;

#include <unordered_map>
#include <numeric>
#include <random>
#include <vector>
#include <cmath>

struct PrimePower{
    int p,e,pe;
    PrimePower(int p,int e,int pe):p(p),e(e),pe(pe){}
};

auto factorize(int n){
    vector<PrimePower> res;
    for(int x=2;x*x<=n;x++){
        int e=0,pe=1;
        for(;n%x==0;n/=x,e++,pe*=x);
        if(e!=0){
            res.emplace_back(x,e,pe);
        }
    }
    if(n>1){
        res.emplace_back(n,1,n);
    }
    return res;
}

int qpow(int a,int b,int P=0){
    int res=1;
    while(b>0){
        if(b&1){
            res = (P==0?res*a:(ll)res*a%P);
        }
        a = (P==0?a*a:(ll)a*a%P);
        b >>= 1;
    }
    return res;
}

int primitiveRoot(PrimePower pp){
    vector<int> exp;
    int phi=pp.pe/pp.p*(pp.p-1);
    for(auto fac:factorize(pp.p-1)){
        exp.push_back(phi/fac.p);
    }
    if(pp.e!=1){
        exp.push_back(phi/pp.p);
    }
    int ans=0;
    bool ok=false;
    while(!ok){
        ans++;
        ok = true;
        for(int b:exp){
            if(qpow(ans,b,pp.pe)==1){
                ok = false;
                break;
            }
        }
    }
    return ans;
}

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

int inv(int a,int P){
    int x,y;
    exgcd(a,P,x,y);
    return (x%P+P)%P;
}

mt19937_64 gen(random_device{}());

int nonResidue(int p,int m,int phi){
    uniform_int_distribution<> dist(1,m-1);
    while(true){
        int c=dist(gen);
        if(gcd(c,m)==1&&qpow(c,phi/p,m)!=1){
            return c;
        }
    }
    return -1;
}

int pethRootModM(int p,int e,int a,int m,int phi){
    if(m==2){
        return 1;
    }
    int s=0,r=phi,pe=qpow(p,e);
    for(;r%p==0;r/=p,s++);
    int q=pe-inv(r,pe),ans=qpow(a,((ll)q*r+1)/pe%phi,m)
        ,eta=nonResidue(p,m,phi);
    unordered_map<int,int> mp;
    int zeta=qpow(eta,r,m),xi=qpow(eta,phi/p,m),B=sqrt((s-e)*p+0.25l)+1
        ,pB=pe/B+1,po0=qpow(xi,pB,m);
    for(int j=1,po1=1;j<=B;j++){
        po1 = (ll)po1*po0%m;
        mp[po1] = j;
    }
    for(int j=0;j<s-e;j++){
        int err=(ll)qpow(ans,pe,m)*inv(a,m)%m,xihj=qpow(err,qpow(p,s-e-j-1),m);
        ll hj=0;
        for(int i=1;i<=pB;i++){
            xihj = (ll)xihj*xi%m;
            if(mp.count(xihj)!=0){
                hj = mp[xihj]*pB-i;
                break;
            }
        }
        ans = (ll)ans*qpow(zeta,phi-hj*qpow(p,j)%phi,m)%m;
    }
    return ans;
}

int kthRootModPe(int k,int a,int pe,int phi){
    a %= pe;
    if(k==0){
        return a==1?0:-1;
    }
    if(a==0){
        return 0;
    }
    int d=gcd(k,phi);
    if(qpow(a,phi/d,pe)!=1){
        return -1;
    }
    a = qpow(a,inv(k/d,phi/d),pe);
    for(int dp=2;dp*dp<=d&&dp*dp*dp*dp<pe;dp++){
        if(d%dp==0){
            int de=0;
            for(;d%dp==0;d/=dp,de++);
            a = pethRootModM(dp,de,a,pe,phi);
        }
    }
    if(d!=1){
        int dp=gcd(d,phi/d),de=0;
        if(dp!=1){
            for(;d%dp==0;d/=dp,de++);
            a = pethRootModM(dp,de,a,pe,phi);
        }
        if(d!=1){
            a = pethRootModM(d,1,a,pe,phi);
        }
    }
    return a;
}

auto calc(int g,int k,int a,int p,int pe){
    int mm=(p==2?pe/4:pe/p*(p-1)),ans=kthRootModPe(k,a,pe,mm);
    if(ans==-1){
        return vector<int>();
    }
    int d=mm/gcd(k,mm),po=qpow(g,d,pe);
    vector<int> res(mm/d);
    for(auto &x:res){
        x = ans;
        ans = (ll)ans*po%pe;
    }
    return res;
}

auto kthRootsModPe(int k,int a,PrimePower pp){
    int p=pp.p,e=pp.e,pe=pp.pe;
    a %= pe;
    if(a==0){
        int d=qpow(p,(e-1)/k+1);
        vector<int> res(pe/d);
        for(int i=0;i<pe/d;i++){
            res[i] = i*d;
        }
        return res;
    }

    int s=0;
    for(;a%p==0;a/=p,s++);
    if(s%k){
        return vector<int>();
    }
    int psk=qpow(p,s/k),pss=qpow(p,s-s/k),pes=qpow(p,e-s);
    vector<int> res;
    if(p!=2){
        int g=primitiveRoot(PrimePower(p,e-s,pes));
        res = calc(g,k,a,p,pes);
    }
    else if(pes==2){
        res.push_back(a);
    }
    else if(k&1){
        bool jd=(a%4==3);
        a = (jd?pes-a:a);
        res = calc(5,k,a,p,pes);
        if(jd){
            for(auto &x:res){
                x = pes-x;
            }
        }
    }
    else{
        if(a%4==3){
            return vector<int>();
        }
        res = calc(5,k,a,p,pes);
        int m=res.size();
        res.reserve(m*2);
        for(int i=0;i<m;i++){
            res.push_back(pes-res[i]);
        }
    }
    int m=res.size();
    res.reserve(m*pss);
    for(int j=1;j<pss;j++){
        for(int i=0;i<m;i++){
            res.push_back(res.end()[-m]+pes);
        }
    }
    for(auto &x:res){
        x *= psk;
    }
    return res;
}

auto kthRootsModM(int k,int a,int m){
    auto facs=factorize(m);
    int m0=0;
    vector<vector<int>> sols;
    for(auto pp:facs){
        sols.push_back(kthRootsModPe(k,a,pp));
        if(sols.back().empty()){
            return vector<int>();
        }
    }
    vector<int> ans;
    for(int i=0;i<(int)facs.size();i++){
        auto pp=facs[i];
        if(i==0){
            m0 = pp.pe;
            ans = sols[i];
        }
        else{
            ll m1=pp.pe*inv(pp.pe,m0),m2=m0*inv(m0,pp.pe);
            m0 *= pp.pe;
            vector<int> _ans;
            _ans.reserve(ans.size()*sols[i].size());
            for(auto x:ans){
                for(auto y:sols[i]){
                    _ans.push_back((m1*x+m2*y)%m0);
                }
            }
            ans.swap(_ans);
        }
    }
    return ans;
}

int main(){
    int T;
    cin >> T;

    while(T--){
        int n,m,k;
        cin >> n >> m >> k;
        auto ans=kthRootsModM(n,k,m);
        cout << ans.size() << "\n";
        if(ans.empty()){
            continue;
        }
        sort(ans.begin(),ans.end());
        for(auto itm:ans){
            cout << itm << " ";
        }
        cout << "\n";
    }

    return 0;
}

基于值域预处理的快速离散对数

跑得很慢，没做什么优化，就是卡着过了

Code
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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
/******************************
 - @GavinCQTD / 2026-04-15 10:35:45
 - "P11175 【模板】基于值域预处理的快速离散对数" From Luogu
 - # https://www.luogu.com.cn/problem/P11175
 - 1500 ms / 512 MB
******************************/

#include <iostream>
#include <cassert>
using namespace std;
using ll = long long;
using ull = unsigned long long;
using ld = long double;

#include <unordered_map>
#include <cmath>

int qpow(int a,int b,int P){
    int res=1;
    while(b>0){
        if(b&1){
            res = (ll)res*a%P;
        }
        a = (ll)a*a%P;
        b >>= 1;
    }
    return res;
}

namespace BSGS{
    unordered_map<int,int> M;
    int B,U,P,g;

    void init(int g,int P0,int B0){
        M.clear();
        B = B0;
        P = P0;
        U = qpow(qpow(g,B,P),P-2,P);
        int w=1;
        for(int i=0;i<B;i++){
            M[w] = i;
            w = (ll)w*g%P;
        }
    }

    int solve(int y){
        int w=y;
        for(int i=0;i<=P/B;i++){
            auto it=M.find(w);
            if(it!=M.end()){
                return i*B+it->second;
            }
            w = (ll)w*U%P;
        }
        return -1;
    }
}

const int MAXN = 1e7+5;

int H[MAXN],P[MAXN],H0,p,h,g,mod;
bool V[MAXN];

int solve(int x){
    if(x<=h){
        return H[x];
    }
    int v=mod/x,r=mod%x;
    if(r<x-r){
        return ((H0+solve(r))%(mod-1)-H[v]+mod-1)%(mod-1);
    }
    else{
        return (solve(x-r)-H[v+1]+mod-1)%(mod-1);
    }
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

    cin >> mod >> g;
    h = sqrt(mod)+1;

    BSGS::init(g,mod,sqrt((ll)mod*sqrt(mod)/log(mod)));
    H0 = BSGS::solve(mod-1)%(mod-1);

    H[1] = 0;
    for(int i=2;i<=h;i++){
        if(!V[i]){
            P[++p] = i;
            H[i] = BSGS::solve(i);
        }
        for(int j=1;j<=p&&P[j]<=h/i;j++){
            int p=P[j];
            H[i*p] = (H[i]+H[p])%(mod-1);
            V[i*p] = true;
            if(i%p==0){
                break;
            }
        }
    }

    int T;
    cin >> T;
    for(int i=1;i<=T;i++){
        int x;
        cin >> x;
        cout << solve(x) << "\n";
    }

    return 0;
}
  • 标题: 数论模板
  • 作者: Gavin
  • 创建于 : 2026-03-11 08:41:00
  • 更新于 : 2026-04-15 11:49:00
  • 链接: https://gavin-blog.pages.dev/2026/数论模板/
  • 版权声明: 本文章采用 CC BY-NC-SA 4.0 进行许可。
目录
数论模板
  1. exgcd 求解逆元
  2. Miller-Rabin 素性测试 & Pollard-Rho
  3. CRT & exCRT
  4. Lagrange Interpolation
  5. Lucas & exLucas
  6. 原根
  7. BSGS & exBSGS
  8. 康托展开&逆展开
  9. 高次剩余
  10. 基于值域预处理的快速离散对数