oi-wiki/math

szu acm

筛法

素数筛法

如果我们想要知道小于等于 n 有多少个素数呢？

埃拉托斯特尼筛法(埃氏筛)

对于任意一个大于 1 的正整数 n，那么它的 x 倍就是合数（x > 1）。
如果我们从小到大考虑每个数，然后同时把当前这个数的所有（比自己大的）倍数记为合数，那么运行结束的时候没有被标记的数就是素数了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 𝑂 (𝑛 log log 𝑛)
vector<int> prime;
bool is_prime[N];

void Eratosthenes(int n)
{
    is_prime[0]=is_prime[1]=0;
    for(int i=2;i<=n;i++)   is_prime[i]=1;
    for(int i=2;i<=n;i++)
    {
        if(is_prime[i])
        {
            prime.push_back(i);
            for(int j=i*i;j<=n;j+=i)    is_prime[j]=0;
        }
    }
}

线性筛法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// O(n)
vector<int> pri;
bool not_prime[N];

void pre(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!not_prime[i])   pri.push_back(i);
        for(int pri_j : pri)
        {
            if(i*pri_j>n)   break;
            not_prime[i*pri_j]=1;
            if(i%pri_j==0)  break;
        }
    }
}

q

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
const ll N =1e8+10;
const ll mod = 1e9+7;

vector<int> pri;
bool not_prime[N];

void pre(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!not_prime[i])   pri.push_back(i);
        for(int pri_j : pri)
        {
            if(i*pri_j>n)   break;
            not_prime[i*pri_j]=1;
            if(i%pri_j==0)  break;
        }
    }
}

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    
    ll n,q;
    cin>>n>>q;
    pre(n);
    while(q--)
    {
        ll k;
        cin>>k;
        cout<<pri[k-1]<<"\n";
    }
    return 0;
}

扩展欧几里得(exgcd)

1
2
3
4
5
6
7
8
9
10
11
12
13
//与欧几里得算法一样，复杂度为 O(log n)
int exgcd(int a,int b,int &x,int &y)
{
    if(!b)  //b=0时,a=gcd(a,b)
    {
        x=1,y=0;
        return a;
    }
    int gcd=exgcd(b,a%b,y,x);   //交换实参x,y
    y-=(a/b)*x;
    return gcd;
}

乘法逆元

如果一个线性同余方程 𝑎𝑥 ≡ 1 (mod 𝑏)，则 𝑥 称为 𝑎 mod 𝑏 的逆元，记作 𝑎 ^ −1

费马小定理求逆元

费马小定理：𝑏 为素数且 gcd (𝑎, 𝑏) = 1 时，有 𝑎 ^ 𝑏−1 ≡ 1 (mod 𝑏)

根据定义 𝑎 ^ −1 = 𝑎 ^ 𝑏−2 ，使用快速幂求解即可

扩展欧几里得求逆元

• 要求 gcd (𝑎, 𝑏) = 1。
• 利用 exgcd 解线性方程 𝑎𝑥 + 𝑏𝑦 = 1。
• 解得 𝑥 显然满足 𝑎𝑥 ≡ 1 (mod 𝑏)
• 根据定义 𝑎 ^ −1 = 𝑥，考虑到 exgcd 的解可能为负数，故最终结果一般写成 ((𝑥%𝑏)+𝑏) % b

线性递推求逆元

• 用于求解 1,2, … , 𝑛 中每个数关于 𝑝 的逆元
• 显然有 1 ^ −1 = 1 mod 𝑝
• 对于求 𝑖 ^ −1，令 𝑘 = [𝑝/𝑖], 𝑗 = 𝑝 mod 𝑖，则有 𝑝 = 𝑘𝑖 + 𝑗，
故能得到 𝑘i+𝑗 ≡ 0 (mod 𝑝)
• 左右同乘 𝑖 ^ −1 × 𝑗 ^ −1，得到 𝑘 * 𝑗 ^−1 + 𝑖 ^−1 ≡ 0 (mod 𝑝)
• 带入 𝑗 = 𝑝 mod 𝑖，得到 𝑖 ^ −1 ≡ −[𝑝/𝑖](𝑝 mod 𝑖)^−1 (mod 𝑝)
• 由于 𝑝 mod 𝑖 < 𝑖，故可通过递推求得 𝑖 ^−1 ：
1 , 𝑖 = 1
𝑖 ^−1 =
−[𝑝/i](𝑝 mod 𝑖)^ −1 (mod 𝑝)

线性求任意n个数的逆元

求任意给定的 𝑛 个数 𝑎1, … , 𝑎𝑛的逆元；
利用前缀积思想，具体参考
szu acm
计算 𝑛 个数的前缀积，记为 𝑠[i]
使用快速幂或扩展欧几里得法计算 𝑠[n] 的逆元，记为 𝑠𝑣[n]
𝑠𝑣[n] 乘上 𝑎𝑛 会与 𝑎𝑛 的逆元抵消,得到 𝑎1𝑎2 … 𝑎𝑛−1 的逆元,即 𝑠𝑣[𝑛−1]
𝑎𝑖^(−1) = 𝑠𝑣[i] × 𝑠[𝑖−1]

q

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
const int N =5e6+10;
const ll mod = 1e9+7;

inline ll read()
{
    ll x=0,tag=1;
    char c=getchar();
    for(; !isdigit(c);c=getchar())
    {
        if(c=='-')  tag=-1;
    }
    for(; isdigit(c);c=getchar())
    {
        x=(x<<1)+(x<<3)+c-'0';
    }
    return x*tag;
}

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

signed main()
{
    ll n,p,k;
    n=read(),p=read(),k=read();
    vector<ll> a(n+1),s(n+1),sv(n+1);

    s[0]=1;
    for(ll i=1;i<=n;i++)
    {
        a[i]=read();
        s[i]=(s[i-1]*a[i])%p;   //前缀积 s[i]=a[1]*a[2]*...*a[i]
    }
    //扩展欧几里得求逆元
    ll t;
    exgcd(s[n],p,sv[n],t);     //求sv[n]=s[i]^(-1)
    sv[n]=(sv[n]%p+p)%p;
    //逆元抵消
    for(ll i=n-1;i;i--)
    {
        sv[i]=sv[i+1]*a[i+1]%p; //逆元的前缀积
    }

    ll nowk=k,ans=0;
    for(ll i=1;i<=n;i++)
    {
        (ans += nowk * s[i-1] % p * sv[i] % p) %=p;
        (nowk *= k) %= p;
    }
    cout<<ans<<"\n";

    return 0;
}

龟速幂

q

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
const int N =5e6+10;

ll x,y,a,b,n,m,l;

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

inline ll mul(ll a,ll b,ll mod)  //龟速乘  求a*b
{
    ll ret=0;
    while(b)
    {
        if(b&1) ret=(ret+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return ret%mod;
}

inline ll Pow(ll a,ll b,ll mod) //快速幂  求a^b
{
    ll ret=1;
    while(b)
    {
        if(b&1) ret=mul(ret,a,mod); // ret*a
        a=mul(a,a,mod); // a*a
        b>>=1;
    }
    return ret%mod;
}

int main()
{
    cin>>n>>m>>l;
//扩展欧几里得求逆元
    exgcd(2,n+1,x,y);       
    x=(x%(n+1)+n+1)%(n+1);  //2 mod(n+1)的逆元

    x=Pow(x,m,n+1);     //x^m
    cout<<mul(l,x,n+1);  // l*x
    return 0;
}

Lucas定理

oi-wiki
C(n,m) mod p = C(n/p , m/p)*C(n mod p , m mod p) mod p

【模板】卢卡斯定理/Lucas定理

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
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
const int N =1e5+10;

ll a[N];
int n,m,p;

ll pow(ll y,int z,int p)    //快速幂 y^z
{
    y%=p;
    ll ans=1;
    while(z)
    {
        if(z&1) ans=ans*y%p;
        z>>=1;
        y=y*y%p;
    } 
    return ans;
}

ll C(ll n,ll m)     //组合C(n,m)
{
    if(m>n) return 0;
    return ((a[n]*pow(a[m],p-2,p))%p*pow(a[n-m],p-2,p)%p);    //费马小定理
}

ll Lucas(ll n,ll m)
{
    if(!m)  return 1;
    return C(n%p,m%p)*Lucas(n/p,m/p)%p;
}

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    
    int t;
    cin>>t;
    while(t--)
    {  
        cin>>n>>m>>p;
        a[0]=1;
        for(int i=1;i<=p;i++)   a[i]=(a[i-1]*i)%p;  //阶乘
        cout<<Lucas(n+m,n)<<endl;
    }

    return 0;
}