Number 数论

已于 2022-03-29 14:36:02 修改
阅读量369 收藏 1
点赞数 2
分类专栏： XCPC模板 文章标签： 线性代数 算法
于 2021-10-19 10:51:13 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/u014711890/article/details/120841840
版权

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 getPrime(){
    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;
            }
        }
    }
}

任意区间素数筛

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);
    }
}

欧拉函数

​ 用处：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;
            }
		}
	}
}

Mobius

Mobius反演

  1. Mobius函数反演的第一种形式,其中f(n)为n的因数和；

f ( n ) = ∑ d ∣ n g ( d )    ⟺    g ( n ) = ∑ d ∣ n μ ( n d ) ⋅ f ( d ) f(n)=\sum_{d|n}^{}g(d){\iff}g(n)=\sum_{d|n}\mu(\frac{n}{d})·f(d) f(n)=d∣n∑​g(d)⟺g(n)=d∣n∑​μ(dn​)⋅f(d)

  1. Mobius函数反演的第二种形式,其中f(n)为gcd为n的倍数的关系求和；

f ( n ) = ∑ n ∣ m g ( m )    ⟺    g ( n ) = ∑ n ∣ m μ ( m n ) ⋅ f ( m ) f(n)=\sum_{n|m}^{}g(m){\iff}g(n)=\sum_{n|m}\mu(\frac{m}{n})·f(m) f(n)=n∣m∑​g(m)⟺g(n)=n∣m∑​μ(nm​)⋅f(m)

线性筛Mobius函数

μ ( i ) = { ( − 1 ) k p 1 p 2 ⋅ ⋅ ⋅ p k ( p i 为 素 数 ) 0 i 有 平 方 因 子 1 i = 1 \mu(i)=\left\{ \begin{matrix} (-1)^{k} && p_1p_2···p_k&(p_i为素数)\\ 0 && i有平方因子\\ 1 && i = 1 \end{matrix} \right. μ(i)=⎩⎨⎧​(−1)k01​​p1​p2​⋅⋅⋅pk​i有平方因子i=1​(pi​为素数)

constexpr int maxn = 1e5 + 7;

int prime[maxn], vis[maxn], mobius[maxn];

void getMobius() {
	mobius[1] = 1;
	for(int i = 2;i < maxn;i++) {
		if(vis[i] == 0) {
			prime[++prime[0]] = i;
			mobius[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) {
				mobius[i * prime[j]] = 0;
				break;
			}
			else {
				mobius[i * prime[j]] = -mobius[i];
			}
		}
	}
}

大数素数测试

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意义下的逆元

c o n d i t i o n : [ g c d ( a , b ) = 1 ] condition:[gcd(a, b)=1] condition:[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;
}

组合数

img

n大m小组合数-O(m)

const int mod = 1e9 + 7;

struct Binom {
	int n;
	vector<ll> inv;
	Binom(int n_ = 0) : n(n_), inv(n) {
		inv[1] = 1;
		for(int i = 2; i < n; i++) inv[i] = ((-(mod / i) * inv[mod % i]) % mod + mod) % mod;
	}
	ll C(ll n, ll m) { // m in a small range
		ll ret = 1;
		for(int i = 0; i < m; i++)
			ret *= (n - i) % mod, ret %= mod, ret *= inv[i + 1], ret %= mod;
		return ret;
	}
};

正常版本-O(logx)

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];
}

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;
}

Exlucas

容斥原理

欧几里得

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;
	}
}

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;
}

矩阵

  • 加速递推， O ( k 3 l o g 2 n ) O(k^{3}log_{2}n) O(k3log2​n)
  • 线段树无脑合并

Matrix Class

struct Matrix { 
	int n, m;
	vector<vector<ll> > a;//If you get TLE, try to instead of vector in array...

	Matrix(int n_ = 0, int m_ = 0) : n(n_), m(m_), a(n + 1, vector<ll>(m + 1, 0)) {}

	Matrix operator + (const Matrix& b) {
		Matrix ret(n, m);
		for(int i = 1;i <= n;i++) 
			for(int j = 1;j <= m;j++) 
				ret.a[i][j] = (a[i][j] + b.a[i][j]) % mod;
		return ret;
	}

	Matrix operator - (const Matrix& b) {
		Matrix ret(n, m);
		for(int i = 1;i <= n;i++) 
			for(int j = 1;j <= m;j++) 
				ret.a[i][j] = (a[i][j] - b.a[i][j] + mod) % mod;
		return ret;
	}

	Matrix operator * (const Matrix& b) {
		Matrix ret(n, b.m);
		for(int i = 1;i <= n;i++) 
			for(int j = 1;j <= b.m;j++) 
				for(int k = 1;k <= m;k++) 
					ret.a[i][j] = (ret.a[i][j] + a[i][k] * b.a[k][j]) % mod;
		return ret;
	}	
};

矩阵快速幂

constexpr int mod = 1e9 + 7;

struct Matrix { 
	int n, m;
	vector<vector<ll> > a;//If you get TLE, try to instead of vector in array...

	Matrix(int n_ = 0, int m_ = 0) : n(n_), m(m_), a(n + 1, vector<ll>(m + 1, 0)) {}

	Matrix unitMatrix() {		//to E
		Matrix E(n, m);
		for(int i = 1;i <= n;i++) E.a[i][i] = 1;
		return E;
	}
	
	Matrix f1_() {		//fib(n) = f(n - 1) + f(n - 2),  f[n] = qpow(A, n - 2);
		Matrix ret(n, m);
		ret.a[1][1] = 1, ret.a[1][2] = 1;
		ret.a[2][1] = 1, ret.a[2][2] = 0;
		return ret;
	}

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

	Matrix operator * (const Matrix& b) {
		Matrix ret(n, b.m);
		for(int i = 1;i <= n;i++) 
			for(int j = 1;j <= b.m;j++) 
				for(int k = 1;k <= m;k++) 
					ret.a[i][j] = (ret.a[i][j] + a[i][k] * b.a[k][j]) % mod;
		return ret;
	}	

	Matrix qpow(Matrix a, ll b) {
		Matrix ans = unitMatrix();
		while(b) {
			if(b & 1) ans = ans * a;
			b >>= 1,  a = a * a;
		}
		return ans;
	}
	void show() {
		for(int i = 1;i <= n;i++) 
			for(int j = 1;j <= m;j++) 
				cout << a[i][j] << " \n"[j == m];
		return ;
	}
};

生成函数

生成函数（Generating Function），又称母函数，是一种形式幂级数，其每一项的系数可以提供关于这个序列的信息。

生成函数有许多不同的种类，但大多可以表示为单一的形式

F ( x ) = ∑ n a n k n ( x ) F(x)=\sum_{n}a_{n}k_{n}(x) F(x)=n∑​an​kn​(x)

其中 k n ( x ) k_{n}(x) kn​(x)被称为核函数，不同的核函数会导出不同的生成函数，拥有不同的性质。举个例子：

  1. 普通生成函数： k n ( x ) = x n k_{n}(x)=x^{n} kn​(x)=xn。
  2. 指数生成函数： k n ( x ) = x n n ! k_{n}(x)=\frac{x^n}{n!} kn​(x)=n!xn​。

对于生成函数 F ( x ) F(x) F(x)，我们用 [ k n ( x ) ] F ( x ) [k_n(x)]F(x) [kn​(x)]F(x)来表示它的第n项核函数对应的系数，也就是 a n a_{n} an​。

普通生成函数

定义

序列 a a a的普通生成函数（Ordinary Generating Function，OGF）定义为形式幂级数：

F ( x ) = ∑ n a n k n ( x ) F(x)=\sum_{n}a_{n}k_{n}(x) F(x)=n∑​an​kn​(x)

a a a既可以是有穷序列，也可以是无穷序列。常见的例子（假设 s s s以 0 0 0为起点）:

  1. 序列 a = < 1 , 2 , 3 > a=<1,2,3> a=<1,2,3>的普通生成函数是 F ( x ) = 1 + 2 x + 3 x 2 F(x)=1+2x+3x^2 F(x)=1+2x+3x2。
  2. 序列 a = < 1 , 1 , 1.... > a=<1,1,1....> a=<1,1,1....>的普通生成函数是 F ( x ) = ∑ n ≥ 0 x n F(x)=\sum_{n\geq0}x^n F(x)=∑n≥0​xn。
  3. 序列 a = < 1 , 2 , 4 , 8... > 的 普 通 生 成 函 数 是 a=<1,2,4,8...>的普通生成函数是 a=<1,2,4,8...>的普通生成函数是 F ( x ) = ∑ n ≥ 0 2 n x n F(x)=\sum_{n\geq0}2^nx^n F(x)=∑n≥0​2nxn。
  4. 序列 a = < 1 , 3 , 5 , 7... > a=<1,3,5,7...> a=<1,3,5,7...>的普通生成函数是 F ( x ) = ∑ n ≥ 0 ( 2 n + 1 ) x n F(x)=\sum_{n\geq0}(2n+1)x^n F(x)=∑n≥0​(2n+1)xn。
封闭形式

在运用生成函数的过程中，我们不会一直使用形式幂级数的形式，而会适时地转化为封闭形式以更好地化简。

例如：序列 a = < 1 , 1 , 1.... > a=<1,1,1....> a=<1,1,1....>的普通生成函数是 F ( x ) = ∑ n ≥ 0 x n F(x)=\sum_{n\geq0}x^n F(x)=∑n≥0​xn，我们可以发现
x ⋅ F ( x ) + 1 = F ( x ) x·F(x) + 1=F(x) x⋅F(x)+1=F(x)
那么解这个方程得到
F ( x ) = 1 x − 1 F(x)=\frac{1}{x-1} F(x)=x−11​
这就是 F ( x ) = ∑ n ≥ 0 x n F(x)=\sum_{n\geq0}x^n F(x)=∑n≥0​xn的封闭形式。

( n + k − 1 ) \binom{n+k-1}{} (n+k−1​)

卷积&多项式

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

原根
多项式快速幂

F ( x ) = ( G ( x ) ) k l n ( F ( x ) ) = k ⋅ l n ( G ( x ) ) F ( x ) = e k ⋅ l n ( G ( x ) ) F(x)=(G(x))^k\newline ln(F(x))=k{·}ln(G(x))\newline F(x)=e^{k·ln(G(x))} F(x)=(G(x))kln(F(x))=k⋅ln(G(x))F(x)=ek⋅ln(G(x))

NTT(Standard)
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 mulInt(int x) const {
		auto b = a;
		for(int i = 0;i < (int)(b.size());i++) 
			b[i] = 1ll * b[i] * x % P;
		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;
}
NTT(Simplify)
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]; 
	}
	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; }
};

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 ;
}

Mess

小数循环节

int Circle(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;
}
Zjkai_
关注
  • 2
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
数论
，，，，，，，，，的博客
08-06 1475
1.同余方程：a≡b(mod m) 给定一个正整数n，如果两个整数a和b满足a-b能被n整除，即(a-b)modn=0，那么就称整数a与b对模n同余，记作a≡b(modn)，同时可成立amodn=b。 如果m是a - b的一个因子，那么就说a，b关于模m同余，即为a≡b(mod m)。 如果有a≡b(mod m)，那么b就是a模m的一个剩余。 有关性质: 传递性 同加减乘（设函数f，...
C语言的register
JCRunner的专栏
03-21 1624
摘自：http://c.biancheng.net/cpp/html/2864.html 一般情况下，变量的值是存储在内存中的，CPU 每次使用数据都要从内存中读取。如果有一些变量使用非常频繁，从内存中读取就会消耗很多时间，例如 for 循环中的增量控制： int i;for(i=0; i1000; i++){ // Some Code} 执行这段代码，CPU 为了获得
浅析C语言的一个关键字——register
鹏哥说编程
07-06 1259
转载自: http://zqwt.012.blog.163.com/blog/static/1204468420107745836163/ 浅析C语言的一个关键字——register 1、register修饰符暗示编译程序相应的变量将被频繁地使用，如果可能的话，应将其保存在CPU的寄存器中，以加快其存储速度。例如下面的内存块拷贝代码， /* Procedure for the ass
c语言 register volatile 变量,c语言之volatile和register关键字
weixin_39989980的博客
05-20 287
//register.cpp#include #include #include int main1(){for(register int i=0;i<10;i++){//由于register变量使用的是硬件CPU中的寄存器，寄存器变量无地址，所以不能使用取地址运算符"&"求寄存器变量的地址。printf("%d\n",&i);}printf("\n");return 0;}...
register int 的作用
斯坦福-MIT-Harvard-CMU
12-01 9470
使用register int申请的变量是存储在CPU中央寄存器中的（寄存器是中央处理器内的组成部分。寄存器是有限存贮容量的高速存贮部件），而使用int申请的变量是存储在内存中。 使用register修饰的变量可以提高它的读写速度，一般用于多层循环中。 使用register之后的运行 ...
register int的使用
热门推荐
LiuKairui的博客
07-19 1万+
今天看到这么一段代码 void floyd(){ for(register int a=1;a&lt;=v;a++){ for(register int b=1;b&lt;=v;b++){ for(register int c=1;c&lt;=v;c++){ if(G[b][a]+G[a][c]&lt;G[b][...
哈代数论number theory
06-27
非常经典的数论书籍
代数数论
01-02
7. **Diophantine分析**：这涉及寻找整数解的代数方程，是数论的一个分支，Lang的另一本书《Number Theory III》也涵盖了这一主题。 8. **阿贝尔函数**：这些函数在代数几何和代数数论中都有应用，它们通常与代数...
A Pathway Into Number Theory
04-03
R. P. Burn的数论入门，非常好的一本数论指导书，影印版
Register int i ( 声明)
zhangfn2011的专栏
08-25 1978
Register声明的作用 是为了提高效率。 它明确要求CPU把变量始终保存在寄存器里面，直至它消亡。 不过现代编译器都很厉害，根本不需要你多此一举。 所以根本就极少用。   补充：         大多数情况下，你声明了也没有用，因为编译器不会照你说得做，
c语言之volatile和register关键字
云守护的专栏
12-28 3654
//register.cpp #include #include #include int main1() { for(register int i=0;i<10;i++) { //由于register变量使用的是硬件CPU中的寄存器，寄存器变量无地址，所以不能使用取地址运算符"&"求寄存器变量的地址。 printf("%d\n",&i);
解题报告（二）E、（BZOJ3513） [MUTC2013] idiots（生成函数 + FFT + 组合计数）
繁凡さん的博客
03-02 369
繁凡出品的全新系列：解题报告系列 —— 超高质量算法题单，配套我写的超高质量题解和代码，题目难度不一定按照题号排序，我会在每道题后面加上题目难度指数（1∼51 \sim 51∼5），以模板题难度 111 为基准。 这样大家在学习算法的时候就可以执行这样的流程： 阅读我的【学习笔记】 / 【算法全家桶】学习算法 ⇒\Rightarrow⇒ 阅读我的相应算法的【解题报告】获得高质量题单 ⇒\Rightarrow⇒ 根据我的一句话题解的提示尝试自己解决问题 ⇒\Rightarrow⇒ 点开我的详细题解链接学.
微积分-微分应用2（平均值定理）
weixin_45911156的博客
07-20 779
平均值定理
线性代数：向量组与向量组等价
最新发布
huanfeng_AI的博客
07-21 330
向量组是由若干个同维数的行向量或列向量组成的集合。这些向量可以是二维的、三维的，甚至是更高维度的。在数学和物理学的许多领域中，向量组扮演着至关重要的角色，特别是在线性代数、矩阵论、以及工程学中。例如，一个由m个n维向量组成的向量组可以表示为A: a₁, a₂, …, aₘ，其中每个aᵢ都是一个n维向量。当m > n时，这样的向量组往往是线性相关的，即存在不全为零的系数k₁, k₂, …, kₘ，使得k₁a₁ + k₂a₂ + … + kₘaₘ = 0。
【18】输出字符菱形
2301_78786609的博客
07-19 372
用‘*’构造一个对角线长n个字符，倾斜放置的菱形。
Leetcode—207. 课程表【中等】
Atticus的博客
07-19 310
运行结果 之后我会持续更新，如果喜欢我的文章，请记得一键三连哦，点赞关注收藏，你的每一个赞每一份关注每一次收藏都将是我前进路上的无限动力 ！！！↖(▔▽▔)↗感谢支持！
编程视角的初等数论
"Elementary Number Theory with Programming" 是一本适合对数论和编程都感兴趣的读者的教材，无论他们是初学者还是有一定基础的程序员或数学爱好者。通过阅读这本书，读者可以增强自己的数学思维，同时提升编程技能...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值