其他分享
首页 > 其他分享> > Number 数论

Number 数论

作者:互联网

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