数论模板
|
521
|
0
|
模板

1273 字
|
57 分钟

数论常用模板和一些经典问题

简单的防爆模板

namespace SimpleModInt {
    constexpr int md = (int)1e9 + 7;
    inline int norm(long long a) {
        return (a % md + md) % md;
    }
    inline int add(int a, int b) {
        a += b;
        if (a >= md) a -= md;
        return a;
    }
    inline int sub(int a, int b) {
        a -= b;
        if (a < 0) a += md;
        return a;
    }
    inline int mul(int a, int b) {
        return (int)((long long)a * b % md);
    }
    inline int powmod(int a, long long b) {
        int res = 1;
        while (b > 0) {
            if (b & 1) res = mul(res, a);
            a = mul(a, a);
            b >>= 1;
        }
        return res;
    }
    inline int inv(int a) {
        a %= md;
        if (a < 0) a += md;
        int b = md, u = 0, v = 1;
        while (a) {
            int t = b / a;
            b -= t * a; swap(a, b);
            u -= t * v; swap(u, v);
        }
        assert(b == 1);
        if (u < 0) u += md;
        return u;
    }
}
using namespace SimpleModInt;

素数

求约数个数

int getFactors(long long n) { // 求 n 的因子个数
    if (n < 2) return 1; // 注意一个优化 我们可以先筛出O(sqrt(border)) 的质因子再求因子个数
    int ans = 1;
    for (long long i = 2; i * i <= n; ++i) {
        int cnt = 1;
        while (!(n % i)) {
            ++cnt;
            n /= i;
        }
        ans *= cnt;
    }
    return (n == 1 ? ans : ans + ans);
}

求约数和

// 求因子和 n必须为正整数 模数可自定义
long long getSigma(long long n, long long mod = 1e18) {
    if (n == 1) return 1;
    long long ans = 1;
    for (long long i = 2; i * i <= n; ++i) {
        long long cnt = 1;
        while (!(n % i)) {
            cnt = (cnt * i + 1) % mod;
            n /= i;
        }
        ans *= cnt;
    }
    return (n == 1 ? ans : ans * (n + 1) % mod);
}

Pollard-Rho && Miller-Rabin

namespace Pollard_Rho {
    typedef long long ll;
    vector<ll> ans; // 存储质因子的数组
    inline ll gcd(ll a, ll b) { ll c; while (b) c = a % b, a = b, b = c; return a; }
    inline ll mulmod(ll x, ll y, const ll z) {
        return (x * y - (ll)(((long double)x * y + 0.5) / (long double)z) * z + z) % z;
    }
    inline ll powmod(ll a, ll b, const ll mo) {
        ll s = 1;
        for (; b; b >>= 1, a = mulmod(a, a, mo)) if (b & 1) s = mulmod(s, a, mo);
        return s;
    }
    bool isPrime(ll p) { // Miller-Rabin O(klog^3(n)) k为素性测试轮数
        const int lena = 10, a[lena] = { 2,3,5,7,11,13,17,19,23,29 };
        if (p == 2) return true;
        if (p == 1 || !(p & 1) || (p == 46856248255981ll)) return false;
        ll D = p - 1; while (!(D & 1)) D >>= 1;
        for (int i = 0; i < lena && a[i] < p; i++) {
            ll d = D, t = powmod(a[i], d, p); if (t == 1) continue;
            for (; d != p - 1 && t != p - 1; d <<= 1) t = mulmod(t, t, p);
            if (d == p - 1) return false;
        }
        return true;
    }
    void reportFactor(ll n) { // 得到一个素因子
        ans.emplace_back(n); // 存储素因子
    }
    ll ran() { return rand(); } // 随机数
    void getFactor(ll n) { // Pollard-Rho O(n ^ 1/4)
        if (n == 1) return;
        if (isPrime(n)) { reportFactor(n); return; }
        while (true) {
            ll c = ran() % n, i = 1, x = ran() % n, y = x, k = 2;
            do {
                ll d = gcd(n + y - x, n);
                if (d != 1 && d != n) { getFactor(d); getFactor(n / d); return; }
                if (++i == k) y = x, k <<= 1;
                x = (mulmod(x, x, n) + c) % n;
            } while (y != x);
        }
    }
}
using namespace Pollard_Rho;

逆元

vector<long long> inv;
void linearInv(int n, long long p) { // 线性求逆元
    inv.resize(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) {
        inv[i] = (p - p / i) * inv[p % i] % p;
    }
}

扩展欧几里得解线性方程

exgcd求 $ax+by=c$ 的最小非负整数解详解：

  1. 求出 $a,b$ 的最大公约数 $g=\gcd(a,b)$ ，根据裴蜀定理检查是否满足 $c\% g=0$ ，不满足则无解；
  2. 调整系数 $a,b,c$ 为 $a’=\frac{a}{g},b’=\frac{b}{g},c’=\frac{c}{g}$ ，这是因为 $ax+by=c$ 和 $a’x+b’y=c’$ 是完全等价的；
  3. 实际上exgcd求解的方程是 $a’x+b’y=1$ ，求解前需要注意让系数 $a’,b’\geq 0$ （举个例子，如果系数 $b’$ 原本 $<0$ ，我们可以翻转 $b’$ 的符号然后令解 $(x,y)$ 为 $(x,-y)$ ，但是求解的时候要把 $y$ 翻回来）；
  4. 我们通过exgcd求出一组解 $(x_0,y_0)$ ，这组解满足 $a’x_0+b’y_0=1$ ，为了使解合法我们需要令 $x_0=c’x_0,y_0=c’y_0$ ，于是有 $a'(c’x_0)+b'(c’y_0)=c”$ ；
  5. 考虑到 $a’x_0+b’y_0=1$ 等价于同余方程 $a’x_0\equiv 1\pmod{b’}$ ，因此为了求出最小非负整数解，我们最后还需要对 $b’$ 取模，
  6. 最后注意特判 $c’=0$ 的情况，如果要求解 $y$ 且系数 $b$ 发生了翻转，将其翻转回来。
long long exgcd(long long a, long long b, long long& x, long long& y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    long long g = exgcd(b, a % b, y, x);
    y -= (a / b) * x;
    return g;
}

ll x, y; // 最小非负整数解
bool solve(ll a, ll b, ll c) { // ax+by=c
    ll g = gcd(a, b);
    if (c % g) return false;
    a /= g, b /= g, c /= g;
    bool flag = false;
    if (b < 0) b = -b, flag = true;
    exgcd(a, b, x, y);
    x = (x * c % b + b) % b;
    if (flag) b = -b;
    y = (c - a * x) / b;
    if (!c) x = y = 0; // ax+by=0
    return true;
}

一定范围内线性方程整数解的个数

exgcd通解：假设我们通过上方的exgcd流程获得了一组解 $(x_0,y_0)$ （没有乘 $c$ ），那么 $a’x+b’y=1$ 的通解就是 $(x_0+b’t,y_0-a’t)$ ，因此 $a’x+b’y=c’$ 的通解是 $(c'(x_0+b’t),c'(y_0-a’t))$ 。

/*
 * 求解 ax+by+c=0 模板CF710D
 * 返回 [xl, xr] x [yl, yr] 的解数
 * 若至少有一组解 则x是一个合法解
 * 可以根据x推出y 但要注意b=0/a=0等特殊情况
 * 特别注意！！！设置边界时的取整问题！！！
 * x/y 向上取整时(x+y-1)/y 向下取整时floor(1.0*x/y)
 * */
ll a, b, c, xl, xr, yl, yr;
ll x, y, d;
ll exgcd(ll a, ll b, ll& x, ll& y) {
    if (!b) return x = 1, y = 0, a;
    ll d = exgcd(b, a % b, x, y), t = x;
    return x = y, y = t - a / b * y, d;
}
ll solve(ll a, ll b, ll c, ll xl, ll xr, ll yl, ll yr) {
    if (xl > xr) return 0;
    if (yl > yr) return 0;
    if (!a && !b) {
        if (c) return 0;
        return (xr - xl + 1) * (yr - yl + 1);
    }
    if (!b) {
        swap(a, b);
        swap(xl, yl);
        swap(xr, yr);
    }
    if (!a) {
        if (c % b) return 0;
        ll y = -c / b;
        if (y < yl || y > yr) return 0;
        return xr - xl + 1;
    }
    d = exgcd((a % abs(b) + abs(b)) % abs(b), abs(b), x, y);
    if (c % d) return 0;
    x = (x % abs(b) + abs(b)) % abs(b) * ((((-c) % abs(b)) + abs(b)) % abs(b) / d) % abs(b / d);
    d = abs(b / d);
    ll kl = (xl - x) / d - 3, kr = (xr - x) / d + 3;
    while (x + kl * d < xl) kl++;
    while (x + kr * d > xr) kr--;
    ll A = (-yl * b - a * x - c) / (a * d), B = (-yr * b - a * x - c) / (a * d);
    if (A > B) swap(A, B);
    kl = max(kl, A - 3);
    kr = min(kr, B + 3);
    while (kl <= kr) {
        ll y = (-c - a * x - a * d * kl) / b;
        if (yl <= y && y <= yr) break;
        kl++;
    }
    while (kl <= kr) {
        ll y = (-c - a * x - a * d * kr) / b;
        if (yl <= y && y <= yr) break;
        kr--;
    }
    if (kl > kr) return 0;
    return kr - kl + 1;
}

筛法

线性素数筛

vector<bool> isPrime; // true 表示非素数  false 表示是素数
vector<int> prime; // 保存素数
int sieve(int n) {
    isPrime.resize(n + 1, false);
    isPrime[0] = isPrime[1] = true;
    for (int i = 2; i <= n; i++) {
        if (!isPrime[i]) prime.emplace_back(i);
        for (int j = 0; j < (int)prime.size() && prime[j] * i <= n; ++j) {
            isPrime[prime[j] * i] = true;
            if (!(i % prime[j])) break;
        }
    }
    return (int)prime.size();
}

线性欧拉函数筛

bool is_prime[SIZE];
int prime[SIZE], phi[SIZE]; // phi[i] 表示 i 的欧拉函数值
int Phi(int n) { // 线性筛素数的同时线性求欧拉函数
    phi[1] = 1; is_prime[1] = true;
    int p = 0;
    for (int i = 2; i <= n; ++i) {
        if (!is_prime[i]) prime[p++] = i, phi[i] = i - 1;
        for (int j = 0; j < p && prime[j] * i <= n; ++j) {
            is_prime[prime[j] * i] = true;
            if (!(i % prime[j])) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
    return p;
}

线性约数个数筛

bool is_prime[SIZE];
int prime[SIZE], d[SIZE], num[SIZE]; // d[i] 表示 i 的因子数  num[i] 表示 i 的最小质因子出现次数
int getFactors(int n) { // 线性筛因子数
    d[1] = 1; is_prime[1] = true;
    int p = 0;
    for (int i = 2; i <= n; ++i) {
        if (!is_prime[i]) prime[p++] = i, d[i] = 2, num[i] = 1;
        for (int j = 0; j < p && prime[j] * i <= n; ++j) {
            is_prime[prime[j] * i] = true;
            if (!(i % prime[j])) {
                num[i * prime[j]] = num[i] + 1;
                d[i * prime[j]] = d[i] / num[i * prime[j]] * (num[i * prime[j]] + 1);
                break;
            }
            num[i * prime[j]] = 1;
            d[i * prime[j]] = d[i] + d[i];
        }
    }
    return p;
}

线性质因子个数筛

bool is_prime[SIZE];
int prime[SIZE], num[SIZE]; // num[i] 表示 i 的质因子数
int getPrimeFactors(int n) { // 线性筛质因子数
    is_prime[1] = true;
    int p = 0;
    for (int i = 2; i <= n; ++i) {
        if (!is_prime[i]) prime[p++] = i, num[i] = 1;
        for (int j = 0; j < p && prime[j] * i <= n; ++j) {
            is_prime[prime[j] * i] = true;
            if (!(i % prime[j])) {
                num[i * prime[j]] = num[i];
                break;
            }
            num[i * prime[j]] = num[i] + 1;
        }
    }
    return p;
}

线性约数和函数筛

bool is_prime[SIZE];
int prime[SIZE], f[SIZE], g[SIZE]; // f[i] 表示 i 的约数和
int getSigma(int n) {
    g[1] = f[1] = 1; is_prime[1] = true;
    int p = 0;
    for (int i = 2; i <= n; ++i) {
        if (!is_prime[i]) prime[p++] = i, f[i] = g[i] = i + 1;
        for (int j = 0; j < p && prime[j] * i <= n; ++j) {
            is_prime[prime[j] * i] = true;
            if (!(i % prime[j])) {
                g[i * prime[j]] = g[i] * prime[j] + 1;
                f[i * prime[j]] = f[i] / g[i] * g[i * prime[j]];
                break;
            }
            f[i * prime[j]] = f[i] * f[prime[j]];
            g[i * prime[j]] = 1 + prime[j];
        }
    }
    return p;
}

线性莫比乌斯函数筛

bool is_prime[SIZE];
int prime[SIZE], mu[SIZE]; // mu[i] 表示 i 的莫比乌斯函数值
int getMu(int n) { // 线性筛莫比乌斯函数
    mu[1] = 1; is_prime[1] = true;
    int p = 0;
    for (int i = 2; i <= n; ++i) {
        if (!is_prime[i]) prime[p++] = i, mu[i] = -1;
        for (int j = 0; j < p && prime[j] * i <= n; ++j) {
            is_prime[prime[j] * i] = true;
            if (!(i % prime[j])) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    return p;
}

欧拉函数

欧拉定理

a^b\equiv
\begin{cases}
a^{b\bmod\varphi(p)},\,&\gcd(a,\,p)=1\\
a^b,&\gcd(a,\,p)\ne1,\,b<\varphi(p)\\
a^{b\bmod\varphi(p)+\varphi(p)},&\gcd(a,\,p)\ne1,\,b\ge\varphi(p)
\end{cases}
\pmod p

$O(\sqrt{n})$ 求欧拉函数

long long phi(long long n) {
    int m = int(sqrt(n + 0.5));
    long long ans = n;
    for (int i = 2; i <= m; 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;
}

$O(\sqrt{n})$ 预处理， $O(\log{n})$ 求欧拉函数

vector<int> prime;  // 求 n 的欧拉函数需要先把 <= ceil(sqrt(n)) 的素数筛出
long long phi(long long n) {
    long long res = n;
    for (int i = 0; i < (int)prime.size(); i++) {
        long long p = prime[i];
        if (p * p > n) break;
        if (n % p == 0) {
            res = res / p * (p - 1);
            while (n % p == 0) n /= p;
        }
    }
    if (n > 1) res = res / n * (n - 1);
    return res;
}

BSGS

long long qpow(long long a, long long b, long long mod) { // 快速幂 a ^ b % mod
    long long ans = 1;
    while (b) {
        if (b & 1) ans = (ans * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return ans;
}

long long BSGS(long long a, long long b, long long m) { // a ^ x = b mod m
    long long n = (long long)sqrt((double)m) + 1;
    map<int, int> MP; // 可以换 unordered_map
    b %= m;
    for (int i = 0; i < n; ++i) MP[b * qpow(a, i, m) % m] = i;
    a = qpow(a, n, m);
    if (!a) return b ? -1 : 1;
    for (int i = 0; i <= n; ++i) {
        long long val = qpow(a, i, m);
        int j = (!MP.count(val) ? -1 : MP[val]);
        if (j >= 0 && i * n >= j) return i * n - j;
    }
    return -1; // 无解
}

中国剩余定理

CRT

// 求解形如 x = ci (mod mi) 的线性方程组 mi 必须是素数
long long CRT(vector<long long>& c, vector<long long>& m) {
    long long M = m[0], ans = 0;
    for (int i = 1; i < (int)m.size(); ++i) M *= m[i];
    for (int i = 0; i < (int)m.size(); ++i) { // Mi * ti * ci
        long long mi = M / m[i];
        long long ti = inv(mi, m[i]); // mi 模 m[i] 的逆元
        ans = (ans + mi * ti % M * c[i] % M) % M;
    }
    ans = (ans + M) % M; // 返回模 M 意义下的唯一解
    return ans;
}

EXCRT

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

long long mulmod(long long x, long long y, const long long z) { // x * y % z 防爆
    return (x * y - (long long)(((long double)x * y + 0.5) / (long double)z) * z + z) % z;
}

// 求解形如 x = ci (mod mi) 的线性方程组
long long EXCRT(vector<long long>& c, vector<long long>& m) {
    long long M = m[0], ans = c[0];
    for (int i = 1; i < (int)m.size(); ++i) { // M * x - mi * y = ci - C
        long long x, y, C = ((c[i] - ans) % m[i] + m[i]) % m[i]; // ci - C
        long long G = exgcd(M, m[i], x, y);
        if (C % G) return -1; // 无解
        long long P = m[i] / G;
        x = mulmod(C / G, x, P); // 防爆求最小正整数解 x
        ans += x * M;
        M *= P; // LCM(M, mi)
        ans = (ans % M + M) % M;
    }
    return ans;
}

迪利克雷卷积

$$
\begin{aligned}
(f*g)(n)=\underset{d|n}{\sum} f(d)g(\frac{n}{d})=\underset{xy=n}{\sum} f(x)g(y)
\end{aligned}
$$

杜教筛

$$
\begin{aligned}
g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S(\lfloor\frac{n}{i}\rfloor)
\end{aligned}
$$

用于计算数论函数 $S(n)=\sum_{i=1}^n f(i)$ 的前缀和。以 $\mu$ 为例：

unordered_map<ll, ll> mapMu;
ll calMu(ll n) { // 计算 sigma(mu(i)) i=1,2,...,n
    if (n <= border) { // 需要先筛出前 n^2/3 项mu[i]复杂度才是 O(n^2/3) 筛其他函数同理
        return mu[n];  // border 设置为 n^2/3
    } else if (mapMu.count(n)) {
        return mapMu[n];
    } else {
        ll res = 0;
        for (ll l = 2, r; l <= n; l = r + 1) {
            r = n / (n / l);
            res += calMu(n / l) * (r - l + 1);
        }
        return mapMu[n] = 1 - res;
    }
}

$x^2+y^2=n$ 整数解数/圆环整点数

正整数解数量：将 $n$ 分解为 $2^x\prod p_i^{c_i}\prod q_i^{d_i}$ ，其中 $p_i$ 为形如 $4k+1$ 的质因子，$q_i$ 为 $4k+3$ 的质因子，解的总数为 $4\prod (c_i+1)$ 。

蔡勒公式

int zeller(int y, int m, int d) { // 蔡勒公式 返回星期几
    if (m <= 2) y--, m += 12;
    int c = y / 100; y %= 100;
    int w = ((c >> 2) - (c << 1) + y + (y >> 2) + 
        (13 * (m + 1) / 5) + d - 1) % 7;
    if (w < 0) w += 7;
    return (w);
}
int getId(int y, int m, int d) { // 返回到公元1年1月1日的天数
    if (m < 3) { y--; m += 12; }
    return 365 * y + y / 4 - y / 100 + y / 400 +
        (153 * (m - 3) + 2) / 5 + d - 307;
}

分数类

template<typename T> struct fraction {
    T numerator;   // 分子
    T denominator; // 分母
    T gcd(T a, T b) { return b == 0 ? a : gcd(b, a % b); }
    constexpr fraction(int numerator_ = 0, int denominator_ = 1)
        : numerator(numerator_), denominator(denominator_) {}
    constexpr fraction(const fraction& k) : numerator(k.numerator), denominator(k.denominator) {}
    constexpr fraction norm() { // 最简化
        T g = gcd(numerator, denominator);
        if (g < 0) g = -g;
        if (denominator < 0) {  // 如果是负数 保证分子<0 分母>0
            return fraction(-numerator / g, -denominator / g);
        } else {
            return fraction(numerator / g, denominator / g);
        }
    }
    constexpr fraction operator+ (const fraction& k) const noexcept {
        return fraction(numerator * k.denominator + denominator * k.numerator, denominator * k.denominator).norm();
    }
    constexpr fraction operator- (const fraction& k) const noexcept {
        return fraction(numerator * k.denominator - denominator * k.numerator, denominator * k.denominator).norm();
    }
    constexpr fraction operator* (const fraction& k) const noexcept {
        return fraction(numerator * k.numerator, denominator * k.denominator).norm();
    }
    constexpr fraction operator/ (const fraction& k) const noexcept {
        return fraction(numerator * k.denominator, denominator * k.numerator).norm();
    }
    constexpr bool operator == (const fraction& k) const noexcept {
        return ((numerator == k.numerator) && (denominator == k.denominator));
    }
    constexpr bool operator != (const fraction& k) const noexcept {
        return ((numerator != k.numerator) || (denominator != k.denominator));
    }
    friend constexpr istream& operator >> (istream& is, fraction& k) noexcept {
        is >> k.numerator >> k.denominator;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const fraction& k) noexcept {
        return os << k.numerator << "/" << k.denominator;
    }
};
using frac = fraction<long long>;
暂无评论

发送评论 编辑评论 


				

			

						

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
						

					

				

				

					

						

							

								
							

							
							
													

					

				

			

			
			

				

					
				

			

				

											

							
							
						

																					

							
							
						

										
					
											
						

	

		
|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
（╯‵□′）╯︵┴─┴
￣﹃￣
(/ω＼)
∠( ᐛ 」∠)＿
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ｀)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ(￣∇￣o)
ヾ(´･ ･｀｡)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò｡)
Σ(っ °Д °;)っ
( ,,´･ω･)ﾉ"(´っω･｀｡)
╮(╯▽╰)╭
o(*////▽////*)q
＞﹏＜
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
	

	

		
颜文字
Emoji
小恐龙
花!
	


									

			

			
			
		

	






上一篇
下一篇