数论基础总结

前言

1. 今天【2021.11.13】总结之前学习的博客，才发现，只有大二下，才系统的学习了很多东西，其他时间虽然也是在学习，但是效率实在不敢恭维。
   1. 什么概念？大二下，一周大概能把【网络流】或者【数论基础】或者【最短路（全）】或者【最小生成树（全）】全给搞定，在之前的话，怕是一个月也不一定学到那么多，那么系统（特别网络流还真得下点决心去学，很烧脑，而且学的还是比较深入了，虽然后面没怎么用到（二分图匹配用到过））。——效率之上的学习方式，一周如果能有平均每天6小时的高效学习时间，那么这一周便不算是虚度。

参考博客

1. 【已更新】《算法竞赛中的初等数论》（ACM / OI / MO）前言、后记、目录索引（十五万字符的数论书）
2. oi-wiki数学部分
3. 神O的数论全家桶

专题训练

1. 传送门：“kuangbin带你飞”专题计划——专题十四：数论基础

素数筛

埃氏筛（拓展性很强）

int p[N], vis[N];
void Prime(int n) {
    int i, j;  //快一点点
    for (i = 2; i <= n; i++) {
        if (vis[i]) continue;
        p[++p[0]] = i;
        for (j = 2 * i; j <= n; j += i) vis[j] = 1;
    }
}

欧拉筛（线性筛）

int p[N], vis[N];
void Prime(int n) {
    int i, j;
    for (i = 2; i <= n; i++) {
        if (!vis[i]) p[++p[0]] = i;
        //注意，进入循环的并非都是素数
        for (j = 1; j <= p[0]; j++) {
            if (i * p[j] > n) break;
            vis[i * p[j]] = 1;
            if (i % p[j] == 0) break;
            //防止重复。比如4*2=8,4*3=12，12还是留在6*2的时候筛
        }
    }
}

欧拉函数

欧拉函数的定义

​​

欧拉函数的一些性质

欧拉函数模板（单个&线性）

//由性质四求某个数的欧拉函数
int euler_phi(int n) {
    int ans = n;
    for (int i = 2; i * i <= n; i++)
        if (n % i == 0) {
            ans = ans / i * (i - 1);
            while (n % i == 0) n /= i;
        }
    if (n > 1) ans = ans / n * (n - 1);
    return ans;
}

//由性质四，线性筛求欧拉函数。时间复杂度和埃氏筛一样
void phi_table(int n, int* phi) {
    for (int i = 2; i <= n; i++) phi[i] = 0;
    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        //素数才进行操作
        if (!phi[i]) {
            for (int j = i; j <= n; j += i) {
                if (!phi[j]) phi[j] = j;  //还没操作的时候另phi[j]=j
                phi[j] = phi[j] / i * (i - 1);  //会遇到j的每一个素数
            }
        }
    }
}

欧拉定理

扩展欧拉定理

1. 例题：Bi-shoe and Phi-shoe LightOJ - 1370 (线性筛+欧拉函数)

欧拉常数（貌似和欧拉函数没多大关系来着？？）

1. 题目：Harmonic Number LightOJ - 1234（调和级数+欧拉常数+打表）
2. 题意：$T\le 10000$组样例，每组样例一个数$1\le n\le 1e8$，求，注意全’/‘不是取整数部分，而是当小数运算。
3. 题解：
   1. 方法一：打表，每隔$100$个数记录一次答案，然后枚举。
   2. 方法二：公式法$H_n=\ln n+\frac{1}{2n}+\gamma$越大误差越小，当n较小的时候建议枚举。其中欧拉常数$\gamma =0.57721566490153286060651209$。
4. 补充：log(1.0*n)，其中log的底默认为自然对数！

//方法一
#include <bits/stdc++.h>
// #define int long long
#define read(x) scanf("%d", &x)
#define print(a, c) printf("%d%c", a, c)
#define dbg(x) cout << #x << "===" << x << endl
#define pb push_back
using namespace std;
const int N = 1e6 + 20;
 
int n;
double dp[N];
void init() {
    int i;
    double x = 0.0;
    dp[0] = x;
    for (i = 1; i <= 1e8; i++) {
        x += 1.0 / i;
        if (i % 100 == 0) dp[i / 100] = x;
    }
    // for (i = 0; i <= 10; i++) printf("%.10lf\n", dp[i]);
}
signed main() {
    int T;
    read(T);
    init();
    for (int _ = 1; _ <= T; _++) {
        read(n);
        double ans = dp[n / 100];
        int i;
        for (i = n - n % 100 + 1; i <= n; i++) ans += 1.0 / i;
        printf("Case %lld: %.10lf\n", _, ans);
    }
    return 0;
}

唯一分解定理

1. 参考博客：唯一分解定理一篇就够了
2. 定义：
   

应用1：求所有因数的个数

1. 例题：Aladdin and the Flying Carpet LightOJ - 1341 （唯一分解定理+线筛+一个数的因数个数）

应用2：求所有因数之和

​​1. 例题：Sigma Function LightOJ - 1336 （Sigma函数+一个数的所有因子和）
2. 例题传送门：Sigma Function LightOJ - 1336

#include <bits/stdc++.h>
#define int long long
#define pb push_back
#define dbg(x) cout << #x << "===" << x << endl
using namespace std;
const int N = 1e6 + 10;
 
int n;
vector<int> v;
void init() {
    v.pb(0);
    for (int i = 1; i <= 1e6; i++) {
        v.pb(i * i), v.pb(2 * i * i);
    }                          //最大2e12
    sort(v.begin(), v.end());  //别忘了
    // for (int i = 0; i <= 10; i++) cout << i << " " << v[i] << endl;
}
signed main() {
    int T;
    scanf("%lld", &T);
    init();
    for (int _ = 1; _ <= T; _++) {
        scanf("%lld", &n);
        int ans = n;
        int pos = upper_bound(v.begin(), v.end(), n) - v.begin();
        --pos;
        // dbg(pos);
        ans -= pos;
        printf("Case %lld: %lld\n", _, ans);
    }
    return 0;
}
/*
input:::
4
3
10
100
1000
output:::
Case 1: 1
Case 2: 5
Case 3: 83
Case 4: 947
*/

应用3：求LCM/GCD

1. 例题：Pairs Forming LCM LightOJ - 1236 （唯一分解定理的应用+线筛+MLE的情况）

扩展欧几里得算法

模板：拓展欧几里得&拓展欧几里得求逆元

int extend_gcd(int a, int b, int &x, int &y) {
    if (a == 0 && b == 0) return -1;  //无gcd。一般不会出现这种情况
    //求逆元的时候需要用到的
    if (b == 0) {
        x = 1, y = 0;
        return a;  // gcd边界
    }
    int g = extend_gcd(b, a % b, y, x);
    y -= a / b * x;  //这两行，真不懂。不管了
    return g;
}
int inv(int a, int M) {
    int X0, Y0, X, Y;
    int g = extend_gcd(a, M, X0, Y0);
    // gcd=1，所以：mod_x=b,mod_y=a
    return (g == 1) ? (X0 % M + M) % M : -1;  // c=1,g=1,X=c/g*X0=1,mod_x=b/1=M
}

洛谷P1516 青蛙的约会（数论同余、扩欧）（现在我们只需要关注ax+by+c=0怎么解得x，y即可！！！）

#include <bits/stdc++.h>
#define int long long
 
using namespace std;
int x, y, m, n, l;
int a, b, c, g, X0, Y0, X, Y;
int extend_gcd(int a, int b, int &x, int &y) {
    if (a == 0 && b == 0) return -1;  //无gcd。一般不会出现这种情况
    if (b == 0) {
        x = 1, y = 0;
        return a;  // gcd边界
    }
    int g = extend_gcd(b, a % b, y, x);
    y -= a / b * x;  //这两行，真不懂。不管了
    return g;
}
signed main() {
    cin >> x >> y >> m >> n >> l;
    a = m - n, b = l, c = y - x;  // c别搞错了
    if (a < 0) a = -a, c = -c;
    g = extend_gcd(a, b, X0, Y0);  // g已经得到，不必再赋值
    int mod_x = b / g, mod_y = a / g;
    X = ((c * X0 / g) % mod_x + mod_x) % mod_x;
    if (c % g)
        puts("Impossible");
    else
        cout << X << endl;
    return 0;
}

乘法逆元（3种方法）

！！！一个数a与模数mod互质的时候才存在逆元

参考博客：逆元的求法总结（3种基本方法+4种实现）

方法1：扩展欧几里得算法

1. 时间复杂度：$O(\log n)$
2. mod可以不为质数
   

int extend_gcd(int a, int b, int &x, int &y) {
    if (a == 0 && b == 0) return -1;  //无gcd。一般不会出现这种情况
    //求逆元的时候需要用到的
    if (b == 0) {
        x = 1, y = 0;
        return a;  // gcd边界
    }
    int g = extend_gcd(b, a % b, y, x);
    y -= a / b * x;  //这两行，真不懂。不管了
    return g;
}
int inv(int a, int M) {
    int X0, Y0, X, Y;
    int g = extend_gcd(a, M, X0, Y0);
    // gcd=1，所以：mod_x=b,mod_y=a
    return (g == 1) ? (X0 % M + M) % M : -1;  // c=1,g=1,X=c/g*X0=1,mod_x=b/1=M
}

方法2：费马小定理/欧拉定理

费马小定理：inv(a)=a^(mod-2)%mod，mod只能为质数。
欧拉定理：inv(a)=a^(phi(mod)-1)%mod，mod可以不为质数，只要互质（即有逆元/求出来的逆元有意义）就行。

方法3：递推求逆元

1. mod只能是质数。
2. 时间复杂度$O(\log m)$。

//递归版本
int inv(int x) { return (x == 1) ? 1 : (mod - mod / x) * inv(mod % x) % mod; }
//数组版本（线性）
inv[1] = 1;
for (int i = 2; i <= n; i++) inv[i] = (mod - mod / i) * inv[mod % i] % mod;

中国剩余定理

1. 详细讲解：oi-wiki_中国剩余定理
2. 情形说明：
   
3. 算法操作：
   

模板：中国剩余定理

P3868 [TJOI2009]猜数字

#include <bits/stdc++.h>
#define int long long
#define read(x) scanf("%lld", &x)
#define print(a, c) printf("%lld%c", a, c)
#define dbg(x) cout << #x << "===" << x << endl
using namespace std;
const int N = 20;

int n, a[N], b[N];

int exgcd(int a, int b, int &x, int &y) {
    if (a == 0 && b == 0) return -1;  //好像没啥卵用
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}
//因为任意两个b[i]互质，所以lcm/b[i]在b[i]意义下的逆元一定存在，因为gcd(lcm/b[i],b[i])=1
int inv(int a, int M) {
    int X0, Y0;
    int g = exgcd(a, M, X0, Y0);
    return (g == 1) ? (X0 % M + M) % M : -1;  //-1不存在逆元
}
//龟速乘，以下所有值不会超过2*lcm，保证了不会超long long
int qmul(int a, int b, int mod) {
    int res = 0;  //注意这里也变成了0
    while (b) {
        if (b & 1) res = (res + a) % mod;
        b >>= 1, a = 2 * a % mod;
    }
    return res;
}
int china() {
    int ans = 0, lcm = 1;
    int M, c;
    for (int i = 1; i <= n; i++) lcm *= b[i];
    for (int i = 1; i <= n; i++) {
        M = lcm / b[i];
        c = M * inv(M, b[i]) % lcm;  // inv(M)是对bi取余，a[i]*c[i]都是对lcm取余
        // ans = (ans + c * a[i] % lcm) % lcm;
        //以上操作。最大可能出现的数达到lcm^2数量级会爆long long，需要使用快速乘
        // inv中的数不会超过max(bi)
        ans = (ans + qmul(c, a[i], lcm)) % lcm;
        // qmul中c在前会快一点，因为一半ai<bi小一点
    }
    return (ans + lcm) % lcm;
}
signed main() {
    read(n);
    int lcm = 1;
    for (int i = 1; i <= n; i++) read(a[i]);
    for (int i = 1; i <= n; i++) read(b[i]);
    for (int i = 1; i <= n; i++) a[i] = (a[i] % b[i] + b[i]) % b[i];
    //这是为何。使符合中国剩余定理的条件ai=N%bi，则ai<bi
    // bi|(N-ai)，看来意思是(N-ai)是bi的整数倍，不管这个整数是多少（<=0也oK）
    print(china(), '\n');
    return 0;
}

解决模指数方程

//[Emoogle Grid UVA - 11916 （BSGS算法，a^x+by=c求x）](https://blog.csdn.net/I_have_a_world/article/details/117330991)
//求满足a^x+p*y=b的最小非负数x。返回-1表示没有答案
int BSGS(int a, int b, int p) {
    //答案x=i*m+j.分成一块正方形，复杂度最高sqrt(p)*log(p)就可以解出x，如果解不出，只能说无解（具有周期性，且余数个数不会超过phi[p]）
    int m = (int)ceil(sqrt(1.0 * p));
    map<int, int> mp;
    int t = b, i, j, inv_a = qpow(a, p - 2, p),
        a_m = qpow(a, m, p);  // a_i表示的是a^m
    mp[b] = 0;                //每一行0~m-1
    for (j = 1; j < m; j++) {
        t = t * inv_a % p;
        if (!mp[t]) mp[t] = j;  //模数相等，只需要最小的i
    }
    t = 1;
    for (i = 0; i < m; i++) {
        if (mp[t]) return i * m + mp[t];
        t = t * a_m % p;
    }
    return -1;  //无解
}

推荐拓展学习：

夜深人静写算法（三）- 初等数论入门

莫比乌斯反演

参考资料

1. peng-ym-莫比乌斯反演

预备知识：整数分块（感觉就是一个思想）

1. 参考博客：peng-ym-整除分块
   

总结

1. 一般形式：${\sum }_{i=1}^n[n/i]$
2. 过程：r=n/(n/l)，然后不断枚举l，求得r，[l,r]部分的值是一样的，不必遍历，直接得到长度*数。
3. 复杂度：由$O(n)$变为$O(\sqrt{n})$

莫比乌斯函数$\mu$（数论函数也可以视作一个数列）

总结

1. 定义：很重要
2. 性质1：很重要
3. 性质2还需要回顾以下欧拉函数才行。

莫比乌斯反演

题目训练

题目1（oi-wiki推荐）：求满足gcd(x,y)=k且x,y分别在区间[a,b],[c,d]内的个数（莫比乌斯反演&容斥定理&整数分块）

1. 题目：P2522 [HAOI2011]Problem b
   
   

2. 题解：

   1. 首先由容斥定理把题目简化为求四个${\sum }_{i=1}^n{\sum }_{j=1}^m[gcd(i,j)=1]$的形式。
   2. 其中$[gcd(i,j)=1]$表示满足内部条件时返回true，否则返回false。
   3. 由莫比乌斯函数性质1：${\sum }_{d|n}\mu (d)=[n=1]$
   4. 带入得到：${\sum }_{i=1}^n{\sum }_{j=1}^m{\sum }_{d|gcd(i,j)}\mu (d)$（即条件转化为具体求和式子）
   5. 变形得到：${\sum }_{d=1}^{min(n,m)}{\sum }_{i=1}^n[d|i]{\sum }_{j=1}^m[d|j]$。
   6. 最后用整数分块解决以下就ok了（不过这里怎么就敢说复杂度为$O(\sqrt{max(n,m})$呢？还是雀食不是根号，但是也大不了多少）。

3. 代码：

/*
4. 莫比乌斯反演
*/
#include <bits/stdc++.h>
using namespace std;
const int N = 5e4 + 5;
int a, b, c, d, k;
int mu[N], p[N], vis[N];

void Moblus(int n) {
    memset(vis, 0, sizeof(vis));
    mu[1] = 1;
    int tot = 0;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) p[tot++] = i, mu[i] = -1;
        for (int j = 0; j < tot; j++) {
            if (i * p[j] > n) break;
            vis[i * p[j]] = 1;
            if (i % p[j] == 0) {
                mu[i * p[j]] = 0;
                break;
            } else
                mu[i * p[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= n; i++) mu[i] += mu[i - 1];
}
int solve(int n, int m) {
    int res = 0;
    for (int i = 1, j; i <= min(n, m); i = j + 1) {
        j = min(
            n / (n / i),
            m / (m / i));  // emmm，不敢多问，这你敢说是$O(\sqrt{n})$复杂度？
            // 得到一个正确的限制条件，然后模式化（无脑）枚举就好emm
        res += (mu[j] - mu[i - 1]) * (n / i) * (m / i);
    }
    return res;
}
//容斥定理
int f(int x, int k) {
    if (x % k == 0) return x / k - 1;
    return x / k;
}
signed main() {
    Moblus(N - 1);
    int T;
    cin >> T;
    while (T--) {
        scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
        int ans = solve(b / k, d / k) - solve(b / k, f(c, k)) -
                  solve(f(a, k), d / k) + solve(f(a, k), f(c, k));
        printf("%d\n", ans);
    }
    return 0;
}

题目2（cf2200）：

NTT（快速数论变换）

1. 数论变换是一种计算卷积（convolution）的快速算法。