抱歉，您的浏览器无法访问本站
本页面需要浏览器支持（启用）JavaScript
了解详情 >

lfyszyの小窝

我是源点，你是终点。我们之间有负权环。|| Joker.

博客

分类

标签

归档

友链

关于

啊吧

开往

数论全家桶

算法

发布于：2024年1月24日

次浏览

闲话

我谔谔，第一稿写了一大部分没有保存。

质数 & 约数

质数判定

试除法

没啥好说，会的都会，复杂度 O(n)O(\sqrt{n})

1
2
3
4
5
6
7
bool isprime(int x)
{
    if(x == 1) return false;
    for(int i = 2; i <= n / i; i ++)
        if(x % i == 0) return false
    return true;
}

Fermat 素性检验

通过反着使用费马小定理，多次换底增加正确性。

但是可以被卡迈尔数完美卡掉。

Miller Rabin 素性检验

【朝夕的ACM笔记】数论-Miller Rabin素数判定 - 知乎 (zhihu.com)

在 Fermat 素性检验基础上，引入二次探测定理（对于质数 pp，x2≡1(modp)x^2\equiv1\pmod p 的解仅可能为 x≡1(modp)x\equiv 1\pmod p 或 x≡p−1(modp)x\equiv p-1\pmod p）以增加正确性。

对于 Fermat 素性检验中，ap−1≡1(modp)a^{p-1}\equiv 1\pmod p，可以解出 ap−12≡1(modp)a^{\frac{p-1}{2}}\equiv 1\pmod p 或 ap−12≡p−1(modp)a^{\frac{p-1}{2}}\equiv p-1\pmod p.

以此类推，即可验证 pp 是否为质数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
int qmul(int a, int b, int p) {return (a * b - (int)((long double) a / p * b) * p + p) % p;}//O(1)龟速乘
int qpow(int a, int b, int p) {int res = 1; while(b) {if(b & 1) res = qmul(res, a, p); a = qmul(a, a, p); b >>= 1;} return res;}
bool miller_rabin(int p)
{
	if(p < 3 || p % 2 == 0) return p == 2;
	int u = p - 1, t = 0;
	while(u % 2 == 0) u /= 2, t ++;
	int used[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
	for(auto a : used)
	{
		if(a >= p) break;
		int v = qpow(a, u, p);
		if(v == 1 || v == p - 1) continue;
		for(int j = 1; j <= t; j ++)
		{
			v = qmul(v, v, p);
			if(v == p - 1 && j != t) {v = 1; break;}
			if(v == 1) return false;
		}
		if(v != 1) return false;
	}
	return true;
}

分解质因数

试除法

没啥好说，会的都会，复杂度 O(n)O(\sqrt{n}).

1
2
3
4
5
6
7
8
9
10
11
12
void divdie(int n)
{
    for(int i = 2;i <= n / i; i ++)
        if(n % i == 0)
        {
            int s = 0;
            while(n % i == 0)
                n /= i, s ++;
            cout << i << " " << s << endl;
        }
    if(n > 1) cout << n << " " << 1 << endl;
}

筛法

比较鸡肋……复杂度比上一个稍快一点点。

Pollard Rho 大数分解

前置：生日悖论（简单来讲就是你选两个来的答案的概率会比一个高）。

生成两个随机数，尝试它们的差与被分解数的 gcd.

利用一种特殊的随机函数 f(x)=(x2+c)mod  pf(x)=(x^2+c)\mod p，由于模数的限制，此函数一定会出现循环，有因为是混循环，长得像 ρ\rho 所以得名 /kk.

实现时，由于 gcd 的 O(log⁡n)O(\log{n}) 的复杂度，我们考虑将差累乘起来，每隔一段时间 gcd 一下，这样便可以大大降低复杂度（由于前人的口胡经验，我们每 127 次 gcd 一下）。

具体可以搭配 Miller Rabin 与 dfs 实现分解质因数等。

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
int qmul(int a, int b, int p) {return (a * b - (int)((long double) a / p * b) * p + p) % p;}//O(1)龟速乘
int f(int x, int c, int p) {return (qmul(x, x, p) + c) % p;}
#define gcd __gcd
int pollar_rho(int p)
{
    if(p % 2 == 0) return 2;
#define step 127
    bool isfind = false;
    int c, r1, r2, pd = 1, tp, res = 0;
    while(!isfind)
    {
        c = rand() % (p - 1) + 1, r1 = r2 = 0;
        do
        {
            for(int i = 1; i <= step; i ++)
            {
                r1 = f(r1, c, p), r2 = f(f(r2, c, p), c, p);
                if(r1 == r2 || !(tp = qmul(pd, abs(r2 - r1), p))) break;
                pd = tp;
            }
            res = gcd(pd, p);
            if(res > 1) isfind = true;
        } while(!isfind && r1 != r2);
    }
    return res;
}

比较乱搞的东西qwq

设 x=q1c1q2c2…qkckx = q_1^{c_1}q_2^{c_2}\dots q_k^{c_k}，则有：

约数个数

(c1+1)(c2+1)…(ck+1)(c_1+1)(c_2+1)\dots (c_k+1)

约数之和

(p10+p11+⋯+p1c2)(p20+p21+⋯+p1c2)…(pk0+pk1+⋯+p1ck)(p_1^0+p_1^1+\dots +p_1^{c_2})(p_2^0+p_2^1+\dots +p_1^{c_2})\dots (p_k^0+p_k^1+\dots +p_1^{c_k})

过于简单（

欧拉函数

φ(n)\varphi(n) 表示 1∼n−1{1}\sim n-1 中，与 nn 互质的数量。

将 nn 分解为 n=q1c1q2c2…qkckn = q_1^{c_1}q_2^{c_2}\dots q_k^{c_k}

φ(n)=n(1−1p1)(1−1p2)…(1−1pk)\varphi(n)=n(1-\dfrac{1}{p_1})(1-\dfrac{1}{p_2})\dots(1-\dfrac{1}{p_k})

那么考虑 nn 与 n′=n×pn'=n\times p ，pp 为质数，有：

φ(n′)={φ(n)×p×(1−1p)=φ(n)×(p−1)  (p∤n)φ(n)×p  (p∣n)\varphi(n')=\begin{cases} \varphi(n)\times p\times (1-\dfrac{1}{p})=\varphi(n)\times (p-1)~~(p\nmid n)\\ \varphi(n)\times p~~(p\mid n) \end{cases}

那么我们就可以在标记的过程中递推求得欧拉函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
const int N = /*max num*/;
bool st[N];
int phi[N];
vector<int> p;
void getprime(int n)
{
    for(int i = 2; i <= n; i ++)
    {
        if(!st[i]) p.push_back(i), phi[i] = i - 1;
        for(auto v : p)
        {
            if(v * i > n) break;
            st[v * i] = 1;
            if(i % v == 0) {phi[i * v] = phi[i] * v; break;}
            else phi[i * v] = phi[i] * (v - 1);
        }
    }
}

同余

性质

  1. a≡a(modp)a\equiv a\pmod p
  2. a≡b(modp)a\equiv b\pmod p 时，b≡a(modp)b\equiv a\pmod p
  3. a≡b(modp)a\equiv b\pmod p 且 b≡c(modp)b\equiv c\pmod p时，a≡c(modp)a\equiv c\pmod p
  4. a≡b(modp)a\equiv b\pmod p 时，a+c≡b+c(modp)a+c\equiv b+c\pmod p
  5. a≡b(modp)a\equiv b\pmod p 时，ac≡bc(modp)ac\equiv bc\pmod p
  6. a≡b(modp)a\equiv b\pmod p 时，ac≡bc(modp)a^c\equiv b^c\pmod p
  7. ac≡bc(modp)ac\equiv bc\pmod p 且 (c,p)=1(c, p)=1 时，a≡b(modp)a\equiv b\pmod p
  8. a+b≡a mod p+b mod p(modp)a+b\equiv a\bmod p+b\bmod p\pmod p
  9. ab≡(a mod p)(b mod p)(modp)ab\equiv (a\bmod p)(b\bmod p)\pmod p
  10. a≡b(modp1)a\equiv b\pmod{p_1} 且 a≡b(modp2)a\equiv b\pmod{p_2} 时，a≡b(mod[p1,p2])a\equiv b\pmod{[p_1,p_2]}

简化剩余系

不会不知道用途在哪，爬了爬了[揩汗]

基本定理

  1. 欧拉定理：若 (a,p)=1(a,p)=1，一定有 aφ(p)≡1(modp)a^{\varphi(p)}\equiv 1\pmod p.
  2. 费马小定理：若 pp 为质数，则对于任意整数 aa，有 ap−1≡1(modp)a^{p-1}\equiv 1\pmod p.
  3. 扩展欧拉定理：若 b≥φ(p)b\ge \varphi(p)，则有 ab≡ab mod φ(p)+φ(p)(modp)a^b\equiv a^{b\bmod\varphi(p)+\varphi(p)}\pmod p
  4. 威尔逊定理：若 pp 为质数，则 (p−1)!≡p−1(modp)(p-1)!\equiv p-1\pmod p.

至于用途……降幂是主要的吧

比如 (222…) mod p(2^{2^{2^{\dots}}})\bmod p 显然满足扩展欧拉定理，那么令 f(a)=2f(a)f(a)=2^{f(a)}，则

2f(a)≡2f(a) mod φ(p)+φ(p)(modp)2^{f(a)}\equiv 2^{f(a)\bmod\varphi(p)+\varphi(p)}\pmod p

然后观察指数，发现 f(a) mod φ(p)f(a)\bmod\varphi(p) 实际上是 2f(a) mod φ(p){2}^{f(a)}\bmod\varphi(p)，那就一眼递归。然后 φ(p)\varphi(p) 是会一直缩小直到 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
// problem: 降幂大法[2]
// date: 2024.1.25-11:19
#include <bits/stdc++.h>
#define int long long
#define FH signed
using namespace std;
const int N = 1e7 + 10;
bool st[N];
int phi[N];
vector<int> p;
void getprime(int n)
{
    for(int i = 2; i <= n; i ++)
    {
        if(!st[i]) p.push_back(i), phi[i] = i - 1;
        for(auto v : p)
        {
            if(v * i > n) break;
            st[v * i] = 1;
            if(i % v == 0) {phi[i * v] = phi[i] * v; break;}
            else phi[i * v] = phi[i] * (v - 1);
        }
    }
}
int qmul(int a, int b, int p) {return (a * b - (int)((long double) a / p * b) * p + p) % p;}
int qpow(int a, int b, int p) {int res = 1; while(b) {if(b & 1) res = qmul(res, a, p); a = qmul(a, a, p); b >>= 1;} return res;}
int dfs(int p)
{
    if(p == 1) return 0;
    return qpow(2, dfs(phi[p]) + phi[p], p);
}
FH main()
{
    // freopen(".in", "r", stdin);
    // freopen(".out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    getprime(N - 1);
    int p; 
    while(cin >> p) cout << dfs(p) << "\n";
    return 0;
} //FH

扩展欧几里得算法

用于求得形如 ax+by=(a,b)ax+by=(a, b) 的方程的一组特解。

觉得这张图讲的特别好

https://s2.loli.net/2024/01/25/dQneFrGzcaubgm1.png

鸽鸽 寒假的以后找机会补吧。

其他

数论分块

对于常数 nn，使得式子

⌊ni⌋=⌊nj⌋\left\lfloor\frac ni\right\rfloor=\left\lfloor\dfrac nj\right\rfloor

成立的最大的满足 i≤j≤ni\leq j\leq n 的 jj 的值为 ⌊n⌊ni⌋⌋\left\lfloor\dfrac n{\lfloor\frac ni\rfloor}\right\rfloor。即值 ⌊ni⌋\left\lfloor\dfrac ni\right\rfloor 所在的块的右端点为 ⌊n⌊ni⌋⌋\left\lfloor\dfrac n{\lfloor\frac ni\rfloor}\right\rfloor。

同时可以证明，这样的块的数量总是在 n\sqrt n 级别的。

对于需要同时满足 ⌊a1i⌋\left\lfloor\dfrac {a_1}i\right\rfloor、⌊a2i⌋⋯⌊ani⌋\left\lfloor\dfrac {a_2}i\right\rfloor\cdots\left\lfloor\dfrac {a_n}i\right\rfloor 时，当前块的右端点为 min⁡j=1n{⌊aji⌋}\min\limits_{j=1}^n\{\left\lfloor\dfrac {a_j}i\right\rfloor\}。

莫比乌斯反演

出于这又臭又长的篇幅，请在这看->莫比乌斯反演学习笔记

拉格朗日插值

额……比较感性的算法，说真的，模板是小学生水平。

给出 nn 个点 (x1,y1),(x2,y2)…(xn,yn)(x_1,y_1),(x_2,y_2)\dots(x_n,y_n)，由这 nn 个点可以确定一个多项式 f(x)f(x)，求 f(k)f(k)。

显然我们可以很轻松的构造他：f(x)=∑i=1n∏j=1,i≠jnx−xjxi−xjyi\displaystyle f(x)=\sum_{i=1}^{n}\prod_{j=1,i\ne j}^{n}\frac{x-x_j}{x_i-x_j}y_i，这是显然的，n2n^2 解决就好了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void exgcd(int a, int b, int &x, int &y) {if(b == 0) {x = 1, y = 0; return ;} 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 gety(int x[], int y[], int n, int k, int mod)
{
    int res = 0;
    for(int i = 1; i <= n; i ++)
    {
        int a = 1, b = 1;
        for(int j = 1; j <= n; j ++) if(i != j)
            a = a * (((k - x[j]) % mod + mod) % mod) % mod,
            b = b * (((x[i] - x[j]) % mod + mod) % mod) % mod;
        res = (res + y[i] * a % mod * inv(b, mod) % mod) % mod;
    }
    return res;
}
相关文章
生成函数初步在车人的要求下，玉帛拿一下午讲了生成函数 这个大坑不知道要填多久 前置知识：求导（真够让人头大） 这是讲解：懒得写了...莫比乌斯反演学习笔记莫比乌斯反演 一句话总结：题面极其简洁并且大多能一眼出算法，代码短但是式子又臭又长。 前置知识：数论分块：数论全家桶 -...2-SAT下文中以 ⟨⟩\langle\rangle⟨⟩ 表示矛盾的集合，[][][] 表示题目给定的的集合 定义 2-SAT，简...20230920模拟赛挂分寄又挂分了……KMP与Z函数0x00 引入 给定主串与模式串 ，问主串中是否存在模式串。 实际上这是一个最为简单的字符串匹配问题，对于暴力匹配来讲，...

更新于：2024年2月23日

oi

数学

评论