ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

Number 数论

2021-10-19 11:02:11  阅读:158  来源: 互联网

标签:return 数论 ll Number int const Mint mod


Number 数论

素数筛

欧拉筛

对if(i % prime[j] == 0) break;的解释

当i % prime[j] == 0时有 i = k * prime[j]; 若j++有 i * prime[j + 1] = k * prime[j] * prime[j + 1] 也是prime[j]的因子,导致重复筛

const int maxn = 1e7 + 8;
vector<int > prime(maxn, 0);
bitset<maxn> vis(0);
void getPrime() {
	for(int i = 2;i < maxn;i++) {
		if(vis[i] == 0) prime[++prime[0]] = i;
		for(int j = 1;j <= prime[0];j++) {
			if(i * prime[j] >= maxn) break;
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) break;
		}
	}
}

埃式筛

void Prime(){
    vector<int > vis(maxn, 0); //初始化都是素数
    vis[0] = vis[1] = 1;    //0 和 1不是素数
    for (int i = 2; i <= maxn; i++) {
        if (!vis[i]) {      //如果i是素数,让i的所有倍数都不是素数
            for (int j = i * i; j <= maxn; j += i) { 
                vis[j] = 1;
            }
        }
    }
}

欧拉函数

​ 用处:1:求1-n内与n互质的数量

​ 2:求合数的逆元

​ 3:欧拉降幂

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

线性筛欧拉函数

const int maxn = 100005;
vector<int > prime(maxn, 0), vis(maxn, 0), phi(maxn, 0);
void getPrime() {
    phi[1] = 1;
	for(int i = 2;i < maxn;i++) {
		if(vis[i] == 0) prime[++prime[0]] = i, phi[i] = i - 1;
		for(int j = 1;j <= prime[0];j++) {
			if(i * prime[j] >= maxn) break;
			vis[i * prime[j]] = 1;
			if(i % prime[j] != 0) phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            //对于任意的a,b若gcd(a, b) == 1 -> phi(a * b) = phi(a) * phi(b)
			else {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
		}
	}
}

任意区间素数打表

const int M=1e6+5,N=(1<<16);
int prime[N+5],is[N+5],tot;
void getPrime(){
    is[1]=1;//mmpaaaaaaaaaaaaaaaa
    for(int i=2;i<=N;++i){
        if(!is[i]) prime[++tot]=i;
        for(int j=1;j<=tot&&(ll)i*prime[j]<=N;++j){
            is[i*prime[j]]=1;
            if(i%prime[j]==0) break;
        }
    }
}
int issp[M];
int main(){
    getPrime();
    int T,cas=0;scanf("%d",&T);
    while(T--){
        m(issp,0);
        int a,b;scanf("%d%d",&a,&b);
        if(a<=N&&b<=N) {
            int ans=0;
            for(int i=a;i<=b;++i) if(!is[i]) ++ans;
            printf("Case %d: %d\n",++cas,ans);
            continue;
        }
        int ans=0;
        if(a<=N) {
            for(int i=a;i<=N;++i) if(!is[i]) ++ans;
            a=N+1;
        }
        for(int i=1;i<=tot;++i) {
            ll l=a/prime[i]*prime[i];//左端点l.
            if(l<a) l+=prime[i];
            if(l==prime[i]) l+=prime[i];
            for(ll j=l;j<=b;j+=prime[i]) issp[j-a]=1;
        }
        for(int j=a;j<=b;++j) if(!issp[j-a]) ++ans;
        printf("Case %d: %d\n",++cas,ans);
    }
}

大数素数测试

Miller_Rabin

#include <bits/stdc++.h>

using namespace std;
using ll = __int128;

const int S = 8; //随机算法判定次数一般 8~10 就够了
// 计算 ret = (a*b)%c a,b,c < 2^63
__int128 mult_mod(__int128 a, __int128 b, __int128 c) {
    a %= c;
    b %= c;
    __int128 ret = 0;
    __int128 tmp = a;

    while (b) {
        if (b & 1) {
            ret += tmp;

            if (ret > c)
                ret -= c;//直接取模慢很多
        }

        tmp <<= 1;

        if (tmp > c)
            tmp -= c;

        b >>= 1;
    }

    return ret;
}
// 计算 ret = (a^n)%mod
__int128 pow_mod(__int128 a, __int128 n, __int128 mod) {
    __int128 ret = 1;
    __int128 temp = a % mod;

    while (n) {
        if (n & 1)
            ret = mult_mod(ret, temp, mod);

        temp = mult_mod(temp, temp, mod);
        n >>= 1;
    }

    return ret;
}
bool check(__int128 a, __int128 n, __int128 x, __int128 t) {
    __int128 ret = pow_mod(a, x, n);
    __int128 last = ret;

    for (int i = 1; i <= t; i++) {
        ret = mult_mod(ret, ret, n);

        if (ret == 1 && last != 1 && last != n - 1)
            return true;//合数

        last = ret;
    }

    if (ret != 1)
        return true;
    else
        return false;
}
//**************************************************
// Miller_Rabin 算法
// 是素数返回 true,(可能是伪素数)
// 不是素数返回 false
//**************************************************
bool MR(__int128 n) {
    if (n < 2)
        return false;

    if (n == 2)
        return true;

    if ((n & 1) == 0)
        return false;//偶数

    __int128 x = n - 1;
    __int128 t = 0;

    while ((x & 1) == 0) {
        x >>= 1;
        t++;
    }

    srand(time(NULL));

    for (int i = 0; i < S; i++) {
        __int128 a = rand() % (n - 1) + 1;

        if (check(a, n, x, t))
            return false;
    }

    return true;
}
inline __int128 read() {
    __int128 x = 0, f = 1;
    char ch = getchar();

    while (ch < '0' || ch > '9') {
        if (ch == '-')
            f = -1;

        ch = getchar();
    }

    while (ch >= '0' && ch <= '9') {
        x = x * 10 + ch - '0';
        ch = getchar();
    }

    return x * f;
}
int main() {
    __int128 n = read();

    if (MR(n))
        puts("Yes");
    else
        puts("No");

    return 0;
}

Miller_Rabin(Easy ver)

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

ll mul(ll a, ll b, ll p) {
	return (__int128)a * b % p;
}

ll qp(ll x, ll y, ll p) {
	ll z = 1;
	for (; y; y >>= 1, x = mul(x, x, p)) {
		if (y & 1) z = mul(z, x, p);
	}
	return z;
}

bool MR(ll x, ll b) {
	for (ll k = x - 1; k; k >>= 1) {
		ll cur = qp(b, k, x);
		if (cur != 1 && cur != x - 1) return 0;
		if (k & 1 || cur == x - 1) return 1;
	}
	return true;
}

bool test(ll x) {
	if (x == 1) return 0;
	static ll p[] = {2, 3, 5, 7, 17, 19, 61};
	for (ll y : p) {
		if (x == y) return 1;
		if (!MR(x, y)) return 0;
	}
	return 1;
}

int main() {
	ll n;
	while (scanf("%lld", &n) != EOF) {
		puts(test(n) ? "Y" : "N");
	}
	return 0;
}

大数分解

逆元

Exgcd求a在b意义下的逆元

i f ( g c d ( a , b ) = = 1 ) if(gcd(a, b)==1) if(gcd(a,b)==1)

int exgcd(int a, int b, int &x, int &y) {
	if(b == 0) {
		x = 1;
		y = 0;
		return a;
	}
	else {
		int gcd = exgcd(b, a % b, x, y);
		int t = x;
		x = y;
		y = t - a / b * y;
		return gcd;
	}
}
int getInv(int a, int b){
	int x, y;
	exgcd(a, b, x, y);
	return (x + b) % b;
}

费马小定理求a在mod p意义下的逆元(p是素数)

a ∗ a p − 2 = 1 ( m o d p ) a * a^{p - 2} = 1{(modp)} a∗ap−2=1(modp)

int qpow(int a, int b, int mod) {
	int ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % mod;
		b >>= 1, a = a * a % mod;
	}
	return ans;
}
int getInv(int x, int mod) {
	return qpow(x, mod - 2, mod);
}

费马小定理求a在mod p意义下的逆元(p不是素数)

a ∗ a p h i ( p ) − 1 = 1 ( m o d p ) a*a^{phi(p) - 1} = 1(modp) a∗aphi(p)−1=1(modp)

int qpow(int a, int b, int mod) {
	int ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % mod;
		b >>= 1, a = a * a % mod;
	}
	return ans;
}
int phi(int x) {
	int ans = x;
	for(int i = 2;i * i <= x;i++) {
		if(x % i == 0) {
			ans = ans * (i - 1) / i;
			while(x % i == 0) x /= i;
		}
	}
	if(x > 1) ans = ans * (x - 1) / x;
	return ans;
}
int getInv(int x, int mod) {
	return qpow(x, phi(mod) - 1, mod);
}

O(n)求逆元

vector<ll> inv(3000006, 0);
void liner_inv(int n, int p) {
	inv[1] = 1;
	for(int i = 2;i <= n;i++) inv[i] = 1ll * (p - p / i) * inv[p % i] % p;
}

O(n + log(mod))求任意n个数的逆元

const ll mod  = 1e9 + 7;
const ll p    = 998244353;
const ll maxn = 1e5 + 7;

vector<ll> a(maxn, 0), pre(maxn, 1), pre_inv(maxn, 0), inv(maxn, 0);
//		   数组元素    前缀积		 前缀积的逆元	   每个数的逆元

ll qpow(ll a, ll b) {
	ll ret = 1;
	while(b) {
		if(b & 1) ret = ret * a % mod;
		b >>= 1, a = a * a % mod;
	}
	return ret;
}

void solve(int n) {//数组长度参数
	for(int i = 1;i <= n;i++) cin >> a[i];
	for(int i = 1;i <= n;i++) pre[i] = pre[i - 1] * a[i] % mod;
	pre_inv[n] = qpow(pre[n], mod - 2);
	for(int i = n;i >= 1;i--) pre_inv[i - 1] = pre_inv[i] * a[i] % mod;
	for(int i = 1;i <= n;i++) inv[i] = pre_inv[i] * pre[i - 1] % mod;
}

小数循环节

int c(int tmp) {
    while(tmp % 2 == 0) tmp /= 2;
	while(tmp % 5 == 0) tmp /= 5;
	int len = 0, e = 1;
	if(tmp == 1) return -1;
	while(1) {
		e = e * 10 % tmp;
		if(e == 1) break;
		len++;
	}
	return len;
}

组合数

C n m = n ! m ! ∗ ( n − m ) ! C_n^m = \frac{n!}{m!*(n - m)!} Cnm​=m!∗(n−m)!n!​

正常版本

const int maxn = 2e5 + 7;
vector<ll> fac(maxn, 1);

void init(ll p) {//if module is different
	for(ll i = 1;i <= 200000;i++) fac[i] = fac[i - 1] * i % p;
}

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

ll getInv(ll x, ll p) {
	if(x == 1) return 1;
	return qpow(x, p - 2, p);
}

ll c(ll n, ll m, ll p) {
	if(!m) return 1;
	if(m > n) return 0;
	return fac[n] * getInv(fac[m], p) % p * getInv(fac[n - m], p) % p;
}

Mint版本

constexpr int mod = 1000000007;
constexpr int maxn = 1e5 + 6;
// assume -mod <= x < 2mod
int norm(int x) {if(x < 0) x += mod; if (x >= mod) x -= mod; return x; }
template<class T>T power(T a, ll b){T res = 1;for (; b; b /= 2, a *= a)if (b & 1)res *= a;return res;}
struct Mint {
    int x;Mint(int x = 0) : x(norm(x)){}
    int val() const {return x;}
    Mint operator-() const {return Mint(norm(mod - x));}
    Mint inv() const {assert(x != 0);return power(*this, mod - 2);}
    Mint &operator*=(const Mint &rhs) { x = ll(x) * rhs.x % mod; return *this;}
    Mint &operator+=(const Mint &rhs) { x = norm(x + rhs.x); return *this;}
    Mint &operator-=(const Mint &rhs) { x = norm(x - rhs.x); return *this;}
    Mint &operator/=(const Mint &rhs) {return *this *= rhs.inv();}
    friend Mint operator*(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res *= rhs; return res;}
    friend Mint operator+(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res += rhs; return res;}
    friend Mint operator-(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res -= rhs; return res;}
    friend Mint operator/(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res /= rhs; return res;}
};

vector<Mint> fac(maxn, 1);

void init() {
	for(ll i = 1;i < maxn;i++) fac[i] = fac[i - 1] * i;
}

Mint c(ll n, ll m) {
	return fac[n] / fac[m] / fac[n - m];
}

容斥原理

CRT

using ll = long long;

ll qpow(ll a, ll b, ll mod) {
	ll ret = 1;
	while(b) {
		if(b & 1) ret = ret * a % mod;
		b >>= 1;
		a = a * a % mod;
	}
	return ret;
}

ll getInv(ll x, ll mod) {
	return qpow(x, mod - 2, mod);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

	ll n, mul = 1, sum = 0;
	cin >> n;
	vector<ll> a(n + 1, 0), mu(n + 1, 0);
	for(int i = 1;i <= n;i++) cin >> a[i] >> mu[i];
	for(int i = 1;i <= n;i++) mul *= a[i];
	for(int i = 1;i <= n;i++) {
		ll lcm = mul / a[i];
		ll val = getInv(lcm, a[i]) * mu[i] * lcm;
		sum += val;
	}
	cout << sum % mul << '\n';

    return 0;
}

Lucas/ExLucas

Lucas

const int maxn = 2e5 + 7;
vector<ll> fac(maxn, 1);

void init(ll p) {//if module is different
	for(ll i = 1;i <= 200000;i++) fac[i] = fac[i - 1] * i % p;
}

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

ll getInv(ll x, ll p) {
	if(x == 1) return 1;
	return qpow(x, p - 2, p);
}

ll c(ll n, ll m, ll p) {
	if(!m) return 1;
	if(m > n) return 0;
	return fac[n] * getInv(fac[m], p) % p * getInv(fac[n - m], p) % p;
}

ll lucas(ll n, ll m, ll p) {
	if(n == 0) return 1;
	return lucas(n / p, m / p, p) * c(n % p, m % p, p) % p;
}

Mint version

constexpr int mod = 1000000007;
const int maxn = 1e5 + 6;
// assume -mod <= x < 2mod
int norm(int x) {if(x < 0) x += mod; if (x >= mod) x -= mod; return x; }
template<class T>T power(T a, ll b){T res = 1;for (; b; b /= 2, a *= a)if (b & 1)res *= a;return res;}
struct Mint {
    int x;Mint(int x = 0) : x(norm(x)){}
    int val() const {return x;}
    Mint operator-() const {return Mint(norm(mod - x));}
    Mint inv() const {assert(x != 0);return power(*this, mod - 2);}
    Mint &operator*=(const Mint &rhs) { x = ll(x) * rhs.x % mod; return *this;}
    Mint &operator+=(const Mint &rhs) { x = norm(x + rhs.x); return *this;}
    Mint &operator-=(const Mint &rhs) { x = norm(x - rhs.x); return *this;}
    Mint &operator/=(const Mint &rhs) {return *this *= rhs.inv();}
    friend Mint operator*(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res *= rhs; return res;}
    friend Mint operator+(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res += rhs; return res;}
    friend Mint operator-(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res -= rhs; return res;}
    friend Mint operator/(const Mint &lhs, const Mint &rhs) {Mint res = lhs; res /= rhs; return res;}
};

vector<Mint> fac(maxn, 1);

void init() {
	for(ll i = 1;i < maxn;i++) fac[i] = fac[i - 1] * i;
}

Mint c(ll n, ll m) {
	return fac[n] / fac[m] / fac[n - m];
}

Mint lucas(ll n, ll m) {
	if(m == 0) return 1;
	if(n == 0) return 1;
	if(n < mod && m < mod) return c(n, m);
	return lucas(n / mod, m / mod) * c(n % mod, m % mod);
}

Gcd

Exgcd

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

矩阵快速幂

实例:求斐波那契数列的第n项

$$
\begin{bmatrix}
{}&{}&{}&{}\
{}&{}&{}&{}\
{}&{}&{}&{}\
{}&{}&{}&{}\
\end{bmatrix}

\begin{bmatrix}
{}&{}&{}\
{}&{}&{}\
{}&{}&{}\
\end{bmatrix}

\begin{bmatrix}
{}&{}\
{}&{}\
\end{bmatrix}
$$

$$
\begin{bmatrix}
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}&{f_{}}\
\end{bmatrix}

\begin{bmatrix}
{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}\
{f_{}}&{f_{}}&{f_{}}\
\end{bmatrix}

\begin{bmatrix}
{f_{}}&{f_{}}\
{f_{}}&{f_{}}\
\end{bmatrix}
$$

$$
\begin{bmatrix}
{f_{}}\
{f_{}}\
{f_{}}\
{f_{}}\
\end{bmatrix}

\begin{bmatrix}
{f_{}}\
{f_{}}\
{f_{}}\
\end{bmatrix}

\begin{bmatrix}
{f_{}}\
{f_{}}\
\end{bmatrix}
$$

$$
f_{i} = a*f_{i - 1} + i:
\begin{bmatrix}
{a}&{1}&{0}\
{0}&{1}&{1}\
{0}&{0}&{1}\
\end{bmatrix}
\cdot
\begin{bmatrix}
{f_{i-1}}\
{{i}}\
{{1}}\
\end{bmatrix}

\begin{bmatrix}
{f_{i}}\
{{i+1}}\
{{1}}\
\end{bmatrix}
$$

$$
f_{i} = f_{i - 1} + f_{i - 2}:
\begin{bmatrix}
{1}&{1}\
{1}&{0}\
\end{bmatrix}
\cdot
\begin{bmatrix}
{f_{n-1}}\
{f_{n-2}}\
\end{bmatrix}

\begin{bmatrix}
{f_{n}}\
{f_{n-1}}\
\end{bmatrix}
$$

#include<bits/stdc++.h>

using namespace std;
using ll = long long;

int n_;       //Size of Matrix
const ll  p = 1e9 + 7;
const int N = 10;

struct matrix { 
	ll a[N][N]{};

	matrix init() {		//to E
		matrix E;
		for(register int i = 0;i < n_;i++) E.a[i][i] = 1;
		return E;
	}
	
	matrix f1_() {		//f(n) = f(n - 1) + f(n - 2),  f[n] = qpow(A, n - 2);
		matrix ret;
		ret.a[0][0] = 1, ret.a[0][1] = 1;
		ret.a[1][0] = 1, ret.a[1][1] = 0;
		return ret;
	}

	matrix f2_() {		//f(n) = f(n - 1) + f(n - 2) + f(n - 3);
		matrix ret;
		ret.a[0][0] = 1, ret.a[0][1] = 1, ret.a[0][2] = 1;
		ret.a[1][0] = 1, ret.a[1][1] = 0, ret.a[1][2] = 0;
		ret.a[2][0] = 0, ret.a[2][1] = 1, ret.a[2][2] = 0;
		return ret;
	}

	matrix f3_() {
		matrix ret;
		//to construct ...
		return ret;
	}

	matrix operator * (const matrix& b) {
		matrix ret;
		for(register int i = 0;i < n_;i++) 
			for(register int j = 0;j < n_;j++) 
				for(register int k = 0;k < n_;k++) 
					ret.a[i][j] = (ret.a[i][j] + a[i][k] * b.a[k][j]) % p;
		return ret;
	}
	matrix qpow(matrix a, ll b) {
		matrix ans = init();
		while(b) {
			if(b & 1) ans = ans * a;
			b >>= 1,  a = a * a;
		}
		return ans;
	}
	void show() {
		for(register int i = 0;i < n_;i++) 
			for(register int j = 0;j < n_;j++) 
				cout << a[i][j] << " \n"[j == n_ - 1];
		return ;
	}
};

ll qpow_(ll a, ll b, ll mod) {
	ll ret = 1;
	while(b){
		if(b & 1) ret = ret * a % mod;
		b >>= 1, a = a * a % mod;
	}
	return ret;
}

void Fib_n() {
	int n; cin >> n;
	if(n == 1 || n == 2) {
		cout << 1 << '\n';
		return ;
	}
	matrix f = f.f1_();
	n_ = 2;
	f = f.qpow(f, n - 2);
	cout << f.a[0][0] + f.a[0][1] << '\n';
}

int main() {
	Fib_n();

    return 0;
}

约瑟夫环

FFT

const int MAXN = 1e6 + 10;
const double Pi = acos(-1.0);
struct complex{
	double x, y;
	complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
}a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);} 
int l, r[MAXN], ans[MAXN];
int limit = 1;
void fft(complex *A, int type) {
	for(int i = 0; i < limit; i++) if (i < r[i]) swap(A[i], A[r[i]]); 
	for (int mid = 1; mid < limit; mid <<= 1) { 
		complex Wn( cos(Pi / mid) , type * sin(Pi / mid) ); 
		for (int R = mid << 1, j = 0; j < limit; j += R) { 
			complex w(1, 0);
			for (int k = 0; k < mid; k++, w = w * Wn) { 
				complex x = A[j + k], y = w * A[j + mid + k]; 
				A[j + k] = x + y;
				A[j + mid + k] = x - y;
			}
		}
	}
}

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);

	string aa, bb;
	cin >> aa >> bb;
	int n = aa.size(), m = bb.size();
	for (int i = 0; i < n; i++) a[i].x = aa[i] ^ 0x30;
	for (int i = 0; i < m; i++) b[i].x = bb[i] ^ 0x30;
	while (limit <= n + m) limit <<= 1, l++;
	for (int i = 0; i < limit; i++) r[i] = ( r[i >> 1] >> 1 ) | ( (i & 1) << (l - 1));
	fft(a, 1);
	fft(b, 1);
	for (int i = 0; i <= limit; i++) a[i] = a[i] * b[i];
	fft(a, -1);
	for (int i = 0; i <= n + m; i++) ans[i] = (int)(a[i].x / limit + 0.5);
	for(int i = n + m;i >= 1;i--) ans[i - 1] += ans[i] / 10, ans[i] %= 10;
	for(int i = 0;i <= n + m - 2;i++) cout << ans[i]; cout << '\n';
}

NTT

#include <bits/stdc++.h>

using namespace std;

constexpr int P = 998244353;
vector<int> rev, roots{0, 1};

int power(int a, int b) {
	int res = 1;
	for(; b; b >>= 1, a = 1ll * a * a % P) if(b & 1) res = 1ll * res * a % P;
	return res;
}
void dft(vector<int> &a) {
	int n = a.size();
	if(int(rev.size()) != n) {
		int k = __builtin_ctz(n) - 1;
		rev.resize(n);
		for (int i = 0; i < n; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
	}
	for(int i = 0; i < n; ++i) if(rev[i] < i) swap(a[i], a[rev[i]]);
	if(int(roots.size()) < n) {
		int k = __builtin_ctz(roots.size());
		roots.resize(n);
		while((1 << k) < n) {
			int e = power(3, (P - 1) >> (k + 1));
			for(int i = 1 << (k - 1); i < (1 << k); ++i) {
				roots[2 * i] = roots[i];
				roots[2 * i + 1] = 1ll * roots[i] * e % P;
			}
			++k;
		}
	}
	for(int k = 1; k < n; k *= 2) {
		for(int i = 0; i < n; i += 2 * k) {
			for(int j = 0; j < k; ++j) {
				int u = a[i + j];
				int v = 1ll * a[i + j + k] * roots[k + j] % P;
				int x = u + v;
				if(x >= P) x -= P;
				a[i + j] = x;
				x = u - v;
				if(x < 0) x += P;
				a[i + j + k] = x;
			}
		}
	}
}
void idft(vector<int> &a) {
	int n = a.size();
	reverse(a.begin() + 1, a.end());
	dft(a);
	int inv = power(n, P - 2);
	for(int i = 0; i < n; ++i) a[i] = 1ll * a[i] * inv % P;
}
struct Poly {
	vector<int> a;
	Poly(){}
	Poly(int a0) { if (a0) a = {a0}; }
	Poly(const vector<int> &a1) : a(a1) {
		while(!a.empty() && !a.back()) a.pop_back();
	}
	int size() const { return a.size(); }
	int operator[](int idx) const { if (idx < 0 || idx >= size()) return 0; return a[idx]; }
	Poly mulxk(int k) const {
		auto b = a;
		b.insert(b.begin(), k, 0);
		return Poly(b);
	}
	Poly modxk(int k) const {
		k = min(k, size());
		return Poly(vector<int>(a.begin(), a.begin() + k));
	}
	Poly divxk(int k) const {
		if (size() <= k) return Poly();
		return Poly(vector<int>(a.begin() + k, a.end()));
	}
	friend Poly operator+(const Poly a, const Poly &b) {
		vector<int> res(max(a.size(), b.size()));
		for (int i = 0; i < int(res.size()); ++i) {
			res[i] = a[i] + b[i];
			if (res[i] >= P) res[i] -= P;
		}
		return Poly(res);
	}
	friend Poly operator-(const Poly a, const Poly &b) {
		vector<int> res(max(a.size(), b.size()));
		for (int i = 0; i < int(res.size()); ++i) {
			res[i] = a[i] - b[i];
			if (res[i] < 0) res[i] += P;
		}
		return Poly(res);
	}
	friend Poly operator*(Poly a, Poly b) {
		int sz = 1, tot = a.size() + b.size() - 1;
		while (sz < tot) sz *= 2;
		a.a.resize(sz);
		b.a.resize(sz);
		dft(a.a);
		dft(b.a);
		for (int i = 0; i < sz; ++i) a.a[i] = 1ll * a[i] * b[i] % P;
		idft(a.a);
		return Poly(a.a);
	}
	Poly &operator+=(Poly b) { return (*this) = (*this) + b; }
	Poly &operator-=(Poly b) { return (*this) = (*this) - b; }
	Poly &operator*=(Poly b) { return (*this) = (*this) * b; }
	Poly deriv() const {
		if (a.empty()) return Poly();
		vector<int> res(size() - 1);
		for (int i = 0; i < size() - 1; ++i) res[i] = 1ll * (i + 1) * a[i + 1] % P;
		return Poly(res);
	}
	Poly integr() const {
		if(a.empty()) return Poly();
		vector<int> res(size() + 1);
		for (int i = 0; i < size(); ++i) res[i + 1] = 1ll * a[i] * power(i + 1, P - 2) % P;
		return Poly(res);
	}
	Poly inv(int m) const {
		Poly x(power(a[0], P - 2));
		int k = 1;
		while(k < m) {
			k *= 2;
			x = (x * (2 - modxk(k) * x)).modxk(k);
		}
		return x.modxk(m);
	}
	Poly log(int m) const { return (deriv() * inv(m)).integr().modxk(m); }
	Poly exp(int m) const {
		Poly x(1);
		int k = 1;
		while (k < m) {
			k *= 2;
			x = (x * (1 - x.log(k) + modxk(k))).modxk(k);
		}
		return x.modxk(m);
	}
	Poly sqrt(int m) const {
		Poly x(1);
		int k = 1;
		while(k < m) {
			k *= 2;
			x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((P + 1) / 2);
		}
		return x.modxk(m);
	}
	Poly mulT(Poly b) const {
		if (b.size() == 0) return Poly();
		int n = b.size();
		reverse(b.a.begin(), b.a.end());
		return ((*this) * b).divxk(n - 1);
	}
	vector<int> eval(vector<int> x) const {
		if (size() == 0) return vector<int>(x.size(), 0);
		const int n = max(int(x.size()), size());
		vector<Poly> q(4 * n);
		vector<int> ans(x.size());
		x.resize(n);
		function<void(int, int, int)> build = [&](int p, int l, int r) {
			if (r - l == 1) {
				q[p] = vector<int>{1, (P - x[l]) % P};
			}
			else {
				int m = (l + r) / 2;
				build(2 * p, l, m);
				build(2 * p + 1, m, r);
				q[p] = q[2 * p] * q[2 * p + 1];
			}
		};
		build(1, 0, n);
		function<void(int, int, int, const Poly &)> work = [&](int p, int l, int r, const Poly &num) {
			if (r - l == 1) {
				if (l < int(ans.size()))
					ans[l] = num[0];
			} 
			else {
				int m = (l + r) / 2;
				work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
				work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
			}
		};
		work(1, 0, n, mulT(q[1].inv(n)));
		return ans;
	}
};

int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);

	int n;
	cin >> n;

	vector<int> f;
	for (int i = 0; i < n; i++) {
		int a;
		cin >> a;
		if (a >= int(f.size())) {
			f.resize(a + 1);
		}
		f[a] = 1;
	}

	int mx = f.size() - 1;

	auto g = f;
	reverse(g.begin(), g.end());  //位置取反  //差的卷积
	auto h = Poly(f) * Poly(g);   //h 数组从mx开始为差值为0,后边以此类推


	return 0;
}

FWT

OR

void fwt_or(vector<ll> &a, int len) {
    for (int mid = 2; mid <= len; mid <<= 1) {
        for (int i = 0; i < len; i += mid) {
            for (int j = i; j < i + (mid >> 1); j++) {
                a[j + (mid >> 1)] += a[j];
			}
		}
	}
	return ;
}
void ifwt_or(vector<ll> &a, int len) {
    for (int mid = 2; mid <= len; mid <<= 1) {
        for (int i = 0; i < len; i += mid) {
            for (int j = i; j < i + (mid >> 1); j++) {
                a[j + (mid >> 1)] -= a[j];
			}
		}
	}
	return ;
}

AND

void fwt_and(vector<ll> &a, int len) {
    for (int mid = 2; mid <= len; mid <<= 1) {
        for (int i = 0; i < len; i += mid) {
            for (int j = i; j < i + (mid >> 1); j++) {
                a[j] += a[j + (mid >> 1)];
			}
		}
	}
	return ;
}
void ifwt_and(vector<ll> &a, int len) {
    for (int mid = 2; mid <= len; mid <<= 1) {
        for (int i = 0; i < len; i += mid) {
            for (int j = i; j < i + (mid >> 1); j++) {
                a[j] -= a[j + (mid >> 1)];
			}
		}
	}
	return ;
}

XOR

void fwt_xor(vector<ll> &a, int len) {
    for (int mid = 2; mid <= len; mid <<= 1) {
        for (int i = 0; i < len; i += mid) {
            for (int j = i; j < i + (mid >> 1); j++) {
                ll x = a[j], y = a[j + (mid >> 1)];
                a[j] = x + y, a[j + (mid >> 1)] = x - y;
            }
		}
	}
	return ;
}
void ifwt_xor(vector<ll> &a, int len) {
    for (int mid = 2; mid <= len; mid <<= 1) {
        for (int i = 0; i < len; i += mid) {
            for (int j = i; j < i + (mid >> 1); j++) {
                ll x = a[j], y = a[j + (mid >> 1)];
                a[j] = (x + y) >> 1, a[j + (mid >> 1)] = (x - y) >> 1;
            }
		}
	}
	return ;
}

标签:return,数论,ll,Number,int,const,Mint,mod
来源: https://blog.csdn.net/u014711890/article/details/120841840

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有