zoukankan html css js c++ java

数论板子（待续）

一、质数：

1、试除法判质数：

bool is_prime(int n){
    if(n < 2) return false;
    for(int i = 2; i <= n / i; i ++){
        if(n % i == 0)
            return false;
    }
    return true;
}

2、试除法分解质因数

void divide(int n){
    for(int i = 2; i <= n / i; i ++){
        if(n % i == 0){
            int s = 0;
            while(n % i == 0){
                n /= i;
                s ++;
            }
            v.push_back({i, s});
        }
    }
    if(n > 1){
        v.push_back({n, 1});
    }
}

3、朴素筛法

把所有数的倍数全部筛掉

int primes[N], cnt;
bool st[N];

void get_primes(int n){
    for (int i = 2; i <= n; i++){
        if (!st[i])
            primes[cnt++] = i;
        for (int j = i + i; j <= n; j += i){
            st[j] = true;
        }
    }
}

4、埃氏筛法

把质数的倍数全部筛掉

void get_primes(int n)
{
    for (int i = 2; i <= n; i++)
    {
        if (!st[i])
        {
            primes[cnt++] = i;
            for (int j = i + i; j <= n; j += i)
            {
                st[j] = true;
            }
        }
    }
}

5、线性筛法：

每个合数只会被最小质因子筛掉

void get_primes(int n)
{
    for (int i = 2; i <= n; i++)
    {
        if (!st[i])
        {
            primes[cnt++] = i;
        }
        for (int j = 0; primes[j] <= n / i; j++)
        {
            st[i * primes[j]] = true;
            if (i % primes[j] == 0)
            {
                break;
            }
        }
    }
}

二、约数

1、试除法求所有约数

void get_d(int n)
{
    for(int i = 1; i <= n / i; i ++){
        if(n % i == 0){
            v.push_back(i);
            if(i != n / i){
                v.push_back(n / i);
            }
        }
    }
}

2、约数个数

设n = p₁^k₁p₂^k₂ p₃^k₃...... p_i^k_i

则ans = (k₁ + 1) * (k₂ + 1) * ... * (k_i + 1)

三、欧拉函数

１、求欧拉函数

int phi(int n){
    int res = n;
    for(int i = 2; i <= n / i; i ++){
        if(n % i == 0){
            res = res / i * (i - 1);
            while(n % i == 0){
                n /= i;
            }
        }
    }
    if(res > 1){
        res = res / n * (n - 1);
    }
    return res;
}

２、欧拉筛

void get_eulers(int n)
{
    for (int i = 2; i <= n; i++)
    {
        if (!st[i])
        {
            primes[cnt++] = i;
            phi[i] = i - 1;
        }
        for (int j = 0; primes[j] <= n / i; j++)
        {
            st[i * primes[j]] = true;
            if (i % primes[j] == 0)
            {
                phi[i * primes[j]] = primes[j] * phi[i];
                break;
            }
            phi[i * primes[j]] = (primes[j] - 1) * phi[i];
        }
    }
}

四、快速幂

ll qmi(ll a, ll b, ll p){
    ll res = 1;
    while(b){
        if(b & 1){
            res = res * a % p;
        }
        a = a * a % p;
        b >>= 1;
    }
    return res;
}

五、欧几里得

１．欧几里得

ll gcd(ll a, ll b){
    return b ? gcd(b, a % b) : a;
}

２．扩展欧几里得

ax + by = c

d = gcd(a, b)

x = x0 + k * b / d

y = y0 - k * a / d

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

六、中国剩余定理

１．mi 两两互质

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

ll china(int n, ll M){
    ll ans = 0;
    for(int i = 0; i < n; i ++){
        ll Mi = M / m[i];
        ll ti, x;
        exgcd(Mi, m[i], ti, x);
        ans += a[i] * Mi * ti;
    }
    return (ans % M + M) % M;
}

2.　不互质

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

ll n, m[N], a[N], m1, m2, a1, a2, x, y, g, c;

ll china()
{
    m1 = m[0], a1 = a[0];
    for (int i = 1; i < n; i++)
    {
        m2 = m[i], a2 = a[i];
        c = a2 - a1;
        g = exgcd(m1, m2, x, y);
        if (c % g)
        {
            a1 = -1;
            break;
        }
        x *= c / g;
        ll t = m2 / g;
        x = (x % t + t) % t;
        a1 += m1 * x;
        m1 *= m2 / g;
    }

    return a1;
}

七、高斯消元

1. 求浮点数解（保证有唯一解）

void gauss()
{
    int c, r;
    for (c = 1, r = 1; c <= n; c++)
    {
        int t = r;
        for (int i = r; i <= n; i++)
        {
            if (fabs(b[i][c]) > fabs(b[t][c]))
            {
                t = i;
            }
        }

        for (int i = c; i <= n + 1; i++)
        {
            swap(b[t][i], b[r][i]);
        }
        for (int i = n + 1; i >= c; i--)
        {
            b[r][i] /= b[r][c];
        }

        for (int i = r + 1; i <= n; i++)
        {
            if (fabs(b[i][c]) > eps)
            {
                for (int j = n + 1; j >= c; j--)
                {
                    b[i][j] -= b[r][j] * b[i][c];
                }
            }
        }
        r++;
    }

    for (int i = n; i >= 1; i--)
    {
        for (int j = i + 1; j <= n; j++)
        {
            b[i][n + 1] -= b[i][j] * b[j][n + 1];
        }
    }
}

2.异或（开关问题）

int gauss(){
    int r, c;
    for(r = 1, c = 1; c <= n; c ++){
        int t = r;
        for(int i = r + 1; i <= n; i ++){
            if(a[i][c]){
                t = i;
                break;
            }
        }
        if(!a[t][c]) continue;
        for(int i = c; i <= n + 1; i ++){
            swap(a[t][i], a[r][i]);
        }
        for(int i = r + 1; i <= n; i ++){
            if(a[i][c]){
                for(int j = n + 1; j >= c; j --){
                    a[i][j] ^= a[r][j];
                }
            }
        }
        r ++;
    }

    int res = 1;
    if(r < n + 1){
        for(int i = r; i <= n; i ++){
            if(a[i][n + 1]) return -1;
            res *= 2;
        }
    }
    return res;
}

八、整数分块

for(int l = 1, r; l <= m; l = r + 1){
     r = m / (m / l);
    res = res + (ll)(sum[r] - sum[l - 1]) * (m / l) * (m / l); 
}

九、快速傅里叶变换 FFT

const int N = 3e6 + 10;
const double PI = acos(-1.0);

int n, m;
ll ans[N];
int rev[N];

struct Complex{
    double x, y;

    Complex(double _x = 0.0, double _y = 0.0){
        x = _x, y = _y;
    }

    Complex operator + (const Complex & b){
        return Complex(x + b.x, y + b.y);
    }

    Complex operator - (const Complex &b){
        return Complex(x - b.x, y - b.y);
    }

    Complex operator * (const Complex &b){
        return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
    }
};

void change(Complex y[], int len){
    
    for(int i = 0; i < len; i ++){
        rev[i] = rev[i >> 1] >> 1;
        if(i & 1) rev[i] |= (len >> 1);
    }

    for(int i = 0; i < len; i ++){
        if(i < rev[i]) swap(y[i], y[rev[i]]);
    }
}

void fft(Complex y[], int len, int on){
    change(y, len);

    for(int h = 2; h <= len; h <<= 1){
        Complex wn(cos(2 * PI / h), sin(2 * PI * on / h));

        for(int j = 0; j < len; j += h){
            Complex w(1, 0);

            for(int k = j; k < j + h / 2; k ++){
                Complex u = y[k];
                Complex t = w * y[k + h / 2];

                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;

            }
        }
    }
    if(on == -1){
        for(int i = 0; i < len; i ++){
            y[i].x /= len;
        }
    }
}

　　FFT求快速幂

void qmi(int b){
    while(b){
        if(b & 1){
            fft(x1, len, 1);
            fft(a, len, 1);
 
            for(int i = 0; i < len; i ++) x1[i] = x1[i] * a[i];
            fft(x1, len, -1);
            fft(a, len, -1);
 
            for(int i = 0; i < len; i ++) mid[i] = (int)(x1[i].x + 0.5);
            for(int i = 0; i < len; i ++){
                mid[i + 1] += mid[i] / 10;
                mid[i] %= 10;
            }
            for(int i = 0; i < len; i ++){
                x1[i] = Complex(mid[i], 0);
            }
 
            for(int i = 0; i < len; i ++) mid[i] = (int)(a[i].x + 0.5);
            for(int i = 0; i < len; i ++){
                mid[i + 1] += mid[i] / 10;
                mid[i] %= 10;
            }
            for(int i = 0; i < len; i ++){
                a[i] = Complex(mid[i], 0);
            }
        }
        fft(a, len, 1);
 
        for(int i = 0; i < len; i ++){
            a[i] = a[i] * a[i];
        }
 
        fft(a, len, -1);
 
        for(int i = 0; i < len; i ++) mid[i] = (int)(a[i].x + 0.5);
        for(int i = 0; i < len; i ++){
            mid[i + 1] += mid[i] / 10;
            mid[i] %= 10;
        }
        
        for(int i = 0; i < len; i ++){
            a[i] = Complex(mid[i], 0);
        }
         
        b >>= 1;
    }
}

十、快速数论变换 NTT

void ntt(ll y[], int len, int on){
    change(y, len);

    for(int h = 2; h <= len; h <<= 1){
        ll gn = qmi(3, (mod - 1) / h, mod);   //3是原根
        if(on == -1) gn = qmi(gn, mod - 2, mod);
        for(int j = 0; j < len; j += h){
            ll g = 1;
            for(int k = j; k < j + h / 2; k ++){
                ll u = y[k];
                ll t = g * y[k + h / 2] % mod;
                y[k] = (u + t) % mod;
                y[k + h / 2] = (u - t + mod) % mod;
                g = g * gn % mod;
            }
        }
    }

    if(on == -1){
        ll inv = qmi(len, mod - 2, mod);
        for(int i = 0; i < len; i ++){
            y[i] = y[i] * inv % mod;
        }
    }
}

NTT求快速幂

void ntt_qmi(ll b, int len){
    int mm = m - 1;
    while(b){
        ntt(a, len, 1);
        if(b & 1){
            ntt(res, len, 1);
            for(int i = 0; i < len; i ++){
                res[i] = res[i] * a[i] % mod;
            }
            ntt(res, len, -1);

            for(int i = mm; i < len; i ++){
                res[i % mm] = (res[i % mm] + res[i]) % mod;
                res[i] = 0;
            }
        }
        for(int i = 0; i < len; i ++){
            a[i] = a[i] * a[i] % mod;
        }
        ntt(a, len, -1);
        for(int i = mm; i < len; i ++){
            a[i % mm] = (a[i % mm] + a[i]) % mod;
            a[i] = 0;
        }
        b >>= 1;
    }
}

三模数NTT+快速幂

#include <map>
#include <set>
#include <array>
#include <queue>
#include <stack>
#include <cmath>
#include <vector>
#include <cstdio>
#include <cstring>
#include <sstream>
#include <iostream>
#include <stdlib.h>
#include <algorithm>
#include <unordered_map>

using namespace std;

typedef long long ll;
typedef pair<int, int> PII;

#define sd(a) scanf("%d", &a)
#define sdd(a, b) scanf("%d%d", &a, &b)
#define slld(a) scanf("%lld", &a)
#define slldd(a, b) scanf("%lld%lld", &a, &b)
#define m1 998244353
#define m2 469762049
#define m3 1004535809

const int N = 3e2 + 10;
const int M = 2e7 + 20;
const int mod = 20170408;
const int INF = 0x3f3f3f3f;
const double PI = acos(-1.0);
const int Mod[] = {998244353, 469762049, 1004535809};

int n, m, p;
int rev[N];
ll vis[N], h[N];
int primes[M], cnt = 0;
bool st[M];

void get(int n)
{
    st[1] = true;
    for (int i = 2; i <= n; i++)
    {
        if (!st[i])
            primes[cnt++] = i;
        for (int j = 0; primes[j] <= n / i; j++)
        {
            st[i * primes[j]] = true;
            if (i % primes[j] == 0)
            {
                break;
            }
        }
    }
}

ll qmi(ll a, ll b, ll p)
{
    ll res = 1;
    while (b)
    {
        if (b & 1)
            res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}

void change(ll y[], int len)
{
    for (int i = 0; i < len; i++)
    {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1)
            rev[i] |= (len >> 1);
    }

    for (int i = 0; i < len; i++)
    {
        if (i < rev[i])
            swap(y[i], y[rev[i]]);
    }
}

void ntt(ll y[], int len, int on, ll MOD)
{
    change(y, len);
    for (int h = 2; h <= len; h <<= 1)
    {
        ll wn = qmi(3, (MOD - 1) / h, MOD);
        if (on == -1)
            wn = qmi(wn, MOD - 2, MOD);

        for (int j = 0; j < len; j += h)
        {
            ll w = 1;
            for (int k = j; k < j + h / 2; k++)
            {
                ll u = y[k];
                ll t = w * y[k + h / 2] % MOD;
                y[k] = (u + t) % MOD;
                y[k + h / 2] = (u - t + MOD) % MOD;
                w = w * wn % MOD;
            }
        }
    }

    if (on == -1)
    {
        ll inv = qmi(len, MOD - 2, MOD);
        for (int i = 0; i < len; i++)
        {
            y[i] = y[i] * inv % MOD;
        }
    }
}

ll mult(ll a, ll b, ll p)
{
    ll res = 0;
    while (b)
    {
        if (b & 1)
            res = (res + a) % p;
        a = (a + a) % p;
        b >>= 1;
    }
    return res;
}

ll A[N], B[N], C[N], D[N];
void mul(ll a[], ll b[], ll res[], ll len)
{
    memcpy(A, a, sizeof(A));
    memcpy(B, a, sizeof(B));
    memcpy(C, a, sizeof(C));
    memcpy(D, b, sizeof(D));

    ntt(A, len, 1, Mod[0]);
    ntt(D, len, 1, Mod[0]);
    for (int i = 0; i < len; i++)
    {
        A[i] = A[i] * D[i] % Mod[0];
    }
    ntt(A, len, -1, Mod[0]);

    memcpy(D, b, sizeof(D));
    ntt(B, len, 1, Mod[1]);
    ntt(D, len, 1, Mod[1]);
    for (int i = 0; i < len; i++)
    {
        B[i] = B[i] * D[i] % Mod[1];
    }
    ntt(B, len, -1, Mod[1]);

    memcpy(D, b, sizeof(D));
    ntt(C, len, 1, Mod[2]);
    ntt(D, len, 1, Mod[2]);
    for (int i = 0; i < len; i++)
    {
        C[i] = C[i] * D[i] % Mod[2];
    }
    ntt(C, len, -1, Mod[2]);

    ll M12 = 1ll * m1 * m2;
    ll inv2 = qmi(m2, m1 - 2, m1);
    ll inv1 = qmi(m1, m2 - 2, m2);
    ll mul2 = 1ll * m2 * inv2 % M12;
    ll mul1 = 1ll * m1 * inv1 % M12;
    ll inv = qmi(M12 % m3, m3 - 2, m3);
    ll m12 = M12 % mod;
    ll c1, c2, c3, c4, q;

    for (int i = 0; i <= (p << 1); i++)
    {
        c1 = A[i], c2 = B[i], c3 = C[i];
        c4 = (mult(c1, mul2, M12) + mult(c2, mul1, M12)) % M12;
        q = ((c3 - c4) % m3 + m3) % m3 * inv % m3;
        res[i] = (q * m12 % mod + c4) % mod;
    }
    for (int i = p; i < len; i++)
    {
        res[i % p] = (res[i % p] + res[i]) % mod;
        res[i] = 0;
    }
}

ll res[N];

void qmi_ntt(ll y[], int len, int n)
{
    memset(res, 0, sizeof(res));
    res[0] = 1;
    while (n)
    {
        if (n & 1)
        {
            mul(res, y, res, len);
        }
        mul(y, y, y, len);
        n >>= 1;
    }
}

ll mid[N], ans[3], ans1, ans2;

void solve()
{
    cin >> n >> m >> p;

    get(m);
    for (int i = 1; i <= m; i++)
    {
        vis[i % p]++;
        if (st[i])
            h[i % p]++;
    }

    int len = 1;
    while (len <= p + p - 1)
        len <<= 1;

    qmi_ntt(vis, len, n);

    ans1 = res[0];

    qmi_ntt(h, len, n);

    ans1 = (ans1 - res[0] + mod) % mod;
    cout << ans1 << "
";
}

int main()
{
#ifdef ONLINE_JUDGE
#else
    freopen("/home/jungu/code/in.txt", "r", stdin);
    // freopen("/home/jungu/桌面/11.21/2/in9.txt", "r", stdin);
#endif
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int T = 1;
    // sd(T);
    // cin >> T;
    while (T--)
    {
        solve();
    }

    return 0;
}

查看全文

相关阅读:
Spring框架的反序列化远程代码执行漏洞分析（转）
JVM调优总结（转）
创造与熟练与微创造
 一次 react-router 中遇到的小坑
 mac中强大的快捷键
 数组思维 -- join的一些用法感悟
 undefined null 与字符串相加
 mobx中的数组需要注意的地方
 更换鼠标垫(鼠标)的心路历程
 Object.defineProperty 中的 writable 和 configurable 和 enumerable 的理解

原文地址：https://www.cnblogs.com/jungu/p/13387553.html