[TOC]

数论模板

一, 线性筛

所有数组定义在外面, 自动置为全0即可

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void init() {
    cnt = 0; //统计素数个数
    is_prime[1] = 1; //0代表是素数, 1不是素数
    for (int i = 2; i <= maxn; i++) {
        if (!is_prime[i]) {
            prime[++cnt] = i; //i是素数, prime从下标[1]开始有记录
        }
        for (int j = 1; j <= cnt; j++) { // 1 到 cnt !
            if (prime[j] * i >= maxn)
                break;
            is_prime[i * prime[j]] = 1;
            if (i % prime[j] == 0)
                break;
        }
    }
}

1.分解质因数

通过素数筛法打表, 然后只能通过优化暴力, 最快了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
int get_primefactors(int n) {
    int sum = 0;
    for (int i = 1; i <= tot && pri[i] * pri[i] <= n; ++i) {
        if (n % pri[i] == 0) {
            fac[++sum] = pri[i];
            while (n % pri[i] == 0) {
                n /= pri[i];
            }
        }
    }
    if (n > 1) {
        fac[++sum] = n;
    }
    return sum;
}

2.区间线性筛

a, b可以很大, 但是b - a 的结果可以接受

模板 : 输出a 到 b之间的素数个数

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
#define ll long long
const int MAXN = 100005;
bool is[MAXN];
int pri[MAXN / 10];
int cnt;
bool ans[MAXN];

void init() {
    memset(is, true, sizeof(is));
    is[0] = is[1] = false;
    cnt = 0;
    for (int i = 2; i < MAXN; i++) {
        if (is[i]) {
            pri[++cnt] = i;
        }
        for (int j = 1; j < cnt && pri[j] * i < MAXN; j++) {
            is[pri[j] * i] = false;
            if (i % pri[j] == 0) {
                break;
            }
        }
    }
}

int solve(ll a, ll b) {
    memset(ans, true, sizeof(ans));
    for (int i = 1; i <= cnt && 1ll * pri[i] * pri[i] <= b; i++) {
        ll s = (a + pri[i] - 1) / pri[i];
        if (s < 2) {
            s = 2;
        }
        s *= pri[i];
        for (; s <= b; s += pri[i]) {
            ans[s - a] = false;
        }
    }
    int res = 0;
    for (int i = 0; i <= b - a; i++) {
        if (ans[i]) {
            res++;
        }
    }
    if (a == 1) {
        res--;
    }
    return res;
}

int main() {
    init();
    int a, b;
    cin >> a >> b;
    cout << solve(a, b);
    // 2 , 3 , 5,  7
}

3.求各种积性函数

YC究极版!

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

#include <bits/stdc++.h>
#define ll long long
using namespace std;

const int N = 1e7 + 5;
bool is[N]; //筛数组
int pri[N]; //素数数组
int phi[N]; //欧拉函数数组 !
int mu[N]; //莫⽐乌斯函数数组
int e[N]; //n的最⼩素因⼦在唯⼀分解中的幂次
int d[N]; //因数个数 !
int fac[N]; //i的最⼩的幂次⼤于1的质因数的幂
int tot; //素数个数

void init() {
    memset(is, true, sizeof(is));
    tot = 0;
    is[0] = is[1] = 0;
    mu[1] = 1;
    phi[1] = 1;
    fac[1] = 1;
    e[1] = 0;
    d[1] = 1;
    for (int i = 2; i < N; i++) {
        if (is[i]) {
            pri[++tot] = i;
            phi[i] = i - 1;
            mu[i] = -1;
            e[i] = 1;
            d[i] = 2;
            fac[i] = i;
        }
        for (int j = 1; j <= tot && pri[j] * i < N; j++) {
            is[pri[j] * i] = 0;
            if (i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                mu[i * pri[j]] = 0;
                d[pri[j] * i] = d[i] / (e[i] + 1) * (e[i] + 2);
                e[pri[j] * i] = e[i] + 1;
                fac[i * pri[j]] = fac[i] * pri[j];
                break;
            }
            else {
                mu[i * pri[j]] = -mu[i];
                phi[i * pri[j]] = phi[i] * (pri[j] - 1);
                e[pri[j] * i] = 1;
                d[pri[j] * i] = d[i] * d[pri[j]];
                fac[i * pri[j]] = pri[j];
            }
        }
    }
}


int main() {
    init();
    int a, b;
    cin >> a >> b;
    cout << d[a] << ' ' << d[b];
    return 0;
}

二, Min_25筛

求十亿(1e9)内所有质数的和

$O(n ^{3/4})$

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
// 求比n小的所有素数的和
const int N = 1000010;
typedef long long LL;

namespace Min25 {

    int prime[N], id1[N], id2[N], flag[N], ncnt, m;

    LL g[N], sum[N], a[N], T, n;

    inline int ID(LL x) {
        return x <= T ? id1[x] : id2[n / x];
    }

    inline LL calc(LL x) {
        return x * (x + 1) / 2 - 1;
    }

    inline LL f(LL x) {
        return x;
    }

    inline void init() {
        T = sqrt(n + 0.5);
        for (int i = 2; i <= T; i++) {
            if (!flag[i]) prime[++ncnt] = i, sum[ncnt] = sum[ncnt - 1] + i;
            for (int j = 1; j <= ncnt && i * prime[j] <= T; j++) {
                flag[i * prime[j]] = 1;
                if (i % prime[j] == 0) break;
            }
        }
        for (LL l = 1; l <= n; l = n / (n / l) + 1) {
            a[++m] = n / l;
            if (a[m] <= T) id1[a[m]] = m; else id2[n / a[m]] = m;
            g[m] = calc(a[m]);
        }
        for (int i = 1; i <= ncnt; i++)
            for (int j = 1; j <= m && (LL)prime[i] * prime[i] <= a[j]; j++)
                g[j] = g[j] - (LL)prime[i] * (g[ID(a[j] / prime[i])] - sum[i - 1]);
    }

    inline LL solve(LL x) {
        if (x <= 1) return x;
        return n = x, init(), g[ID(n)];
    }

}

int main() {
    LL n; scanf("%lld", &n);
    printf("%lld\n", Min25::solve(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
299
300
301
302
303
304
305
306
307
308
309
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int maxn = 2000000+100;

/******************************
f()函数中(31-37行) 填函数在质数幂次处的表达式
pow_sum()函数中(38-43行) 填幂和函数（如果需要更高次的话可以在这里添加）
202-205行按要求填写
f_p[][0/1/2/3/...]分别代表质数个数/质数和/质数平方和/质数三次方和/...根据自己需要添加
例：如果该函数在质数处表达式为f(p) = p^2+3*p+1
则表明需要质数个数/质数和/质数平方和
即f_p[][0],f_p[][1],f_p[][2]
******************************/

ll poww(ll a,ll b){
    ll res = 1;
    ll base = a;
    while(b){
        if(b&1){
            res *= base;
            //res %= mod;
        }
        base *= base;
        //base %= mod;
        b>>=1;
    }
    return res;
}

inline ll f(ll p,int e){
    if(p==1||e==0) return 1;
    ///return f(p^e)
    ll res = poww(p,e);
    return res*res+3*res+1;

}
ll pow_sum(ll n,int k){
    ///return sum(i^k),i from 1 to n.
    if(k==0) return n;
    if(k==1) return n*(n+1)/2;
    if(k==2) return n*(n+1)*(2*n+1)/6;
}
ll f_p[maxn][3];///F_prime(id(n/i))
ll n;
int n_2; ///(int)sqrt(n)
int n_3; ///(int)pow(n,1.0/3.0)
int n_6; ///(int)pow(n,1.0/6.0)
ll val_id[maxn]; ///give the id, return the id-th number like 'n/i' ,(val_id[1] = 1)
int val_id_num; ///how many numbers like 'n/i'
int val_id_num_3; ///how many numbers like 'n/i' below n/n_3;
int p[200000+100];
bool isp[maxn];
int p_sz_2; ///pi(n_2)
int p_sz_3; ///pi(n_3)
int p_sz_6; ///pi(n_6)
void init(){
    n_2 = (int)sqrt(n);
    n_3 = (int)pow(n,1.0/3.0);
    n_6 = (int)pow(n,1.0/6.0);
    val_id_num = 0;
    for(ll i=1;i<=n;){
        val_id[++val_id_num] = i;
        if(i==n) break;
        i = n/(n/(i+1));
    }
    memset(isp,1,sizeof isp);
    isp[1] = 0;
    for(int i=2;i<=n_2;i++){
        if(isp[i]){
            p[++p_sz_2] = i;
            if(i<=n_3) p_sz_3++;
            if(i<=n_6) p_sz_6++;
        }
        for(int j=1;j<=p_sz_2&&p[j]*i<=n_2;j++){
            isp[i*p[j]] = 0;
            if(i%p[j]==0) break;
        }
    }
}
inline int get_id(ll k){ ///give a number like 'n/i', return the id of it
    if(k>n_2) return val_id_num-n/k+1;
    else return k;
}
ll c[maxn];
int lowbit(int n){return n & (-n);}
void add(int x,ll d){
    while(x<maxn){
        c[x]+=d;
        x+=lowbit(x);
    }
}
ll sum(int x){
    ll ans=0;
    while(x){
        ans+=c[x];
        x-=lowbit(x);
    }
    return ans;
}

struct node{
    int k_max;
    ll val;
    ll f_val;
};
void update_bfs(int k,int type){
    queue<node> q;
    while(!q.empty()) q.pop();
    int e = 1;
    for(ll i=p[k];i<n/n_3;i*=p[k]){
        node st;
        st.k_max = k;
        st.val = i;
        if(type==-1)st.f_val = f(p[k],e);
        else st.f_val = poww(i,type);
        q.push(st);
        e++;
    }
    while(!q.empty()){
        node hd = q.front();
        q.pop();
        if((hd.val!=p[hd.k_max]&&type>=0)||type==-1) {//if(type==-1)cout << "****" << hd.val << "****" << hd.f_val << endl;
            ll w = n/hd.val;
            w = n/w;//cout << hd.val << "[" << w<<" , " << val_id[val_id_num] << "]" << endl;
            if(type==-1){
                add(get_id(w),hd.f_val);
                add(val_id_num+1,-1ll*hd.f_val);
            }
            else{
                add(get_id(w),-1ll*hd.f_val);
                add(val_id_num+1,hd.f_val);
            }
        }
        for(int i=hd.k_max+1;hd.val*p[i]<n/n_3&&i<=p_sz_2;i++){
            ll res = p[i];
            for(int e=1;;e++){
                if(hd.val*res<n/n_3){
                    node nxt;
                    nxt.k_max = i;
                    nxt.val = hd.val*res;
                    if(type==-1) nxt.f_val = hd.f_val*f(p[i],e);
                    else nxt.f_val = hd.f_val*poww(res,type);
                    q.push(nxt);
                }
                else break;
                res *= p[i];
            }
        }
    }
}
void get_f_p(ll n,int times){
    for(int i=1;i<=val_id_num;i++){
        for(int j=0;j<=times;j++){
            f_p[i][j] = pow_sum(val_id[i],j)-1;
        }
    }
    int now;
    //for(now=1;now<=p_sz_2;now++){
    for(now=1;p[now]<=n_6;now++){
        for(int j=val_id_num;j>=1;j--){
            ll w = val_id[j]/p[now];
            if(w<p[now]) break;
            ll val=1;
            for(int k = 0;k<=times;k++){
                f_p[j][k] = f_p[j][k] - val*(f_p[get_id(w)][k]-f_p[p[now-1]][k]);
                val *= p[now];
            }
        }
    }
    int nnow = now;
    int val = 1;
    for(int tt = 0;tt<=times;tt++){
        now = nnow;
        memset(c,0,sizeof c);
        add(1,f_p[1][tt]);
        for(int i=2;val_id[i]<n/n_3;i++){
            add(i,f_p[i][tt] - f_p[i-1][tt]);
        }
        for(;p[now]<=n_3;now++){
            for(int j=val_id_num;j>=1;j--){
                ll w = val_id[j]/p[now];
                if(val_id[j]<n/n_3) break;
                if(w<p[now]) break;
                if(w<n/n_3) f_p[j][tt] = f_p[j][tt] - (sum(get_id(w)) - sum(p[now-1]))*poww(p[now],tt);
                else f_p[j][tt] = f_p[j][tt] - (f_p[get_id(w)][tt]-sum(p[now-1]))*poww(p[now],tt);
            }
            update_bfs(now,tt);
        }
        for(int i=1;i<=val_id_num&&val_id[i]<n/n_3;i++)
            f_p[i][tt] = sum(i);
        for(;now<=p_sz_2;now++){
            for(int j=val_id_num;j>=1;j--){
                ll w = val_id[j]/p[now];
                if(val_id[j]<n/n_3) break;
                if(w<p[now]) break;
                f_p[j][tt] -= (f_p[get_id(w)][tt]-f_p[p[now-1]][tt])*poww(p[now],tt);
            }
        }
    }

    for(int i=1;i<=val_id_num;i++){
        ///if f(p) = p^2+3p+1,then write:f_p[i][0] = f_p[i][2] + 3*f_p[i][1] + f_p[i][0];
        f_p[i][0] = f_p[i][2] + 3*f_p[i][1] + f_p[i][0];
    }

}
ll F[2000000+100];
void get_f_3(ll n){ ///V(F_{pi(n^(1/3))+1},n)
    ll q = p[p_sz_3+1];
    for(int now=1;now<=val_id_num;now++){
        if(val_id[now]<q){
            F[now] = 1;
        }
        else if(val_id[now]<q*q){
            F[now] = 1+(f_p[now][0]-f_p[q-1][0]);
        }
        else{
            F[now] = 1+(f_p[now][0]-f_p[q-1][0]);
            for(int pp=p_sz_3+1;p[pp]<=(int)(sqrt(val_id[now]))&&pp<=p_sz_2;pp++){
                F[now] += f(p[pp],2) + (f(p[pp],1))*(f_p[get_id(val_id[now]/p[pp])][0]-f_p[get_id(p[pp])][0]);
            }
        }
    }
}
void get_f_6(ll n){ ///V(F_{pi(n^(1/6))+1},n)
    memset(c,0,sizeof c);
    add(1,F[1]);
    for(int i=2;val_id[i]<n/n_3;i++){
        add(i,F[i] - F[i-1]);
    }
    for(int k=p_sz_3;k>p_sz_6;k--){
        int now = val_id_num;
        for(;val_id[now]>=n/n_3;now--){
            int e = 1;
            ll _p = p[k];
            while(val_id[now]/_p){
                if(val_id[now]/_p>=n/n_3){
                    F[now] += F[get_id(val_id[now]/_p)]*f(p[k],e);
                }
                else{
                    F[now] += sum(get_id(val_id[now]/_p))*f(p[k],e);
                }
                _p *= p[k];
                e++;
            }
        }
        if(k==1) break;
        //cout << "******" << p[k] << "******" << n/n_3 << endl;
        update_bfs(k,-1);///bfs to update [lpf(i)==P{k-1}]f(i)
    }
    for(int i=1;i<=val_id_num&&val_id[i]<n/n_3;i++)
        F[i] = sum(i);
}
void get_f(ll n){
    for(int k=p_sz_6;k>=1;k--){
        for(int now = val_id_num;now>=1;now--){
            int e = 1;
            ll _p = p[k];
            while(val_id[now]/_p){
                    F[now] += F[get_id(val_id[now]/_p)]*f(p[k],e);
                _p *= p[k];
                e++;
            }
        }
    }
}
int main(){//n = 1000000000; //1e10:455052511,0.83s/0.58s 1e12:37607912018 9.224s/5.105s
    cin >> n;
    init();
    get_f_p(n,2);
    get_f_3(n);
    get_f_6(n);
    get_f(n);
    for(int i=1;i<=val_id_num;i++){
        cout << val_id[i] << " : " << F[i] << endl;
    }

}

```  

## 三, 唯一分解定理(算术基本定理)

暴力分解即可

代码没有

## 四, GCD/EXGCD

用途: **ax + by = gcd(a,b) 方程求解** 

```C++
#include <bits/stdc++.h>
using namespace std;
int a, b, x, y；

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

int main() {
    scanf("%d %d", &a, &b);
    EXGCD(a,b,x,y);
    printf("%d %d", x, y);
    return 0;
}

使用STL的代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#include <bits/stdc++.h>
using namespace std;
#define pir pair<int,int>    
int a, b, x, y;

pir EXGCD(int a,int b) {
    if (!b) return make_pair(1, 0);
    pir tmp = EXGCD(b, a%b);
    return make_pair(tmp.second ,tmp.first - a/b*tmp.second);
}

int main() {
    scanf("%d %d", &a, &b);
    pir ans = EXGCD(a,b);
    printf("%d %d",ans.first,ans.second);
    return 0;
}

五, 中国剩余定理(孙子定理)

1.CRT(保证$m_1,m_2…m_n$之间两两互质，求最小的x)

用途: 同余方程组求解!

image-20201006204034516

我们求出了任意一个满足条件的x之后，只需要对其模M就是最终的答案(因为MM是所有数的lcm)

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
#include<cstdio>
#define ll long long
//扩展欧几里得算法 
void gcd(ll a,ll b,ll &d,ll &x,ll &y)
{
    if(b==0){
        d=a;
        x=1,y=0;
    }
    else{//else不能省略 
        gcd(b,a%b,d,y,x);
        y-=(a/b)*x;
    }
}
//中国剩余定理 
ll China(int n,ll *m,ll *a)
{
    ll M=1,d,y,x=0;
    for(int i=0;i<n;i++) M*=m[i];
    for(int i=0;i<n;i++){
        ll w=M/m[i];
        gcd(m[i],w,d,d,y);
        x=(x+y*w*a[i])%M;
    }
    return (x+M)%M;
}

2. EXCRT($m_1,m_2…m_n$之间不再限制互质，求最小的x)

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
#include<cstdio>
#include<algorithm>
#include<iostream>
using namespace std;
typedef long long ll;
const int N=100000+10;
int n;
ll a[N],r[N];
ll exgcd(ll a,ll b,ll &x,ll &y){
    if(b==0){
        x=1,y=0;return a;
    }    
    ll ret=exgcd(b,a%b,y,x);
    y-=(a/b)*x;
    return ret;
}
ll excrt(){
    ll M=a[1],R=r[1],x,y,d;
    for(int i=2;i<=n;i++){
        d=exgcd(M,a[i],x,y);
        if((R-r[i])%d) return -1;
        x=(R-r[i])/d*x%a[i];
        R-=M*x;
        M=M/d*a[i];
        R%=M;
    }
    return (R%M+M)%M;
}
int main()
{
    while(scanf("%d",&n)!=EOF){
        for(int i=1;i<=n;i++){
            scanf("%lld%lld",&a[i],&r[i]);
        }
        printf("%lld\n",excrt());
    }
    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

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int S = 8;

//计算ret = (a*b)%c a, b, c < 2^63
ll mult_mod(ll a, ll b, ll c) {
    a %= c;
    b %= c;
    ll ans = 0;
    ll tmp = a;
    while (b) {
        if (b & 1) {
            ans += tmp;
            if (ans > c) {
                ans -= c;
            }
        }
        tmp <<= 1;
        if (tmp > c) {
            tmp -= c;
        }
        b >>= 1;
    }
    return ans;
}
//计算 ret = (a^n) % mod
ll pow_mod(ll a, ll n, ll mod) {
    ll ans = 1;
    ll tmp = a % mod;
    while (n) {
        if (n & 1) {
            ans = mult_mod(ans, tmp, mod);
        }
        tmp = mult_mod(tmp, tmp, mod);
        n >>= 1;
    }
    return ans;
}

//通过 a^(n-1) = 1(mod n) 判断n是不是素数
//n - 1 = x * 2 ^ t 中间使用二次判断
//是合数返回ture, 不一定是合数返回false

bool check(ll a, ll n, ll x, ll t) {
    ll ans = pow_mod(a, x, n);
    ll last = ans;
    for (int i = 1; i <= t; ++i) {
        ans = mult_mod(ans, ans, n);
        if (ans == 1 && last != 1 && last != n - 1) {
            return true;
        }
        last = ans;
    }
    if (ans != 1) {
        return true;
    }
    return false;
}

bool Miller_Rabin(ll n) { //是素数返回 ture, 否则返回false
    if (n < 2) {
        return false;
    }
    if (n == 2) {
        return true;
    }
    if ((n & 1) == 0) {
        return false;
    }
    ll x = n - 1;
    ll t = 0;
    t = (ll)log2(x & (-x));
    x = x / (x & (-x));
    srand(time(NULL));
    for (int i = 0; i < S; ++i) {
        ll a = rand() % (n - 1) + 1;//产生随机数a, 并控制其在0 ~ n-1之间
        if (check(a, n, x, t)) {
            return false;
        }
    }
    return true;
}
int main() {
    int t, ans;
    ll tem;
    scanf("%d", &t);
    ans = 0;
    while (t--) {
        scanf("%lld", &tem);
        if (Miller_Rabin(tem)) {
            ans++;
        }
    }
    cout << ans << endl;
    return 0;
}

增: Pallord_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
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 <bits/stdc++.h>
#define ll long long

#define INF 0x3f3f3f3f
#define maxn 10000+10
#define cle(a) memset(a, 0, sizeof(a))
const double eps = 1e-5;
using namespace std;

const int S = 30;//随机算法判定次数，S越大，判错概率越小
//计算 (a*b)%c.   a,b都是long long的数，直接相乘可能溢出的
//  a,b,c <2^63
ll mult_mod(ll a, ll b, ll c)
{
    a %= c;
    b %= c;
    ll ret = 0;
    while (b)
    {
        if (b & 1) { ret += a; ret %= c; }
        a <<= 1;
        if (a >= c)a %= c;
        b >>= 1;
    }
    return ret;
}
ll pow_mod(ll x, ll n, ll mod)
{
    if (n == 1)return x % mod;
    x %= mod;
    ll tmp = x;
    ll ret = 1;
    while (n)
    {
        if (n & 1) ret = mult_mod(ret, tmp, mod);
        tmp = mult_mod(tmp, tmp, mod);
        n >>= 1;
    }
    return ret;
}
//以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(ll a, ll n, ll x, ll t)
{
    ll ret = pow_mod(a, x, n);
    ll last = ret;
    for (int i = 1; i <= t; i++)
    {
        ret = mult_mod(ret, ret, n);
        if (ret == 1 && last != 1 && last != n - 1) return true;//合数
        last = ret;
    }
    if (ret != 1) return true;
    return false;
}
// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;
bool Miller_Rabin(ll n)
{
    if (n < 2)return false;
    if (n == 2)return true;
    if ((n & 1) == 0) return false;//偶数
    ll x = n - 1;
    ll t = 0;
    while ((x & 1) == 0) { x >>= 1; t++; }
    for (int i = 0; i < S; i++)
    {
        ll a = rand() % (n - 1) + 1;//rand()需要stdlib.h头文件
        if (check(a, n, x, t))
            return false;//合数
    }
    return true;
}
ll factor[100];//质因数分解结果（刚返回时是无序的）
int tol;//质因数的个数。数组小标从0开始
ll gcd(ll a, ll b)
{
    if (a == 0)return 1;//???????
    if (a < 0) return gcd(-a, b);
    while (b)
    {
        ll t = a % b;
        a = b;
        b = t;
    }
    return a;
}
ll Pollard_rho(ll x, ll c)
{
    ll i = 1, k = 2;
    ll x0 = rand() % x;
    ll y = x0;
    while (1)
    {
        i++;
        x0 = (mult_mod(x0, x0, x) + c) % x;
        ll d = gcd(y - x0, x);
        if (d != 1 && d != x) return d;
        if (y == x0) return x;
        if (i == k) { y = x0; k += k; }
    }
}
//对n进行素因子分解
void findfac(ll n)
{
    if (Miller_Rabin(n))//素数
    {
        factor[tol++] = n;
        return;
    }
    ll p = n;
    while (p >= n)p = Pollard_rho(p, rand() % (n - 1) + 1);
    findfac(p);
    findfac(n / p);
}
int main() {
    ll n;
    int t;
    cin >> t;
    while (t--) {
        cin >> n;
        if (n == 1) {
            //cout << "no" << '\n';
            
            continue;
        }
        tol = 0;
        findfac(n);
        sort(factor, factor + tol); //素因子排序
        
        for (int i = 0; i < tol; i++) { //输出所有素因子!!! 有重复的哦: 100 -> 2 2 5 5
            cout << factor[i] << ' ';
        }
        //质因子
    }
    return 0;
}

七, 乘法逆元

1.by EXGCD

O(logN)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 1 #include<bits/stdc++.h>
 2 typedef long long ll;
 3 ll exgcd(ll a,ll b,ll &x,ll &y) {
 4     if(!b) {
 5         x=1,y=0;
 6         return a;
 7     }
 8     ll res=exgcd(b,a%b,y,x);
 9     y-=a/b*x;       ///x=x1,y=x1-a/b*y1   x1,y1代表下一状态
10     return res;
11 }
12 int main()
13 {
14     ll a,p,x,y;  ///扩展欧几里得计算a的逆元（mod p）
15     scanf("%lld%lld",&a,&p);
16     ll d=exgcd(a,p,x,y);
17     printf(d==1?"%lld":"-1",(x+p)%p);///最大公约数不为1，逆元不存在，输出-1
18     return 0;
19 }

2.费马小定理

O(logP) 模为素数p的情况下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 1 #include<bits/stdc++.h>
 2 typedef long long ll;
 3 ll quickpowmod(ll a,ll b,ll mod) {
 4     ll ans=1;
 5     while(b) {
 6         if(b&1) ans=(ans*a)%mod;
 7         b>>=1;
 8         a=(a*a)%mod;
 9     }
10     return ans;
11 }
12 int main()
13 {
14     ll a,p;  ///费马小定理计算a的逆元（mod p）
15     scanf("%lld%lld",&a,&p);
16     ll inva=quickpowmod(a,p-2,p);
17     printf("%lld",inva);
18     return 0;
19 }

3.逆元打表

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

八, 矩阵快速幂

模板是求解 : F（n）= F（n-1）+F（n-2）+F（n-3）

具体题目具体分析, 推算矩阵

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD = 1000000007;

struct mat {
    ll a[3][3];
};

mat mat_mul(mat a,mat b) {
    mat res;
    memset(res.a,0,sizeof(res.a));
    for(int i = 0; i<=2 ; i++) {
        for(int j = 0; j<=2 ; j++) {
            for(int k = 0; k<=2 ; k++) {
                res.a[i][j] += ( a.a[i][k] * b.a[k][j])%MOD;
            }
        }
    }
    return res;
}

ll mat_pow(ll n) {
    mat c,res;
    memset(c.a,0,sizeof(c.a));
    memset(res.a,0,sizeof(res.a));
    c.a[0][0] = 1;
    c.a[0][1] = 1;
    c.a[0][2] = 1;
    c.a[1][0] = 1;
    c.a[2][1] = 1;
    for(int i=0; i<3 ; i++) {
        res.a[i][i] = 1;
    }
    while(n) {
        if(n&1) {
            res = mat_mul(res,c);
        }
        c = mat_mul(c,c);
        n >>= 1;
    }
    return (res.a[0][0])%MOD;
}

int main() {
    ll n;
    while(scanf("%lld",&n)!=EOF) {

//    n -= 1;

        ll ans = mat_pow(n);

        printf("%lld\n",ans);

    }
}

九, 卢卡斯定理

1.Lucas (模数p是质数)

给定$n,m,p (1 <= n,m,p <= 10^5)$ , 且p为素数
求 $C^n_{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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include<iostream>
using namespace std;
typedef long long LL;
const LL N=1e5+2; //1e5 !!!

LL a[N];
void init(LL p)
{
	a[1]=1;
	for(int i=2;i<=p;++i)a[i]=a[i-1]*i%p;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
	if(!b){
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b,y,x);
	y-=a/b*x;
}
LL ksm(LL x,LL n,LL mod)
{
	LL ans=1;
	while(n){
		if(n&1)ans=ans*x%mod;
		n>>=1;
		x=x*x%mod;
	}
	return ans;
}
LL C(LL n,LL m,LL p)
{
	if(n==m||m==0)return 1;
	if(n<m)return 0;
	if(m*2>n)m=n-m;						  /*C(n,m)=c(n,n-m)*/
	return a[n]*ksm(a[m]*a[n-m],p-2,p)%p;
    /*求(a[m]*a[n-m])在(mod p)意义下的乘法逆元*/										/*拓展欧几里得与费马小定理均可*/ 
	/*LL x,y;
	exgcd(a[m]*a[n-m],p,x,y);
	return (a[n]*x%p+p)%p;*/ 
}
LL lucas(LL n,LL m,LL p)
{
	if(!m)return 1;
	return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main()
{
	ios::sync_with_stdio(false);
	LL T,n,m,p;
	cin>>T;
	while(T--){
		cin>>n>>m>>p;
		init(p);
		cout<<lucas(n+m,m,p)<<endl;
	}
	return 0;
}

2. EX-Lucas (模数p不一定是质数的时候)

求$C^m_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
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
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N=1e5+9;
LL A[N],M[N];
LL ksm(LL x,LL n,LL mod)
{
	LL ans=1;
	while(n){
		if(n&1)ans=ans*x%mod;
		n>>=1,x=x*x%mod;
	}
	return ans;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
	if(!b)x=1,y=0;
	else exgcd(b,a%b,y,x),y-=a/b*x;
}
LL inv(LL a,LL p)
{
	LL x,y;
	exgcd(a,p,x,y);
	return (x+p)%p?x:x+p;
}
LL get(LL n,LL pi,LL p)									/*求(与pi互素后的n!)%M[i]*/ 
{
	if(!n)return 1;
	LL ans=1;
	if(n/p){											/*判断有无循环节 */ 
		for(LL i=2;i<=p;++i)if(i%pi)ans=ans*i%p;
		ans=ksm(ans,n/p,p);
	}
	for(LL i=2;i<=n%p;++i)if(i%pi)ans=ans*i%p;			/*循环节剩余部分*/ 
	return ans*get(n/pi,pi,p)%p;
}
LL exlucas(LL n,LL m,LL pi,LL p)						/*求A[i]*/ 
{
	LL nn=get(n,pi,p);									/*求(与pi互素后的n)%M[i]*/ 
	LL mm=get(m,pi,p);									/*求(m!与pi互素后的m!)%M[i]*/ 
	LL nm=get(n-m,pi,p);								/*求(与pi互素后的(n-m)!)%M[i]*/ 
	LL k=0;												/*含质因数pi的数量*/ 
	for(LL i=n;i;i/=pi)k+=i/pi;
	for(LL i=m;i;i/=pi)k-=i/pi;
	for(LL i=n-m;i;i/=pi)k-=i/pi;
	return nn*inv(mm,p)*inv(nm,p)*ksm(pi,k,p)%p;
}
LL crt(LL len,LL Lcm)
{
	LL ans=0;
	for(LL i=1;i<=len;++i){
		LL Mi=Lcm/M[i];
		ans=((ans+A[i]*inv(Mi,M[i])*Mi)%Lcm+Lcm)%Lcm;
	}
	return ans;
}
int main()
{
	ios::sync_with_stdio(false);
	LL n,m,P,num;
	while(cin>>n>>m>>P){
		if(n<m){
			cout<<0<<endl;
			continue;
		}
		num=0;
		memset(A,0,sizeof(A));
		memset(M,0,sizeof(M));
		for(LL x=P,i=2;i<=P;++i)
			if(x%i==0){
				M[++num]=1;
				while(x%i==0){
					M[num]*=i;
					x/=i;
				}
				A[num]=exlucas(n,m,i,M[num])%P;
			} 
		cout<<crt(num,P)<<endl;
	}
	return 0;
}

十, 多项式算法

FFT(根据复数, 答案是小数, 非精确)

NTT

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


#include <bits/stdc++.h>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 4000005, maxm = 100005, INF = 1000000000;

inline int read() {
    int out = 0, flag = 1; char c = getchar();
    while (c < 48 || c > 57) { if (c == '-') flag = -1; c = getchar(); }
    while (c >= 48 && c <= 57) { out = (out << 3) + (out << 1) + c - 48; c = getchar(); }
    return out * flag;
}

const int G = 3, P = (119 << 23) + 1;
int n, m, L, R[maxn];
int A[maxn], B[maxn];

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

void NTT(int* a, int f) {
    for (int i = 0; i < n; i++) if (i < R[i]) swap(a[i], a[R[i]]);
    for (int i = 1; i < n; i <<= 1) {
        int gn = qpow(G, (P - 1) / (i << 1));
        for (int j = 0; j < n; j += (i << 1)) {
            int g = 1;
            for (int k = 0; k < i; k++, g = 1ll * g * gn % P) {
                int x = a[j + k], y = 1ll * g * a[j + k + i] % P;
                a[j + k] = (x + y) % P; a[j + k + i] = (x - y + P) % P;
            }
        }
    }
    if (f == 1) return;
    int nv = qpow(n, P - 2); reverse(a + 1, a + n);
    for (int i = 0; i < n; i++) a[i] = 1ll * a[i] * nv % P;
}

int main() {
    n = read();
    m = read();
    for (int i = 0; i <= n; i++) A[i] = read();
    for (int i = 0; i <= m; i++) B[i] = read();
    m = n + m;                                //m = n + m!
    for (n = 1; n <= m; n <<= 1)
        L++;
    for (int i = 0; i < n; i++)
        R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
    NTT(A, 1);
    NTT(B, 1);
    for (int i = 0; i < n; i++)
        A[i] = 1ll * A[i] * B[i] % P;
    NTT(A, -1);
    for (int i = 0; i <= m; i++)
        printf("%d ", A[i]);
    return 0;
}

十一, 康托展开/逆康托展开

康托展开, 某排列映射一个值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//对前 10 个自然数(0 ~ 9)的阶乘存入表
//以免去对其额外的计算
const int fact[10] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};
/**
 * @brief 康拓展开
 *
 * @param[in] permutation 输入的一个全排列
 * @param[out] num 输入的康拓映射，即是第几个全排列
 */
int contor(const vector<int>& permutation) {
    int num = 0;
    int len = permutation.size();
    for (int i = 0; i < len; ++i) {
        int cnt = 0; // 在 i 之后，比 i 还小的有几个
        for (int j = i + 1; j < len; ++j)
            if (permutation[i] > permutation[j]) ++cnt;
        num += cnt * fact[len - i - 1];
    }
    return num + 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
//对前 10 个自然数(0 ~ 9)的阶乘存入表
//以免去对其额外的计算
const int fact[10] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};
/**
 * @brief 逆康拓展开
 *
 * @param[in] bits 给定全排列的使用数字个数
 * @param[in] num 给定全排列的次位
 * @param[out] permutation 输出对应的全排列
 */
vector<int> revContor(int bits, int num) {
    num = num - 1; //有 num - 1 个排列比目标序列要小
    vector<bool> vis(bits + 1, false);
    vector<int> permutation(bits, -1);

    int n, residue = num;
    for (int i = 0; i < bits; ++i) {
        n = residue / (fact[bits - i - 1]);
        residue = residue % (fact[bits - i - 1]);

        for (int j = 1; j <= bits; ++j) {
            if (!vis[j] && !(n--)) {
                vis[j] = true;
                permutation[i] = j;
                break;
            }
        }
    }
    return permutation;
}

十二, 日期计算

吉姆拉尔森公式

1
2
3
4
5
6
7
8
int CaculateWeekDay(int y, int m, int d) { //年月日, 注意: 1/2月要当成上一年的13/14月
    if (m == 1 || m == 2) {
        m += 12;
        y--;
    }
    int iWeek = (d + 2 * m + 3 * (m + 1) / 5 + y + y / 4 - y / 100 + y / 400) % 7;
    return iWeek; // (0 - 6)
}

十三, 容斥原理

递归/暴力

若在n较大的情况下，应该采用递归进行计算，位运算计算循环次数太多，反而很慢。
在集合1，2，3…600,中，取出能被2，3，5，整除的数的个数

位操作:

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
#include<bits/stdc++.h>
using namespace std;

int p[] = { 2,3,5 };
int cnt = 3;

int cal(int n = 600) {
    int res = 0;
    for (int i = 1; i < (1 << cnt); i++) {
        int t = i, tmp = 1, k = 0;
        int len = 0;
        while (t) {
            if (t & 1) {
                tmp *= p[k];
                len++;
            }
            t >>= 1;
            k++;
        }
        if (len & 1) res += n / tmp;
        else res -= n / tmp;
    }
    return res;
}


int main() {
    cout << cal() << endl;
    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
#include<bits/stdc++.h>
using namespace std;

int a[] = { 2,3,5 };

int b = 600;
int sum = 0;
int n = 3;


void dfs(int i, int num, int x, int mu) {  //i表示第几个元素，num表示现在一共用了几个，x表示最多能用几个，mu表示当前取了的数的乘积
    if (num == x) {
        sum += b / mu;                //b就表示600
        return;
    }

    if (i == n) return;             //一共有n个，比如上面的2，3，5，一共有3个，
    dfs(i + 1, num + 1, x, mu * a[i]);    //a中存储的就是2，3，5，然后取或者不取
    dfs(i + 1, num, x, mu);
}

int rong() {
    int s = 0;
    for (int i = 1; i <= n; i++) {
        sum = 0;
        dfs(0, 0, i, 1);
        if (i & 1) s += sum;        //容斥定理
        else s -= sum;
    }
    return s;                    //s为能被2，3，5整除的数的个数。
}


void dfs(int i, int mu, int num) {
    mu *= a[i];
    if (num % 2) sum += b / mu;
    else sum -= b / mu;
    for (int j = i + 1; j < 3; j++) {
        dfs(j, mu, num + 1);
    }
}

int main() {
    printf("%d\n", rong());
    sum = 0;
    for (int i = 0; i < 3; i++) {
        dfs(i, 1, 1);
    }
    printf("%d\n", sum);
    return 0;
}

从 m 中颜色中选择 k 种, 给 n 个盒子染色, 相邻盒子颜色不同

打表求 组合数(取模)

用不超过 k 种颜色涂完 n 个盒子的方案数为 $ x_k=k(k−1)^{n−1} $

n 个盒子，用 m 种颜色中选出 正好 k 种颜色，且相邻盒子颜色不同的方案数为:
$$
ans = C_{m}^{k} t m p=C_{m}^{k} x_{k}-C_{m}^{k-1} x_{k-1}+C_{m}^{k-2} x_{k-2}+\ldots+(-1)^{k} x_{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
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
#include<stdio.h>
#include<cstdio>
#include<iostream>
using namespace std;

#define ll long long

const ll mod=1e9+7;
const int maxn=1e6+10;

int T;
ll p[maxn], invp[maxn], inv[maxn];

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

void init(){
    p[0] = 1;
    for(int i = 1; i <= 1000000; i++)
        p[i] = p[i - 1] * i % mod;
    invp[1000000] = qp(p[1000000], mod - 2);
    for(int i = 1000000 - 1; i>=0; i--)
        invp[i] = invp[i + 1] * (i + 1) % mod;
    inv[1] = 1;
    for(int i = 2; i <= 1000000; i++)
        inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
}

ll C(ll n, ll m){
    return p[n] * invp[m] % mod * invp[n - m] % mod;
}

int main()
{
    init();
    scanf("%d",&T);
    while(T--)
    {
        ll n,m,k;
        scanf("%lld%lld%lld",&n,&m,&k);
        ll Cmk = 1;
        for(int i=1; i<=k; i++)
            Cmk = Cmk * (m - k + i) % mod * inv[i] % mod;
        ll tmp = 0;
        for(int i=0; i<k; i++)
            if (i % 2 == 0)
                tmp = (tmp + C(k, k-i) * (k - i) % mod * qp(k - i - 1, n - 1) % mod + mod) % mod;
            else
                tmp = (tmp - C(k, k-i) * (k - i) % mod * qp(k - i - 1, n - 1) % mod + mod) % mod;
        ll ans = Cmk * tmp % mod;
        printf("%lld\n", ans);
    }
    return 0;
}

十四, 欧拉降幂

输入n, m
长度n序列, m是mod
q个查询, 每次输入l, r

求:
$$ w[l]^{w[l+…]^{w[r]}}$$

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

#include <bits/stdc++.h>
//#include <hash_map>
//#include <unordered_map> //没快多点
using ll = long long;
using namespace std;

#define MOD(a,b) a>=b?a%b+b:a  //重定义取模，按照欧拉定理的条件
const int maxn = 1e6 + 10;

int is_prime[maxn];
int prime[maxn];
int phi[maxn];
int cnt;

void init() {
    cnt = 0; //统计素数个数
    is_prime[1] = 1; //0代表是素数, 1不是素数
    phi[1] = 1;
    for (int i = 2; i <= maxn - 10; i++) {
        
        if (!is_prime[i]) {
            prime[++cnt] = i; //i是素数, prime从下标[1]开始有记录
            phi[i] = i - 1;
        }
        for (int j = 1; j <= cnt; j++) { // 1 到 cnt !
            if (prime[j] * i >= maxn)
                break;
            is_prime[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else {
                phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            }
        }
    }
}

unordered_map<int, int> ma;
int n;
ll w[maxn], m, l, r, x;

ll qpow(ll a, ll b, ll mod) { //快速幂
    ll res = 1;
    while (b) {
        if (b & 1) res = MOD(res * a, mod);
        a = MOD(a * a, mod);
        b >>= 1;
    }
    return res;
}

ll getphi(ll x) {  //求互质个数
    if (x <= maxn - 10) {
        return phi[x];
    }
    ll ans = x;
    if (ma.count(x)) return ma[x];
    for (int i = 1; 1ll*prime[i] * prime[i] <= x; i++) {
        if (x % prime[i] == 0) {
            ans = ans / prime[i] * (prime[i] - 1);
            while (x % prime[i] == 0) {
                x /= prime[i];
            }
        }
        /*
        if (x % i == 0) {
            ans = ans / i * (i - 1);
            while (x % i == 0)
                x /= i;
        }
        */
    }
    if (x > 1) ans = ans / x * (x - 1); //*****
    ma[x] = ans;
    return ans;
}
ll solve(int ceng, ll mod) { //每层的幂是w[i]
    if (ceng == r || mod == 1)
        return MOD(w[ceng], mod);
    return
        qpow(w[ceng], solve(ceng + 1, getphi(mod)), mod);
}

int main() {
    init();
    //cin >> n >> m;
    scanf("%d %lld", &n, &m);
    for (int i = 1; i <= n; i++) {
        scanf("%d", &w[i]);
    }
    int q;
    //cin >> q;
    scanf("%d", &q);
    while (q--) {
        //cin >> l >> r;
        scanf("%lld %lld", &l, &r);
        //r = rr;
        cout << solve(l, m) % m << '\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
const double EPS = 1E-9;
int n;
vector<vector<double> > a(n, vector<double>(n));

double det = 1;
for (int i = 0; i < n; ++i) {
  int k = i;
  for (int j = i + 1; j < n; ++j)
    if (abs(a[j][i]) > abs(a[k][i])) k = j;
  if (abs(a[k][i]) < EPS) {
    det = 0;
    break;
  }
  swap(a[i], a[k]);
  if (i != k) det = -det;
  det *= a[i][i];
  for (int j = i + 1; j < n; ++j) a[i][j] /= a[i][i];
  for (int j = 0; j < n; ++j)
    if (j != i && abs(a[j][i]) > EPS)
      for (int k = i + 1; k < n; ++k) a[j][k] -= a[i][k] * a[j][i];
}

cout << det;

十六, 杜教BM

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long long ll;
typedef vector<int> VI;
const int maxn = 10005;
const ll mod = 1e9 + 7;
const int INF = 0x3f3f3f3f;
const double eps = 1e-9;
ll fast_mod(ll a, ll n, ll Mod)
{
    ll ans = 1;
    a %= Mod;
    while (n)
    {
        if (n & 1)
            ans = (ans * a) % Mod;
        a = (a * a) % Mod;
        n>>= 1;
    }
    return ans;
}
namespace linear_seq
{
ll res[maxn], base[maxn], num[maxn], md[maxn]; // 数组大小约 10000
vector<int> vec;
void mul(ll *a, ll *b, int k)
{
    for (int i = 0; i <2 * k; i++)
        num[i] = 0;
    for (int i = 0; i < k; i++)
    {
        if (a[i])
            for (int j = 0; j < k; j++)
                num[i + j] = (num[i + j] + a[i] * b[j]) % mod;
    }
    for (int i = 2 * k - 1; i>= k; i--)
    {
        if (num[i])
            for (int j = 0; j <vec.size(); j++)
                num[i - k + vec[j]] = (num[i - k + vec[j]] - num[i] * md[vec[j]]) % mod;
    }
    for (int i = 0; i < k; i++)
        a[i] = num[i];
}
ll solve(ll n, VI a, VI b)
{
    ll ans = 0, cnt = 0;
    int k = a.size();
    assert(a.size() == b.size());
    for (int i = 0; i < k; i++)
        md[k - 1 - i] = -a[i];
    md[k] = 1;
    vec.clear();
    for (int i = 0; i < k; i++)
        if (md[i])
            vec.push_back(i);
    for (int i = 0; i < k; i++)
        res[i] = base[i] = 0;
    res[0] = 1;
    while ((1LL << cnt) <= n)
        cnt++;
    for (int p = cnt; p>= 0; p--)
    {
        mul(res, res, k);
        if ((n>> p) & 1)
        {
            for (int i = k - 1; i>= 0; i--)
                res[i + 1] = res[i];
            res[0] = 0;
            for (int j = 0; j < vec.size(); j++)
                res[vec[j]] = (res[vec[j]] - res[k] * md[vec[j]]) % mod;
        }
    }
    for (int i = 0; i < k; i++)
        ans = (ans + res[i] * b[i]) % mod;
    if (ans < 0)
        ans += mod;
    return ans;
}
VI BM(VI s)
{
    VI B(1, 1), C(1, 1);
    int L = 0, m = 1, b = 1;
    for (int i = 0; i < s.size(); i++)
    {
        ll d = 0;
        for (int j = 0; j < L + 1; j++)
            d = (d + (ll)C[j] * s[i - j]) % mod;
        if (d == 0)
            m++;
        else if (2 * L <= i)
        {
            VI T = C;
            ll c = mod - d * fast_mod(b, mod - 2, mod) % mod;
            while (C.size() < B.size() + m)
                C.push_back(0);
            for (int j = 0; j < B.size(); j++)
                C[j + m] = (C[j + m] + c * B[j]) % mod;
            L = i + 1 - L, B = T, b = d, m = 1;
        }
        else
        {
            ll c = mod - d * fast_mod(b, mod - 2, mod) % mod;
            while (C.size() < B.size() + m)
                C.push_back(0);
            for (int j = 0; j < B.size(); j++)
                C[j + m] = (C[j + m] + c * B[j]) % mod;
            m++;
        }
    }
    return C;
}
int gao(VI a, ll n)
{
    VI c = BM(a);
    c.erase(c.begin());
    for (int i = 0; i < c.size(); i++)
        c[i] = (mod - c[i]) % mod;
    return solve(n, c, VI(a.begin(), a.begin() + c.size()));
}
} // namespace linear_seq
int main()
{
    //freopen("in.txt", "r", stdin);
    // 填数字的时候带上模数之后的
    ll t, n;
    scanf("%d", &t);
    while (t--)
    {
        scanf("%lld", &n);
        // 在 vi 中填入自己推的数 8 个以上 越多越好
        printf("%lld\n", linear_seq::gao(VI{3, 9, 20, 46, 106, 244, 560, 1286, 2956, 6794}, n-1));
    }
    return 0;
}

十七, 拉格朗日插值

高数的应用: 根据导数, 控制函数大小
这里的应用: 给出n个点P(xi, yi)，将过这个点的最多n - 1次的多项式记为f(x)，求f(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
#include <cstdio>

const int maxn = 2010;
using ll = long long;
ll mod = 998244353;
ll n, k, x[maxn], y[maxn], ans, s1, s2;

ll powmod(ll x, ll n) {
  ll ret = 1ll;
  while (n) {
    if (n & 1) ret = ret * x % mod;
    x = x * x % mod;
    n >>= 1;
  }
  return ret;
}

ll inv(ll x) { return powmod(x, mod - 2); }

int main() {
  scanf("%lld%lld", &n, &k);
  for (int i = 1; i <= n; i++)
    scanf("%lld%lld", x[i], y[i]);
  for (int i = 1; i <= n; i++) {
    s1 = y[i] % mod;
    s2 = 1ll;
    for (int j = 1; j <= n; j++)
      if (i != j) s1 = s1 * (k - x[j]) % mod, s2 = s2 * (x[i] - x[j]) % mod;
    ans += s1 * inv(s2) % mod;
  }
  printf("%lld\n", (ans % mod + mod) % mod);
  return 0;
}

十八, Catalan 数列

常常用来求解这样一类问题:

  1. 有2n个人排成一行进入剧场。入场费5元。其中只有n个人有一张5元钞票，另外n人只有10元钞票，剧院无其它钞票，问有多少中方法使得只要有 10 元的人买票，售票处就有 5 元的钞票找零?
  2. 在圆上选择2n个点，将这些点成对连接起来使得所得到的n条线段不相交的方法数？
  3. n个不同的数依次进栈，求不同的出栈结果的种数?

这一类, 有两种选择, 每种选择只能选择有限次数, 并且两种选择的可选择次数相同, 问多少种方案的问题.

1
2
3
4
5
6
7
8
9
10
11
12
#include <iostream>
using namespace std;
int n;
long long f[25];
int main() {
  f[0] = 1;
  cin >> n;
  for (int i = 1; i <= n; i++) f[i] = f[i - 1] * (4 * i - 2) / (i + 1);
  // 这里用的是常见公式2
  cout << f[n] << endl;
  return 0;
}

十九, 常见的期望和方差

0-1分布

1 次实验
P{x = 0} = p
P{x = 1} = 1 - p

E(x) = p;
D(x) = pq;

二项分布 B(n, p)

n次独立重复实验
P{x = k} = p

E(x) = np;
D(x) = npq;

几何分布

在n次独立重复实验中, 前k-1次皆失败，第k次成功的概率
P{x = k} = p

E(x) = 1/p;
D(x) = (1 - p) / (p * p);