算法数学笔记-一、数论

目录

• 一、初等数论
  • 素数与约数的一些常用结论
  • 扩展欧几里得与裴蜀定理
    • • 扩展欧几里得
      • 裴蜀定理
  • 扩展欧拉定理与费马小定理
    • • 费马小定理
      • 扩展欧拉定理
  • 逆元
    • • 扩展欧几里得
      • 线型推法
      • 求任意n个数的逆元
  • 阶、原根、二次剩余
    • • 阶：
      • 原根
      • 二次剩余
  • 类欧几里得算法
  • 扩展中国剩余定理
  • 扩展Lucas
  • 威尔逊定理及推论
  • 升幂定理
  • 扩展BSGS
  • 狄利克雷卷积计算
  • 莫比乌斯反演，欧拉反演 与 狄利克雷卷积
  • 狄利克雷前缀和
    • • (倒)狄利克雷前缀和
      • (倒)狄利克雷后缀和
  • 杜教筛
  • 积性函数在组合数处的取值
    • • 莫队算法\(O(N ^{3 / 2})\)
      • 一阶判别式\(O(\frac{N^{3/2}}{\ln n})\)
  • PN筛
  • min25筛
  • Miller Rabin素数测试 & Pollard Rho因数分解

一、初等数论

素数与约数的一些常用结论

素数数量，对1~n，约有素数 π(n) ≈ n/ln n

在 N! 中质因子p的个数为\(\lfloor \frac{N}{p^1} \rfloor + \lfloor \frac{N}{p^2} \rfloor + \lfloor \frac{N}{p^3} \rfloor + \cdots \cdots + \lfloor \frac{N}{p^{\lfloor log_p N \rfloor}} \rfloor\)

\(N = p_1^{c_1} ·p_2^{c_2}·p_2^{c_2}\cdots p_k^{c_k}\)

N 的正约数的个数为 \((c_1 + 1)(c_2 + 1)(c_3 + 1) \cdots \cdots (c_k + 1)\)

N的正约数的和为 \((1 + p_1^1 + p_1^2 \cdots p_1^{c_1}) (1 + p_2^1 + p_2^2 \cdots p_2^{c_2})\cdots\cdots(1 + p_k^1 + p_k^2 + p_k ^ 3 \cdots p_k^{c_k})\)

扩展欧几里得与裴蜀定理

扩展欧几里得

解同余方程ax + by = gcd(a, b)

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

求aX + bY = s的通解:

exgcd(a, b, x, y);

\(X = s / gcd * x + k_1 * (b / gcd); (k_1 \in Z)\)

\(Y = s / gcd * y + k_2 * (a / gcd); (k_2 \in Z)\)

裴蜀定理

设a，b是不全为零的整数，若存在ax + by = c，则必有 gcd(a, b) | c

可扩展至多个数字（eg. ax + by + cz = d 必有 gcd(a, b, c) | d）

扩展欧拉定理与费马小定理

费马小定理

​ \(若p为质数，则有a^p \equiv a \ (mod\ p)\)

扩展欧拉定理

\[a^b \equiv \begin{cases} a^{b \ mod \ \varphi(m)},&gcd(a, m) = 1 \\ a^{b \ mod \ \varphi (m) \ +\ \varphi(m)},&gcd(a, m) \ne 1 \end{cases} \ ( mod \ m)\ \hspace{9cm} \\ \\满足a，m互质,则满足a^x \equiv1\ (mod\ m)的最小正整数x_0,则有x_0\ |\ \varphi(m)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \]

逆元

扩展欧几里得

\(\ \ \ \ \ ax \equiv 1 (mod \ m) \implies ax +my = 1, \ \ \ \ \ \ a^{-1} = x\)

线型推法

inv[1] = 1;
inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;

求任意n个数的逆元

计算前缀积s[i], 再计算s[n]的逆元sinv[n], sinv[i] = sinv[i + 1] * a[i + 1], inv[i] = sinv[i] * s[i - 1];

vec GetInvSum(vec a){
	
	int n = a.size() - 1;
	vec s(n + 1), sinv(n + 1);
	
	s[0] = a[0];
	for(int i = 1; i <= n; i ++)
		s[i] = 1ll * s[i - 1] * a[i] % mod;
	sinv[n] = GetInv(s[n]);
	
	for(int i = n - 1; i >= 0; i --)
		sinv[i] = 1ll * sinv[i + 1] * a[i + 1] % mod;
	
	for(int i = 1; i <= n; i ++)
		sinv[i] = 1ll * sinv[i] * s[i - 1] % mod;
		
	return sinv;
}

阶、原根、二次剩余

阶：

因此满足\(a^n \equiv 1 \ (mod \ m)\)同余式 的最小正整数n存在，这个n称作a模m的阶，记作\(\delta_m(a) = n\)。

性质：

\(a^1, a^2 , a^3 \ \cdots \ a^{\delta_m(a)} \ \ (mod \ m)\)两两互不同余

若\(a^m \equiv 1 \ ( mod \ m)\)则\(\delta_m(a)\ | \ n\); 若\(a^p \equiv a^q \ (mod \ m)\)则有\(p \equiv q \ (mod \ \delta_m(a))\)

！！！WA： 若\(gcd(a, m) == gcd(b, m) == 1\)，则\(\delta_m(ab) = \delta_m(a) \ * \ \delta_m(b)\)

原根

\(m \in N^+, a \in Z, gcd(a, m) == 1, 且 \delta_m(a) == \varphi(m), 则a为模m的原根\)

原根的判定：\(m >= 3, gcd(a, m)==1, 若对 \varphi(p)的每个素因数p，都有a^{\frac{\varphi(m)}{p}} \not\equiv 1(\%m), 则a是模m的原根\)

原根的个数：\(若m有原根，则有\varphi(\varphi(m))个原根\)

原根存在定理：\(m当且仅当m=2,4,p^{\alpha},2p^{\alpha}时有原根(p为奇素数, \alpha \in N^+)\)

最小原根数量级：\(若m 有原根，则其最小原根不超过m^{\frac{1}{4}}级别\)

模2的原根为1.

\(若g为\%m的最小原根，则\% m的其它原根由g^k \%m得到（gcd(k, \varphi(m) == 1)）\)

\(gcd(a, p) == gcd(b, p) == 1,则 存在c\in Z,使\delta_p(c) = lcm(\delta(a), \delta(b))\)

\(存在g为 \% p的原根，使g^{p-1} \not\equiv 1 (\% p^2)\)

inline long long getG(long long x){

	static long long q[10051];
	long long top = 0, flg = 1;
	for(register int i = 2; i * i <= x - 1; i ++){
		if((x - 1) % i == 0){
			q[++ top] = i;
			if(i * i != x - 1)
				q[++ top] = (x - 1) / i;
		}
	}	
	
	for(register int i = 2; ; i ++){
		flg = 1;
		for(register int j = 1; j <= top; j ++)
			if(quickpow(i, q[j], x) == 1){
				flg = 0;
				break;
			}
		if(flg)	return i;
	}
	

二次剩余

namespace Shengyu2{

	LL I, p, mod;

	struct Complex2{ // 此complex为二次剩余专用 
		LL real, fals;
		
		Complex2(double a = 0, double b = 0){
			real = a, fals = b;
		}
		
		Complex2 operator * (const Complex2 & x){
			LL a, b;
			a = (real * x.real + fals * x.fals % mod * I) % mod;
			b = (real * x.fals + fals * x.real) % mod;
			return Complex2(a > 0 ? a : a + mod, b > 0 ? b : b + mod);
		}
		
		Complex2 operator % (const LL & mod){
			real = real % mod;
			fals  = fals  % mod;
			return Complex2(real > 0 ? real : real + mod, fals > 0 ? fals : fals + mod);			
		}	
	};	
		
	Complex2 quickpow(Complex2 x, LL b, LL mod){
	
		if(b == 0) return Complex2(1, 0);
		if(b == 1) return x % mod;
		
		Complex2 l = quickpow(x, b >> 1, mod);
		l = (l * l) % mod;
		if(b & 1)
		l = (l * x) % mod;
		return l;
	}
	
	bool IsShengyu2(LL x, LL p){
		return (::quickpow(x, (p - 1) >> 1, p) == 1) || (x == 0);
	}

	pair<LL, LL> GetShengyu2(long long n, LL MOD){ // ling : get n 的二次剩余 
		p = mod = MOD;
		
		if(n == 0) return make_pair(0, 0);
		if(IsShengyu2(n, p) == false) 
			return make_pair(-1, -1); // (-1,-1) 表示无解 
		
		srand(time(0));
		long long a = rand();
		
		while(IsShengyu2((a * a - n + mod) % mod, p)) a = rand();
		
		Complex2 ans(a, 1);
		I = ((a * a - n) % mod + mod) % mod;
		
		ans = quickpow(ans, (p + 1) >> 1, mod);
		LL ans1 = ans.real, ans2 = mod - ans.real;
		// 理论上n的二次剩余有两个, 有时候会ans1 == ans2 
		return make_pair(min(ans1, ans2), max(ans1, ans2));
	}
}

类欧几里得算法

常用结论：

\(m = \lfloor \frac{an +b}{c} \rfloor \\ f(a, b, c, n) = \sum^n_{i=0} \lfloor \frac{ai + b}{c} \rfloor \\g(a, b, c, n) = \sum^n_{i = 0} \lfloor \frac{ai + b}{c} \rfloor ^ 2 \\ h(a, b,c,d) = \sum^n_{i = 0} i\lfloor \frac{ai + b}{c} \rfloor\)

\[f = \begin{cases} (n + 1) \lfloor \frac{b}{c} \rfloor & a = 0 \\ \frac {n(n+ 1)}{2} \lfloor \frac{a}{c} \rfloor + (n + 1) \lfloor \frac{b}{c} \rfloor + f(a \% c , b\% c, c , n) & max(a, b) \geq c \\ nm - f(c, c - b - 1, a, m - 1) & otherwise \end{cases} \]

\[g= \begin{cases} (n + 1) \lfloor \frac{b}{c} \rfloor ^ 2 & a = 0 \\ g(a \%c, b \% c, c, n) + 2 \lfloor \frac{a}{b} \rfloor h(a \% c, b\%c, c, n) + 2\lfloor \frac{b}{c} \rfloor f(a \%c, b\%c, c, n) \\ + \frac{n(n + 1)(2n + 1)}{6}\lfloor \frac{a}{c} \rfloor ^2+ n(n + 1)\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c} \rfloor + (n+1) \lfloor \frac{b}{c} \rfloor ^ 2 & max(a, c) \geq c \\ nm(m + 1) - f(a, b, c, n) - 2h(c, c - b - 1, a, m - 1) - 2f(c, c - b - 1, a, m - 1) & otherwise \end{cases} \]

\[h = \begin{cases} \frac{n(n + 1)}{2} \lfloor \frac{b}{c} \rfloor & a = 0 \\ h(a \% c, b \%c , c, n) + \frac{n(n + 1)(2n + 1)}{6}\lfloor \frac{a}{c} \rfloor + \frac{n(n + 1)}{2} \lfloor \frac{b}{c} \rfloor & max(a, b) \geq c \\ \frac{1}{2} [mn(n + 1) - g(c, c - b - 1, a, m - 1) - f(c, c - b - 1, a, m - 1)] & ortherwise \end{cases} \]

扩展中国剩余定理

\[求解\ \ \ x\equiv ai_i \ (mod \ mi_i)\ \hspace{13cm} \]

long long exCRT(int n, int mi[], int ai[]){
	for(int x, y, i = 2; i <= n; i ++){
		long long gcd = exGcd(mi[i], -mi[i - 1], x, y);
		long long lcm = abs(mi[i] / gcd * mi[i - 1]);
		x = quickmul(x, (ai[i - 1] - ai[i]) / gcd, lcm);
		ai[i] = ((quickmul(mi[i], x, lcm) + ai[i]) % lcm + lcm) % lcm;
		mi[i] = lcm;
	}
	
	return ai[n];
}

扩展Lucas

long long Vp_fact(long long n, long long q){ // 1~n中 质因子p 的个数
	long long re = 0;
	for(long long i = q; i <= n; i *= q)
		re += n / i;
		
	return re;
}

long long Askfact(long long n, long long q, long long mod){// 具体这个函数是干啥的我也不知道
	if(n <= q)
		return fact[n];
	else return Askfact(n / q, q, mod) * quickpow(fact[mod], n / mod, mod) % mod * fact[n % mod] % mod;
    // 此处的fact[i] 表示 1~i 中除了q的倍数之外的数的积
}
long long inv(long long x, long long mod){
	long long re = 0, k = 0;
	exGcd(x, mod, re, k);
	return (re % mod + mod) % mod;
}

long long multi_Lucas(long long n, long long m, long long mod, long long q, long long k){
	
	fact[0] = 1;
	for(long long i = 1; i <= mod; i ++)
		if(i % q == 0)
			fact[i] = fact[i - 1];
		else fact[i] = fact[i - 1] * i % mod;
    
	long long x = Vp_fact(n, q), y = Vp_fact(m, q), z = Vp_fact(n - m, q);
	long long XX = Askfact(n, q, mod), YY = Askfact(m, q, mod), ZZ = Askfact(n - m, q, mod);
	
	return XX * inv(YY, mod) % mod * inv(ZZ, mod) % mod * quickpow(q, x - y - z, mod) % mod;
}// ! ! 这里的inv一定要用exGcd()求解, 因为mod不一定是质数

long long exLucas(long long n, long long m, long long mod){
	
	long long P = mod;
	long long q[33], c[33], cnt = 0;
	long long mi[33], ai[33];
	
	for(long long i = 2; i * i <= P; i ++)
		if(P % i == 0){
			q[++ cnt] = i; c[cnt] = 0;
			while(P % i == 0)
				c[cnt] ++, P /= i;
		}
	if(P > 1) q[++ cnt] = P, c[cnt] = 1;
	
	for(long long i = 1; i <= cnt; i ++)
		mi[i] = quickpow(q[i], c[i]), ai[i] = multi_Lucas(n, m, mi[i], q[i], c[i]);
		
	return exCRT(cnt, mi, ai);
}

威尔逊定理及推论

对于素数p， 有\(\ (p - 1)! \equiv -1 \ (mod \ p)\)

对于整数n，令\((n !)_p\)表示所有小于等于 n 但不能被 p 整除的正整数的乘积

\[\begin{array}{l} (n!)_p = n! / (\lfloor n / p \rfloor ! \ p^{\lfloor n / p \rfloor} \ \ \ (p为奇素数)) \hspace{17cm} \\ (p!)_p = (p - 1) ! \equiv -1\ (mod \ p)\\ \\ (p^q ! )_p = \begin{cases} 1, & if(p = 2\ \&\& \ q \geq 3) \\ -1& otherwise &\end{cases} (mod \ p^q) \end{array} \]

推论

对于素数 p 和正整数 q 和非负整数 n 有

\((n!)_p \ \equiv (\ (p^q!)_p\ )^{\lfloor n / p^q \rfloor} (\ (n\% p^q)!\ )_p \ \ (mod \ p^q)\)

若存在\(1 \leq r \leq p - 1 \ 使(-1)^r r! \equiv 1\ (mod \ p)则(p - r - 1 )! + 1 \equiv 0 \ (mod \ p)\)

升幂定理

记\(v_p(n)\)为正整数 n 中 p 的个数， 即\(p^{v_p(n)}\) 恰好正整除 n

1. p 为奇素数
   \(a \equiv b \ (mod \ p)\)且 a、b 不是 p 的倍数，\(v_p(a^n - b^n) = v_p(a - b) + v_p(n)\)

2. p = 2

   a、 b 为素数

   n为奇数时\(v_2(a^b - b^n) = v_2(a - b)\)

   n为偶数时\(v_2(a^n - b^n) = v_2(a - b) + v_2(a + b) + v_2(n) - 1\)

扩展BSGS

long long exBSGS(long long a, long long b, long long p){
	
	map<long long ,long long > Map_bsgs;
	
	a %= p, b %= p;
	if(b == 1 || p == 1) return 0;
	
	long long cnt = 0, add = 1;
	for(long long d = Gcd(a, p); (d ^ 1); d = Gcd(a, p)){
		if(b % d) return -1;
		
		b /= d, p /= d, cnt ++;
		add = (a / d) * add % p;
		if(add == b) return cnt;
	}
	
	Map_bsgs.clear();
	long long t = ceil(sqrt(p));
	long long a_t = 1; 
	for(long long i = 0; i <= t - 1; i ++, a_t = a_t * a % p)
		Map_bsgs[b * a_t % p] = i;
	for(long long i = 0, now = add; i <= t; i ++, now = now * a_t % p)
		if(Map_bsgs.find(now) != Map_bsgs.end())
			if(i * t - Map_bsgs[now] >= 0) 
				return i * t - Map_bsgs[now] + cnt;
			
	return -1;
	
}

狄利克雷卷积计算

在差不多线型的复杂度求狄利克雷卷积\(A_n = \sum_{d|n}f(d)g(\frac{n}{d})\) （f 和 g 都是数论积性函数）

实际时间消耗大概是线筛的两倍，1.5e7/s

#define f(x)
#define g(x) 

g(1) = 1;
f(1) = 1;

for(int i = 2; i <= n; i ++){
		if(vis[i] == 0){
			lpn[i] = i; prime[++ cnt] = i;
			A[i] = (g(i) + f(i)) % mod;
		}
		
		for(int j = 1; j <= cnt && prime[j] * i <= n; j ++){
			vis[prime[j] * i] = 1;
			
			if(i % prime[j] == 0){
			
				int x = i * prime[j], p = prime[j];
				lpn[x] = lpn[i] * p;
	
				if(lpn[x] != i * prime[j])
					A[x] = 1ll * A[lpn[x]] * A[x / lpn[x]] % mod;
				else{
					A[i * prime[j]] = 0;	
					// 注意 ↓ 这里的 LL 不要写成 int， 如果可以的话全开LL比较方便
					for(long long now = 1; now <= lpn[x]; now *= p)
						A[i * p] = (A[i * p] + 1ll * f(now) * g(lpn[x] / now)) % mod;	
					A[i * p] = 1ll * A[i * p] * A[i * p / lpn[x]] % mod;			
				}
				break;	
			}
			else
				lpn[i * prime[j]] = prime[j], A[i * prime[j]] = 1ll * A[i] * A[prime[j]] % mod;
		}
	}

莫比乌斯反演，欧拉反演 与 狄利克雷卷积

莫反结论

\[[gcd(i, j)= 1] = \sum_{d|i\and d|j} \mu(d)\\ \]

\[gcd(i, j) = \sum_{d|i \and d|j} \varphi(d)\\ \]

\[f(n) = \sum_{d|n}g(d) \ \ \ \ \ g(n) = \sum_{d|n} \mu(d)f(\frac{n}{d})\\ \]

\[f(d) = \sum_{d|n}g(n) \ \ \ \ \ g(d) = \sum_{d|n} \mu(\frac{n}{d})f(n)\\ \]

\[\sum_{j=1}^n gcd(i,j)=\sum_{j=1}^n\sum_{d|i\and d|j}\varphi(d) = \sum_{d|i}\varphi(d)\lfloor \frac n d\rfloor \\ \]

\[\sum_{j=1}^n gcd(0,j)=\sum_{d=1}^n\varphi(d)\lfloor \frac n d\rfloor = \sum_ {j=1}^nj =\frac{n*(n+1)}2 \\ \]

\[\sum_{i=1}^n \sum_{j=1}^m gcd(i, j) =\sum_{i=1}^n \sum_{j=1}^m \sum_{d|i\and |j} \varphi(d) = \sum_{d = 1}^{min(n,m)} \varphi(d) \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \\ \]

\[\sum_{i=1}^n\sum_{j=1}^mgcd(i,j,k)=\sum_{i=1}^n\sum_{j=1}^m\sum_{d|i\and d|j\and d|k}\varphi(d)=\sum_{d|k}\varphi(d) \lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor \\ \]

\[\sum_{i = 1} ^ {n}\sum_{j = 1}^m [ gcd(i,j) = 1 ] = \sum_{i = 1} ^ {n}\sum_{j = 1}^m \sum_{x|i\and x|j}\mu(x) = \sum_{x = 1}^{min(n,m)} \mu(x) \lfloor \frac{n}{x}\rfloor \lfloor \frac{m}{x} \rfloor \\ \]

\[\sum_{k=1}^n[gcd(i,k)=d] = \sum_{k=1}^{n/d}[gcd(\frac i d,k)=1]=\sum_{k=1}^{n/d}\sum_{g|\frac i d}[g|k] \mu(g) = \sum_{g|\frac i d}\mu(g)\lfloor\frac{n/d} g \rfloor \]

\[\sum_{i = 1} ^ {n}\sum_{j = 1}^m [ gcd(i,j) = d ] = \sum_{i = 1} ^ {n/d}\sum_{j = 1}^{m/d} \sum_{x|i\and x|j}\mu(x) = \sum_{x = 1}^{min(n,m)/d} \mu(x) \lfloor \frac{n}{dx}\rfloor \lfloor \frac{m}{dx} \rfloor = \sum_{d|T}^{min(n,m)}\mu(\frac{T}{d})\lfloor \frac{n}{T}\rfloor\lfloor \frac{m}{T}\rfloor\\ \]

\[\prod_{i=1}^n\prod_{j=1}^mgcd(i,j)=\prod_{d=1}^{min(n,m)}d^{\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=d]}=\prod_{d=1}^{min(n,m)}\prod_{d|T}d^{\mu(\frac{T}{d})\lfloor \frac{n}{T}\rfloor\lfloor \frac{m}{T}\rfloor} \]

欧拉反演

\[n = \sum _{d|n} \varphi(d) \\ \sum _ {i = 1} ^ n gcd(i, n) = \sum _{d|n} \frac{n}{d}\varphi(d) \]

常用数论函数

狄利克雷卷积单位元:\(\epsilon (n) = [n == 1]\)

$ \ Id_k(n) = n^k$

恒等函数:$Id(n) = n $

除数函数:\(\sigma_k (n) = \sum _{d|n} d^k\)； 当k == 1时，即为除数和 函数\(\sigma(n)\)；当k==0时，即为除数个数函数\(\sigma_0(n) = d(n)\)

常值函数:\(I(n) = 1(n) = 1\)

\[1(n) * Id_k(n) = \sigma_k(n) \\ 1(n) * \varphi(n) = \sum_{d|n} \varphi(d) = n = Id(n) \\ \mu(n) * Id(n) = \varphi(n) \\ d(ij) = \sum_{x|i} \sum_{y|j}[gcd(x, y) == 1] \\ \mu(n) * 1(n) = \epsilon(n) \]

逆元

\(f(n)* g(n) = \epsilon(n)\)，若f是积性函数，则\(g = f^{-1}\)也是积性函数

\[f^{-1} = \begin{cases} \frac{1}{f(n)} & n = 1 \\ -\frac{1}{f(1)}\sum_{d|n,d\not = 1} f(d)f^{-1}(\frac{n}{d}) & ortherwise \end{cases} \]

狄利克雷前缀和

复杂度O(n log log n)

(倒)狄利克雷前缀和

\[b[n] = \sum _{d | n}a[d] \]

b[i] = a[i];// 已知 a 求 b
for(int i = 1; i <= cnt; i ++)
	for(int j = 1; j * Prime[i] <= N; j ++)
		b[j * Prime[i]] += b[j];

a[i] = b[i]; // 已知 b 求 a
for(int i = cnt; i >= 1; i --)
	for(int j = N / Prime[i]; j >= 1; j --)
		a[j * Prime[i]] -= a[j];

(倒)狄利克雷后缀和

\[b[d] = \sum_{d | n} a[n] \]

b[i] = a[i]; // 已知 a 求 b
for(int i = 1; i <= cnt; i ++)
    for(int j = N / Prime[i]; j >= 1; j --)
            b[j] += b[j * Prime[i]];

a[i] = b[i]; // 已知 b 求 a
for(int i = cnt; i >=1; i --)
	for(int j = 1; j * Prime[i] <= N; j ++)
		a[j] -= a[j * Prime[i]];

杜教筛

\(设S(n) = \sum _{i=1}^n f(i)\)

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

积性函数在组合数处的取值

定义\(v_p(n) = \sum_{k >= 1} \lfloor \frac{n}{p_k} \rfloor\) 即\(v_p(n)\) 是1~n中质因子p的个数, 记\(f_p(\alpha) = f(p^{\alpha})\)于是有

\[f( C_{n}^{m}) = \prod_{p \in P} f_p(v_p(n) - v_p(m) - v_p(n - m)) \]

其中\(v_p(n) - v_p(m) - v_p(n - m) = \sum _{k >= 1} (\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor)\)

又因\(\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor \in \{0, 1\}\)

所以 \(v_p(n) - v_p(m) - v_p(n - m) <= log_pn\)

\(p1 : p <= \sqrt n\)

\(p2 : p > \sqrt n\)

莫队算法\(O(N ^{3 / 2})\)

1.\(开一个数组维护C_n^m中个p2数 (O(n\sqrt n))\)

2.\(另外预处理p1的个数 (O(\frac{\sqrt n}{\ln n} Q))\)

莫队维护过程中：

每次n，m的变化，对于新的n,m，只需修改n,m,n- m这三个数的质因子\(p2\)的贡献，显然某个数最多有1个\(p2\)所以是\(O(1)\)的

一阶判别式\(O(\frac{N^{3/2}}{\ln n})\)

1.\(对于p1 预处理 O(\frac{\sqrt n}{\ln n} Q)\)

2.p2:

因为：

\(\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor \in \{0, 1\}\)

\([\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor = 1] * f_p(1) = f_p(1)^{\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor}\)

\([\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor = 0] * f_p(0) = f(1) = 1 = f_p(1)^{\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor}\)

所以：\(\prod_{p2} f_p(\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor) \\ = \prod _{p2} f_p(1)^{\lfloor\frac{n}{p_k}\rfloor - \lfloor\frac{m}{p_k}\rfloor - \lfloor\frac{n - m}{p_k}\rfloor}\)

令\(S(n) = \prod _{p2} f_p(1) ^{\lfloor \frac n p \rfloor}\)

于是\(p2\)的贡献既是\(\frac {S(n)}{S(m)S(n - m)}\)

S预处理：每次只对(0 / 1)个p2 修改\(O(n)\)

PN筛

PN筛求积性函数f的前缀和：

构造一个g函数

​ 1.g是积性函数

​ 2.g的前缀和易求

​ 3.对于质数p，g(p) == f(p)

#include<bits/stdc++.h>

using namespace std;

#define LLL __int128
 
const LLL mod = 4179340454199820289, W = 5e6;

LLL sum_g[W + 1], prime[W + 1], Inv[111];
bool vis[W + 1];
LLL sum, n, m, T, cnt, tot, tot_PN, ans, inv6, inv2;
int yy;

map<LLL , LLL >Map_sum_g;
vector <LLL > h[W];

LLL exGcd(LLL a, LLL b, LLL &x, LLL &y){
	if(b == 0) {x = 1, y = 0; return a; };
	LLL d = exGcd(b, a % b, x, y);
	LLL z = x; x = y; y = z - y * (a / b);
	return d;
}

LLL inv(LLL x, LLL mod){
	LLL re = 0, k = 0;
	exGcd(x, mod, re, k);
	return (re % mod + mod) % mod;
}

void getprime(LLL n){
	
	tot = 0;
	for(LLL i = 2; i <= n; i ++){
		if(vis[i] == false)
			prime[++ tot] = i;
		for(LLL j = 1; j <= tot && prime[j] * i <= n; j ++){
			vis[i * prime[j]] = true;
			if(i % prime[j] == 0) break;
		}
	}	
}

bool ok(LLL x){
	return 0 < x && x <= n;
}

LLL Getf(LLL p, LLL c, LLL pc){
	return pc * Inv[c] % mod;
}

LLL Getg(LLL p, LLL c, LLL pc){
	return pc % mod;
}

LLL GetG(LLL n, LLL d, LLL x){
		return x * (x + 1) / 2 % mod; 
}

void dfsPN(LLL c, LLL now, LLL hh){
	
	ans = (ans + hh * GetG(n, now, n / now)) % mod;	
	for(LLL i = c + 1, p = prime[c + 1]; i <= tot && ok(now * p) && ok(now * p * p); i ++, p = prime[i])
		for(LLL k = 2, pk = p * p; ok(now * pk); k ++, pk *= p)
			dfsPN(i, now * pk, hh * h[p][k] % mod);
	
}

inline LLL quickpow(LLL x, long long b, LLL mod){
	LLL j = 1;
	; 
	for(; b; b >>= 1, x = x * x % mod)
		if(b & 1) 
			j = j * x % mod;
			
	return j % mod;
}

void work(){
	
	long long t;
	cin >> t;
	n = t;
	
	ans = 0; 
	dfsPN(0, 1, 1);

	long long ANS = ans * quickpow(n, mod - 2, mod) % mod % mod;
	cout << ANS << endl;
	
}

signed main(){
	
	n = 1e12 + 1;
	
	Inv[1] = 1;
	for(int i = 2; i <= 100; i ++)
		Inv[i] = (mod - mod / i) * Inv[mod % i] % mod;
	
	inv6 = inv(6, mod);
	inv2 = inv(2, mod);
	getprime(W - 1);
	
	for(LLL u = 1; u <= tot && prime[u] * prime[u] <= n; u ++){
		LLL p = prime[u];
		h[p].push_back(1);
		for(LLL c = 1, p_c = p; p_c <= n; p_c *= p, c ++){
			h[p].push_back(Getf(p, c, p_c));
			for(LLL i = 1, p_i = p; i <= c; i ++, p_i *= p)
				h[p][c] = ((h[p][c] - Getg(p, i, p_i) * h[p][c - i]) % mod + mod) % mod;
		}
	}	
	
	int T = 1;
	cin >> T;
	while(T --) work();
}

min25筛

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

#define int long long

int quickpow(int x, int b, int mod);

const int N = 2e6+100, mod = 1e9+7, inv6 = quickpow(6, mod - 2, mod), inv2 = quickpow(2, mod - 2, mod);

int Q, tot_prime, n, sw, ans = 0, lenth = 2, w[N];
int g[11][N], id1[N], id2[N], gp[11][N], prime[N], vis[N], K[11], c[11], G[N], F_prime[N];
int kk;

inline int getid(int x){
	return (x <= Q ? id1[x] : id2[n / x]);
}

int quickpow(int x, int b, int mod){

	if(b == 0) return 1;
	if(b == 1) return x % mod;

	int l = quickpow(x, b >> 1, mod);
	l = l * l % mod;
	if(b & 1) return l * x % mod;
	else return l;
}

inline int f_p(int p){//真实的f(p)的表达式 
// 注意取模 
	p %= mod;
	return p * (p - 1) % mod;
}
 
int f_p_k(int p_k, int p, int k){//真实的f(p_k)的表达式 
// 注意取模 
	p_k %= mod;
	if(k == 1) return f_p(p);
	else return p_k * (p_k - 1) % mod;
}

int S(int x, int y){
// S(x, y) : 1 ~ x 中 minp > prime[y]的数字的多项式的和 
	if(prime[y] >= x) return 0;

	int p = getid(x);
	int re = (G[p] - F_prime[y] + mod) % mod;

	for(int i = y + 1; i <= tot_prime && prime[i] * prime[i] <= x; i ++){
		int p_k = prime[i];

		for(int k = 1; p_k <= x; k ++, p_k *= prime[i]){
			int o = p_k % mod;
			int add = 0;
			add = f_p_k(p_k, prime[i], k);
			re = (re + add * (S(x / p_k, i) + (k != 1))) % mod;
		}
	}
	return re;
}

int sum_mi(int n, int k, int mod){
	if(k == 0) return n % mod;
	else if(k == 1) return(n + 1) * n % mod * inv2 % mod;
	else if(k == 2) return(n + 1) * n % mod * (2 * n + 1) % mod * inv6 % mod;
	else if(k == 3) return quickpow(n % mod * (n + 1) % mod * inv2 % mod, 2, mod);
}

void work(){

	int ans = 0;
	cin >> n; 
	Q = sqrt(n);
	
	lenth = 2;
	c[1] = 2, K[1] = 1;
	c[2] = 1, K[2] = -1;
	
	//这里的是f(p)（不是f(p_k)）的多项式，但只用于预处理，与真实的f(p)无直接关系 
	
	
	{// clear
		tot_prime = sw = 0;
		for(int i = 1; i <= Q; i ++)
			vis[i] = 0;
		for(int i = 1; i <= 3 * Q; i ++){
			for(int j = 1; j <= lenth; j ++)
				g[j][i] = 0;	
			G[i] = 0;		
		}
	}
	
	
	prime[0] = 1;
	for(int i = 2; i <= Q; i ++){
		if(!vis[i]) prime[++ tot_prime] = i;
		for(int j = 1; j <= tot_prime && i * prime[j] <= Q; j ++){
			vis[i * prime[j]] = 1;
			if(i % prime[j] == 0) break;
		}
	}


	for(int i = 1; i <= tot_prime; i ++)
		for(int j = 1; j <= lenth; j ++)
			gp[j][i] = (gp[j][i - 1] + quickpow(prime[i], c[j], mod)) % mod;
	
	// gp[j][i] 表示prime[1] ~ prime[i] 的单项式的和 i <= tot_prime (prime <= sqrt(n)) 


	for(int l = 1, r; l <= n; l = r + 1){
		r = n / (n / l);

		w[++ sw] = n / r;
		int x = w[sw] % mod;
		for(int i = 1; i <= lenth; i ++)
			g[i][sw] = (sum_mi(x, c[i], mod) - 1 + mod) % mod;

		if(n / r <= Q)
			id1[n / r] = sw;
		else id2[r] = sw;
	}
	
	
/*
	查询1 ~ n / x的单项式的和 --> g[j][sw]
	if(n / x <= Q) : sw = id1[n / x];
	if(n / x >  Q) : sw = id2[x]; 
	w[sw] = n / x

	y = n / x
	查询1 ~ y的单项式的和 --> g[j][getid(y)]
	w[getid(y)] = y
*/ 

	
	for(int i = 1; i <= tot_prime; i ++)
		for(int j = 1, squ = prime[i] * prime[i]; j <= sw && w[j] >= squ; j ++){
			int p = w[j] / prime[i];
			p = getid(p);

			for(int v = 1; v <= lenth; v ++)
				g[v][j] = (g[v][j] - quickpow(prime[i], c[v], mod) * (g[v][p] - gp[v][i - 1] + mod) % mod + mod) % mod;
	}
	
	// 查询1 ~ n / x的所有 质数 的单项式的和 --> g[j][sw]
	 
	 
	for(int j = 1; j <= sw; j ++)
		for(int i = 1; i <= lenth; i ++)
			G[j] = (G[j] + K[i] * g[i][j]) % mod;
			
	// 询1 ~ n / x的所有 质数 的多项式的和 --> G[sw]	
	
/*	
	for(int j = 1; j <= sw - 1; j ++)
	G[j] = (G[j] + 2) % mod;
	这里是G[n / x]中f(p)由多项式之和 向 真实表达式之和 的转化 
*/
	for(int i = 1; i <= tot_prime; i ++)
		F_prime[i] = (F_prime[i - 1] + f_p(prime[i])) % mod;
	//F_prime[i] 小于等于prime[i]的质数的真实表达式之和 , prime[i] <= Q
	
	
	cout << ((S(n, 0) + 1) % mod + mod) % mod;

	return ;
}

signed main(){
	
	work();
}
// 100000000000

Miller Rabin素数测试 & Pollard Rho因数分解

namespace Fctr{
	
	#define int long long 
		
	set<int > Prime;
	
	mt19937 rnd(random_device{}());
	
	int qmul(int x, int y, int p){
		int r = x * y - p * (int)(1.L / p * x * y);
		return r - p * (r >= p) + p * (r < 0);
	}
	
	int quickpow(int x, int y, int p){
		int r = 1;
		for(; y; y >>= 1, x = qmul(x, x, p))
			if(y & 1) r = qmul(r, x, p);
		return r;
	}
	
	bool MR(int n){
		
		if(n < 4) return n > 1;
		
		int k = __builtin_ctz(n - 1), m = n - 1 >> k;
		for(int i = 0; i <= 4; i ++){
			int e = quickpow(rnd() % (n - 3) + 2, m, n), j = k;
			for(; e != 1 && e != n - 1 && j --; e = qmul(e, e, n));
			if(e != n - 1 && j != k)
				return false;
		}
		return true;
	}
	
	int f(int x, int n){
		return qmul(x, x, n) + 1;
	}
	
	int PR(int n){
		
	  int p = 2, q;
	  for(int x = 0, y = 0, i = 0; i ++ % 255 || __gcd(p, n) == 1;){ // !!!
	  	 x = f(x, n), y = f(f(y, n), n);
	    if(x == y){
	   		x = rnd() % n;
			y = f(x, n);	
		}
		
		q = qmul(p,abs(x-y),n);
	    if(q)p=q;
		}
		return __gcd(p,n);
	}
	
	void fctr(int n){
		if(MR(n)){
			Prime.insert(n);
			return ;
		}
		
		int e = PR(n);
		fctr(n / e);
		fctr(e);
	}
}

*2022/10/5 新增积性函数在组合数处的取值
*2022/10/16 修改Miller Rabin素数测试 & Pollard Rho因数分解模板，速度大幅度优化
*2022/11/24 新增了莫比乌斯反演和欧拉反演的诸多公式
*2022/11/24 改进了min_25筛模板