math.go

package copypasta

import (
	. "fmt"
	"math"
	"math/big"
	"math/bits"
	"math/rand"
	"sort"
)

/* 数论 组合数学 博弈论 积分 插值 随机算法

一些不等式及其证明 https://www.luogu.com.cn/blog/chinesepikaync/oi-zhong-kuai-yong-dao-di-yi-suo-fou-deng-shi-ji-ji-zheng-ming

https://en.wikipedia.org/wiki/List_of_recreational_number_theory_topics

https://euler.stephan-brumme.com/toolbox/

todo 期望与贡献 https://codeforces.com/blog/entry/62690 https://codeforces.com/blog/entry/62792
     一类概率期望问题的杀器：势函数和鞅的停时定理 https://www.cnblogs.com/TinyWong/p/12887591.html https://codeforces.com/blog/entry/87598 最后一题

NOTE: a%-b == a%b
NOTE: 对于整数来说有
       ax≤by  =>  x≤⌊by/a⌋
       ax<by  =>  x<⌈by/a⌉
       ax≥by  =>  x≥⌈by/a⌉
       ax>by  =>  x>⌊by/a⌋
NOTE: ⌊⌊x/n⌋/m⌋ = ⌊x/(n*m)⌋
NOTE: ⌈⌈x/n⌉/m⌉ = ⌈x/(n*m)⌉

AP: Sn = n*(2*a1+(n-1)*d)/2
GP: Sn = a1*(q^n-1)/(q-1), q!=1
       = a1*n, q==1
∑i*q^(i-1) = n*q^n - (q^n-1)/(q-1)

https://oeis.org/A001787 n*2^(n-1) = ∑i*C(n,i)   number of ones in binary numbers 1 to 111...1 (n bits)
https://oeis.org/A000337 ∑i*2^(i-1) = (n-1)*2^n+1
https://oeis.org/A036799 ∑i*2^i = (n-1)*2^(n+1)+2 = A000337(n)*2

https://oeis.org/A027992 a(n) = 2^n*(3n-1)+2 = The total sum of squares of parts in all compositions of n （做 https://codeforces.com/problemset/problem/235/B 时找到的序列）
https://oeis.org/A271638 a(n) = (13*n-36)*2^(n-1)+6*n+18 = 	The total sum of the cubes of all parts of all compositions of n

https://en.wikipedia.org/wiki/Faulhaber%27s_formula
https://oeis.org/A000330 平方和 = n*(n+1)*(2*n+1)/6
https://oeis.org/A000537 立方和 = (n*(n+1)/2)^2

https://oeis.org/A061168 Σfloor(log2(i)) = Σ(bits.Len(i)-1)

∑∑|ai-aj|
= 2*∑(i*ai-preSum(i-1)), i=[0,n-1], a 需要排序
https://www.luogu.com.cn/blog/DPair2005/solution-cf340c
https://codeforces.com/problemset/problem/340/C

https://oeis.org/A333885 Number of triples (i,j,k) with 1 <= i < j < k <= n such that i divides j divides k https://ac.nowcoder.com/acm/contest/7613/A

https://oeis.org/A000295 Eulerian numbers: Sum_{k=0..n} (n-k)*2^k = 2^n - n - 1
	Number of permutations of {1,2,...,n} with exactly one descent
	Number of partitions of an n-set having exactly one block of size > 1
	a(n-1) is the number of subsets of {1..n} in which the largest element of the set exceeds by at least 2 the next largest element
		For example, for n = 5, a(4) = 11 and the 11 sets are {1,3}, {1,4}, {1,5}, {2,4}, {2,5}, {3,5}, {1,2,4}, {1,2,5}, {1,3,5}, {2,3,5}, {1,2,3,5}
	a(n-1) is also the number of subsets of {1..n} in which the second smallest element of the set exceeds by at least 2 the smallest element
		For example, for n = 5, a(4) = 11 and the 11 sets are {1,3}, {1,4}, {1,5}, {2,4}, {2,5}, {3,5}, {1,3,4}, {1,3,5}, {1,4,5}, {2,4,5}, {1,3,4,5}

https://oeis.org/A064413 EKG sequence (or ECG sequence)
a(1) = 1; a(2) = 2; for n > 2, a(n) = smallest number not already used which shares a factor with a(n-1)

https://oeis.org/A002326 least m > 0 such that 2n+1 divides 2^m-1
https://leetcode-cn.com/problems/minimum-number-of-operations-to-reinitialize-a-permutation/

http://oeis.org/A003136 Loeschian number: numbers of the form x^2 + xy + y^2
https://en.wikipedia.org/wiki/Loeschian_number
https://www.bilibili.com/video/BV1or4y1A76q

数的韧性 https://en.wikipedia.org/wiki/Persistence_of_a_number 乘法: https://oeis.org/A003001 加法: https://oeis.org/A006050

Smallest number h such that n*h is a repunit (111...1), or 0 if no such h exists
https://oeis.org/A190301 111...1
https://oeis.org/A216485 222...2
相关题目 https://atcoder.jp/contests/abc174/tasks/abc174_c  快速算法见 https://img.atcoder.jp/abc174/editorial.pdf

	Least k such that the decimal representation of k*n contains only 1's and 0's
	https://oeis.org/A079339
	0's and d's (2~9): A096681-A096688

	a(n) is the least value of k such that k*n uses only digits 1 and 2. a(n) = -1 if no such multiple exists
	https://oeis.org/A216482

	a(n) is the smallest positive number such that the decimal digits of n*a(n) are all 0, 1 or 2
	https://oeis.org/A181061

椭圆曲线加密算法 https://ac.nowcoder.com/acm/contest/6916/C

挑战 2.6 节练习题
2429 分解 LCM/GCD = a*b 且 gcd(a,b)=1 且 a+b 最小
1930 https://www.luogu.com.cn/problem/UVA10555 https://www.luogu.com.cn/problem/SP1166 floatToRat
3126 https://www.luogu.com.cn/problem/UVA12101 https://www.luogu.com.cn/problem/SP1841 BFS
3421 质因数幂次和 可重排列
3292 https://www.luogu.com.cn/problem/UVA11105 在 Z={4k+1} 上筛素数
3641 https://www.luogu.com.cn/problem/UVA11287 Carmichael Numbers https://oeis.org/A002997
4.1 节练习题（模运算的世界）
1150 https://www.luogu.com.cn/problem/UVA10212
1284
2115
3708
2720
GCJ Japan 2011 Final B

CF tag https://codeforces.com/problemset?order=BY_RATING_ASC&tags=number+theory
CF tag https://codeforces.com/problemset?order=BY_RATING_ASC&tags=combinatorics

*/

func numberTheoryCollection() {
	const mod int64 = 1e9 + 7 // 998244353
	pow := func(x, n, p int64) (res int64) {
		x %= p
		res = 1
		for ; n > 0; n >>= 1 {
			if n&1 == 1 {
				res = res * x % p
			}
			x = x * x % p
		}
		return
	}

	/* GCD LCM 相关

	TIPS: 一般 LCM 的题目都需要用 LCM=x*y/GCD 转换成研究 GCD 的性质

	GCD 与质因子 https://codeforces.com/problemset/problem/264/B
	数组中最小的 LCM(ai,aj) https://codeforces.com/problemset/problem/1154/G
	分拆与 LCM  https://ac.nowcoder.com/acm/contest/5961/D https://ac.nowcoder.com/discuss/439005

	todo https://codeforces.com/contest/1462/problem/D 的 O(nlogn) 解法

	*/

	gcd := func(a, b int64) int64 {
		for a != 0 {
			a, b = b%a, a
		}
		return b
	}

	// 例题 https://nanti.jisuanke.com/t/A1633
	gcdPrefix := func(a []int64) []int64 {
		n := len(a)
		gp := make([]int64, n+1)
		for i, v := range a {
			gp[i+1] = gcd(gp[i], v)
		}
		return gp
	}
	gcdSuffix := func(a []int64) []int64 {
		n := len(a)
		gs := make([]int64, n+1)
		for i := n - 1; i >= 0; i-- {
			gs[i] = gcd(gs[i+1], a[i])
		}
		return gs
	}

	lcm := func(a, b int64) int64 { return a / gcd(a, b) * b }

	// 前 n 个数的 LCM https://oeis.org/A003418 a(n) = lcm(1,...,n) ~ exp(n)
	// 1, 2, 6, 12, 60, 60, 420, 840, 2520, 2520, 27720, 27720, 360360, 360360, 360360, 720720, 12252240, 12252240, 232792560, 232792560, 232792560, 232792560, 5354228880, 5354228880, 26771144400, 26771144400, 80313433200, 80313433200
	//     相关题目 https://atcoder.jp/contests/arc110/tasks/arc110_a
	//             https://codeforces.com/contest/1485/problem/D
	// a(n)/a(n-1) = https://oeis.org/A014963
	//     前缀和 https://oeis.org/A072107 https://ac.nowcoder.com/acm/contest/7607/A
	// LCM(2, 4, 6, ..., 2n) https://oeis.org/A051426
	// Mangoldt Function https://mathworld.wolfram.com/MangoldtFunction.html
	// a(n) 的因子个数 d(lcm(1,...,n)) https://oeis.org/A056793
	//     这同时也是 1~n 的子集的 LCM 的种类数

	// GCD 性质统计相关
	// NOTE: 对于一任意非负序列，前 i 个数的 GCD 是非增序列，且至多有 O(logMax) 个不同值
	//       应用：https://codeforces.com/problemset/problem/1210/C
	// #{(a,b) | 1<=a<=b<=n, gcd(a,b)=1}   https://oeis.org/A002088
	//     = ∑phi(i)
	// #{(a,b) | 1<=a,b<=n, gcd(a,b)=1}   https://oeis.org/A018805
	//     = 2*(∑phi(i))-1
	//     = 2*A002088(n)-1
	// #{(a,b,c) | 1<=a,b,c<=n, gcd(a,b,c)=1}   https://oeis.org/A071778
	//     = ∑mu(i)*floor(n/i)^3
	// #{(a,b,c,d) | 1<=a,b,c,d<=n, gcd(a,b,c,d)=1}   https://oeis.org/A082540
	//     = ∑mu(i)*floor(n/i)^4

	// GCD 求和相关
	// ∑gcd(n,i) = ∑{d|n}d*phi(n/d)          https://oeis.org/A018804 https://www.luogu.com.cn/problem/P2303
	//     更简化的公式见小粉兔博客 https://www.cnblogs.com/PinkRabbit/p/8278728.html
	// ∑n/gcd(n,i) = ∑{d|n}d*phi(d)          https://oeis.org/A057660
	// ∑∑gcd(i,j) = ∑phi(i)*(floor(n/i))^2   https://oeis.org/A018806   https://www.luogu.com.cn/problem/P2398
	// ∑∑∑gcd(i,j,k) = ∑phi(i)*(floor(n/i))^3   https://ac.nowcoder.com/acm/contest/7608/B

	// LCM 求和相关
	// ∑lcm(n,i) = n*(1+∑{d|n}d*phi(d))/2 = n*(1+A057660(n))/2   https://oeis.org/A051193
	// ∑lcm(n,i)/n = A051193(n)/n = (1+∑{d|n}d*phi(d))/2 = (1+A057660(n))/2   https://oeis.org/A057661
	// ∑∑lcm(i,j)   https://oeis.org/A064951

	// 统计数组的所有子序列的 GCD 的不同个数，复杂度 O(Clog^2C)
	// LC1819/周赛235D https://leetcode-cn.com/problems/number-of-different-subsequences-gcds/
	countDifferentSubsequenceGCDs := func(a []int, gcd func(int, int) int) (ans int) {
		const mx int = 4e5 //
		has := [mx + 1]bool{}
		for _, v := range a {
			has[v] = true
		}
		for i := 1; i <= mx; i++ {
			g := 0
			for j := i; j <= mx; j += i {
				if has[j] {
					g = gcd(g, j)
				}
			}
			if g == i {
				ans++
			}
		}
		return
	}

	// 最简分数
	type pair struct{ x, y int64 }
	frac := func(a, b int64) pair { g := gcd(a, b); return pair{a / g, b / g} }

	// 类欧几里得算法
	// ∑⌊(ai+b)/m⌋, i in [0,n-1]
	// todo https://www.luogu.com.cn/blog/AlanWalkerWilson/Akin-Euclidean-algorithm-Basis
	// todo https://www.luogu.com.cn/blog/Shuchong/qian-tan-lei-ou-ji-li-dei-suan-fa
	//
	// 模板题 https://atcoder.jp/contests/practice2/tasks/practice2_c
	//       https://www.luogu.com.cn/problem/P5170
	// todo https://codeforces.com/problemset/problem/1182/F
	//  https://codeforces.com/problemset/problem/1098/E
	floorSum := func(n, m, a, b int64) (res int64) {
		if a < 0 {
			a2 := a%m + m
			res -= n * (n - 1) / 2 * ((a2 - a) / m)
			a = a2
		}
		if b < 0 {
			b2 := b%m + m
			res -= n * ((b2 - b) / m)
			b = b2
		}
		for {
			if a >= m {
				res += n * (n - 1) / 2 * (a / m)
				a %= m
			}
			if b >= m {
				res += n * (b / m)
				b %= m
			}
			yMax := a*n + b
			if yMax < m {
				break
			}
			n = yMax / m
			b = yMax % m
			m, a = a, m
		}
		return
	}

	sqCheck := func(a int64) bool { r := int64(math.Round(math.Sqrt(float64(a)))); return r*r == a }
	cubeCheck := func(a int64) bool { r := int64(math.Round(math.Cbrt(float64(a)))); return r*r*r == a }
	// 平方数开平方
	sqrt := func(a int64) int64 {
		r := int64(math.Round(math.Sqrt(float64(a))))
		if r*r == a {
			return r
		}
		return -1
	}
	// 立方数开立方
	cbrt := func(a int64) int64 {
		r := int64(math.Round(math.Cbrt(float64(a))))
		if r*r*r == a {
			return r
		}
		return -1
	}

	// 返回差分表的最后一个数
	// return the bottom entry in the difference table
	// 另一种做法是用公式 ∑(-1)^i * C(n,i) * a_i, i=0..n-1
	bottomDiff := func(a []int) int {
		for ; len(a) > 1; a = a[:len(a)-1] {
			for i := 0; i+1 < len(a); i++ {
				a[i] = a[i+1] - a[i]
			}
		}
		return a[0]
	}

	/* 质数 质因子分解 */

	// n/2^k http://oeis.org/A000265

	// 质数表 https://oeis.org/A000040
	// primes[i]%10 http://oeis.org/A007652
	// 10-primes[i]%10 http://oeis.org/A072003
	// p-1 http://oeis.org/A006093
	// p+1 http://oeis.org/A008864
	// p^2+p+1 http://oeis.org/A060800 = sigma(p^2)
	primes := []int{
		2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
		101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199,
		211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293,
		307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
		401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
		503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599,
		601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691,
		701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
		809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887,
		907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, /* #=168 */
		1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
		1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193,
		1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297,
		1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399,
		1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
		1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597,
		1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699,
		1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
		1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
		1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, /* #=303 */
	}

	// map{小于 10^n 的素数个数: 小于 10^n 的最大素数} http://oeis.org/A006880 http://oeis.org/A003618   10^n-a(n): http://oeis.org/A033874
	primes10 := map[int]int64{
		4:         7,
		25:        97,
		168:       997, // 1e3
		1229:      9973,
		9592:      99991,
		78498:     999983, // 1e6
		664579:    9999991,
		5761455:   99999989,
		50847534:  999999937, // 1e9
		455052511: 9999999967,
	}

	// 大于 10^n 的最小素数 http://oeis.org/A090226 http://oeis.org/A003617    a(n)-10^n: http://oeis.org/A033873
	primes10_ := []int64{
		2,
		11,
		101,
		1009, // 1e3
		10007,
		100003,
		1000003, // 1e6
		10000019,
		100000007,
		1000000007, //1e9
		10000000019,
	}

	/* 质数性质统计相关

	质数的幂次组成的集合 {p^k} https://oeis.org/A000961
		补集 https://oeis.org/A024619
		Exponential of Mangoldt function https://oeis.org/A014963

	质数前缀和 http://oeis.org/A007504
	a(n) ~ n^2 * log(n) / 2
	a(n)^2 - a(n-1)^2 = A034960(n)
	EXTRA: divide odd numbers into groups with prime(n) elements and add together http://oeis.org/A034960
		仍然是质数的前缀和 http://oeis.org/A013918 对应的前缀和下标 http://oeis.org/A013916

	质数前缀积 prime(n)# https://oeis.org/A002110
	the least number with n distinct prime factors
	2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870, /9/
	6469693230, 200560490130, 7420738134810, 304250263527210, 13082761331670030, 614889782588491410

	质数差分 Gap https://oeis.org/A001223
		Positions of records https://oeis.org/A005669
		Values of records https://oeis.org/A005250
		Gap 均值 https://oeis.org/A286888 a(n)= floor((prime(n) - 2)/(n - 1))
		相关题目 https://www.luogu.com.cn/problem/P6104 https://class.luogu.com.cn/classroom/lgr69
		Numbers whose distance to the closest prime number is a prime number http://oeis.org/A160666

	任意质数之差 https://oeis.org/A030173
	非任意质数之差 https://oeis.org/A007921

	质数的逆二项变换 Inverse binomial transform of primes https://oeis.org/A007442

	合数前缀和 https://oeis.org/A053767

	合数前缀积 Compositorial number https://oeis.org/A036691

	不与质数相邻的合数 http://oeis.org/A079364

	半素数 https://oeis.org/A001358 也叫双素数/二次殆素数 Semiprimes (or biprimes): products of two primes
		https://en.wikipedia.org/wiki/Semiprime
		https://en.wikipedia.org/wiki/Almost_prime
		非平方半素数 https://oeis.org/A006881 Squarefree semiprimes: Numbers that are the product of two distinct primes.

	绝对素数 https://oeis.org/A003459 各位数字可以任意交换位置，其结果仍为素数
		https://en.wikipedia.org/wiki/Permutable_prime

	哥德巴赫猜想 - 偶数分拆的最小质数 Goldbach’s conjecture https://oeis.org/A020481
	由质数分布可知选到一对质数的概率是 O(1/ln^2(n))
	https://en.wikipedia.org/wiki/Goldbach%27s_conjecture
		Positions of records https://oeis.org/A025018
		Values of records https://oeis.org/A025019
		1e9 内最大的为 a(721013438) = 1789
		2e9 内最大的为 a(1847133842) = 1861

	勒让德猜想 - 在两个相邻平方数之间，至少有一个质数 Legendre’s conjecture
	https://en.wikipedia.org/wiki/Legendre%27s_conjecture
	Number of primes between n^2 and (n+1)^2 https://oeis.org/A014085
	Number of primes between n^3 and (n+1)^3 https://oeis.org/A060199

	伯特兰-切比雪夫定理 - n ~ 2n 之间至少有一个质数 Bertrand's postulate
	https://en.wikipedia.org/wiki/Bertrand%27s_postulate
	Number of primes between n and 2n (inclusive) https://oeis.org/A035250
	Number of primes between n and 2n exclusive https://oeis.org/A060715

	Least k such that H(k) > n, where H(k) is the harmonic number Σ{i=1..k} 1/i
	https://oeis.org/A002387
	https://oeis.org/A004080

		a(n) = smallest prime p such that Σ{primes q = 2, ..., p} 1/q exceeds n
		5, 277, 5_195_977, 1801241230056600523
		https://oeis.org/A016088 pi
		https://oeis.org/A046024 i

	a(n) = largest m such that the harmonic number H(m)= Σ{i=1..m} 1/i is < n
	https://oeis.org/A115515

		a(n) = largest prime p such that Σ{primes q = 2, ..., p} 1/q does not exceed n
		3, 271, 5_195_969, 1801241230056600467
		https://oeis.org/A223037

	Exponent of highest power of 2 dividing n, a.k.a. the binary carry sequence, the ruler sequence, or the 2-adic valuation of n
	a(n) = 0 if n is odd, otherwise 1 + a(n/2)
	http://oeis.org/A007814

	https://oeis.org/A000043 Mersenne exponents: primes p such that 2^p - 1 is prime. Then 2^p - 1 is called a Mersenne prime

	*/

	// 判断一个数是否为质数
	isPrime := func(n int64) bool {
		for i := int64(2); i*i <= n; i++ {
			if n%i == 0 {
				return false
			}
		}
		return n >= 2
	}
	// https://www.luogu.com.cn/problem/U82118
	isPrime = func(n int64) bool { return big.NewInt(n).ProbablyPrime(0) }

	// 判断质数+求最大质因子
	// 先用 Pollard-Rho 算法求出一个因子，然后递归求最大质因子
	// https://zhuanlan.zhihu.com/p/267884783
	// https://www.luogu.com.cn/problem/P4718
	pollardRho := func(n int64) int64 {
		if n == 4 {
			return 2
		}
		if isPrime(n) {
			return n
		}
		abs := func(x int64) int64 {
			if x < 0 {
				return -x
			}
			return x
		}
		mul := func(a, b int64) (res int64) {
			for ; b > 0; b >>= 1 {
				if b&1 == 1 {
					res = (res + a) % n
				}
				a = (a + a) % n
			}
			return
		}
		for {
			c := 1 + rand.Int63n(n-1)
			f := func(x int64) int64 { return (mul(x, x) + c) % n }
			for t, r := f(0), f(f(0)); t != r; t, r = f(t), f(f(r)) {
				if d := gcd(abs(t-r), n); d > 1 {
					return d
				}
			}
		}
	}
	{
		max := func(a, b int64) int64 {
			if a > b {
				return a
			}
			return b
		}
		cacheGPF := map[int64]int64{}
		var gpf func(int64) int64
		gpf = func(x int64) (res int64) {
			if cacheGPF[x] > 0 {
				return cacheGPF[x]
			}
			defer func() { cacheGPF[x] = res }()
			p := pollardRho(x)
			if p == x {
				return p
			}
			return max(gpf(p), gpf(x/p))
		}
	}

	// 预处理: [2,mx] 范围内的质数
	// 埃拉托斯特尼筛法 Sieve of Eratosthenes
	// 质数个数 π(n) https://oeis.org/A000720
	//         π(10^n) https://oeis.org/A006880
	//         4, 25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534, /* 1e9 */
	//         455052511, 4118054813, 37607912018, 346065536839, 3204941750802, 29844570422669, 279238341033925, 2623557157654233, 24739954287740860, 234057667276344607,
	sieve := func() {
		const mx int = 1e6
		primes := []int{}
		pid := [mx + 1]int{-1, -1}
		for i := 2; i <= mx; i++ {
			if pid[i] == 0 {
				primes = append(primes, i)
				pid[i] = len(primes)
				for j := 2 * i; j <= mx; j += i {
					pid[j] = -1
				}
			}
		}

		// EXTRA: pi(n), the number of primes <= n https://oeis.org/A000720
		pi := [mx + 1]int{}
		for i := 2; i <= mx; i++ {
			pi[i] = pi[i-1]
			if pid[i] > 0 {
				pi[i]++
			}
		}
	}

	// 线性筛 欧拉筛
	// 每个合数都从其 LPF 访问到
	// 参考 https://oi-wiki.org/math/sieve/ 以及进阶指南 p.136-137
	// https://www.luogu.com.cn/problem/P3383
	sieveEuler := func() {
		const mx int = 1e7
		primes := []int{}
		pid := [mx + 1]int{-1, -1}
		for i := 2; i <= mx; i++ {
			if pid[i] == 0 {
				primes = append(primes, i)
				pid[i] = len(primes)
			}
			for _, p := range primes {
				if p*i > mx {
					break
				}
				pid[p*i] = -1
				if i%p == 0 {
					break
				}
			}
		}
	}

	// 一般线性筛的模板
	// 记 f(n) 为积性函数
	// 其满足 1. f(p) = p ...
	//       2. f(p^(k+1)) = f(p^k) ... p
	//       3. f(x*p) = f(x) ... p (p 不是 x 的因子)
	// 一个典型的例子见下面 σ(n) 的线性筛求法
	sieveEulerTemplate := func() []int {
		const mx int = 1e7
		f := make([]int, mx+1)
		f[1] = 1 //
		vis := make([]bool, mx+1)
		primes := []int{}
		for i := 2; i <= mx; i++ {
			if !vis[i] {
				// 1: p
				f[i] = i
				primes = append(primes, i)
			}
			for _, p := range primes {
				v := p * i
				if v > mx {
					break
				}
				vis[v] = true
				if i%p == 0 {
					// 2: p^(k+1) <- p^k
					f[v] = f[i] * p
					break
				}
				// 3: x*p <- x 且 x 的质因子是没有 p 的
				f[v] = f[i] * p
			}
		}
		return f
	}

	// 区间筛法
	// 预处理 [2,√R] 的所有质数，去筛 [L,R] 之间的质数

	// 质因数分解 prime factorization
	// 返回分解出的质数及其指数
	// https://mathworld.wolfram.com/PrimeFactorization.html
	// todo 更高效的算法 - Pollard's Rho
	// n 的质因数分解中 2 的幂次 http://oeis.org/A007814
	// n 的质因数分解中非 2 的幂次之和 http://oeis.org/A087436
	type factor struct {
		p int64
		e int
	}
	factorize := func(x int64) (factors []factor) {
		for i := int64(2); i*i <= x; i++ {
			e := 0
			for ; x%i == 0; x /= i {
				e++
			}
			if e > 0 {
				factors = append(factors, factor{i, e})
			}
		}
		if x > 1 {
			factors = append(factors, factor{x, 1})
		}
		return
	}
	primeDivisors := func(x int) (primes []int) {
		for i := 2; i*i <= x; i++ {
			if x%i == 0 {
				for x /= i; x%i == 0; x /= i {
				}
				primes = append(primes, i)
			}
		}
		if x > 1 {
			primes = append(primes, x)
		}
		return
	}

	// 阶乘的质因数分解 Finding Power of Factorial Divisor
	// 见进阶指南 p.138
	// https://cp-algorithms.com/algebra/factorial-divisors.html
	// https://codeforces.com/contest/1114/problem/C
	powerOfFactorialPrimeDivisor := func(n, p int64) (k int64) {
		for pp := p; ; pp *= p {
			k += n / pp
			if pp > n/p { // 注意判断方式，防止爆 int64
				break
			}
		}
		return
	}

	// 预处理: [2,mx] 的质因数分解的系数和 bigomega(n) or Omega(n) https://oeis.org/A001222
	// https://en.wikipedia.org/wiki/Prime_omega_function
	// a(n) depends only on prime signature of n (cf. https://oeis.org/A025487)
	// So a(24) = a(375) since 24 = 2^3 * 3 and 375 = 3 * 5^3 both have prime signature (3, 1)
	//
	// 		Omega(n) - omega(n) https://oeis.org/A046660
	//
	// 另一种写法 https://math.stackexchange.com/questions/1955105/corectness-of-prime-factorization-over-a-range
	// 性质：Omega(nm)=Omega(n)+Omega(m)
	// 前缀和 https://oeis.org/A022559 = Omega(n!) ~ O(nloglogn)
	// EXTRA: https://oeis.org/A005361 Product of exponents of prime factorization of n
	//        https://oeis.org/A135291 Product of exponents of prime factorization of n!
	primeExponentsCountAll := func() {
		const mx int = 1e6
		Omega := [mx + 1]int{} // int8
		primes := []int{}
		for i := 2; i <= mx; i++ {
			if Omega[i] == 0 {
				Omega[i] = 1
				primes = append(primes, i)
			}
			for _, p := range primes {
				if p*i > mx {
					break
				}
				Omega[p*i] = Omega[i] + 1
			}
		}

		// EXTRA: 前缀和，即 Omega(n!) https://oeis.org/A022559
		for i := 3; i <= mx; i++ {
			Omega[i] += Omega[i-1]
		}
	}

	/* 因子/因数/约数

	n 的因子个数 d(n) = Π(ei+1), ei 为第 i 个质数的系数 https://oeis.org/A000005 d(n) 也写作 τ(n)
		Positions of records (高合成数，反素数) https://oeis.org/A002182
		Values of records https://oeis.org/A002183
		相关题目：范围内的最多约数个数 https://www.luogu.com.cn/problem/P1221

		max(d(i)), i=1..10^n https://oeis.org/A066150
			方便估计复杂度 - 近似为开立方
			4, 12, 32, 64, 128, /5/
	        240, 448, 768, 1344, /9/
			2304, 4032, 6720, 10752, 17280, 26880, 41472, 64512, 103680, 161280, /19/

			上面这些数对应的最小的 n https://oeis.org/A066151
			6, 60, 840, 7560, 83160, 720720, 8648640, 73513440, 735134400,
			6983776800, 97772875200, 963761198400, 9316358251200, 97821761637600, 866421317361600, 8086598962041600, 74801040398884800, 897612484786617600

		d(n) 前缀和 = Σ{k=1..n} floor(n/k) https://oeis.org/A006218
	               = 见后文「数论分块/除法分块」

		n+d(n) https://oeis.org/A062249
		n-d(n) https://oeis.org/A049820   count https://oeis.org/A060990   前缀和 https://oeis.org/A161664
		n*d(n) https://oeis.org/A038040   前缀和 https://oeis.org/A143127 = Sum_{i=1..floor(√n)}i*(i+floor(n/i))*(floor(n/i)+1-i) - 平方和(floor(√n))
												https://atcoder.jp/contests/abc172/tasks/abc172_d
		n*n*d(n) https://oeis.org/A034714   前缀和 https://oeis.org/A319085
		d(n)|n https://oeis.org/A033950 refactorable numbers / tau numbers
			n/d(n) https://oeis.org/A036762
		n%d(n) https://oeis.org/A054008
		a(1)=1, a(n+1)=a(n)+d(a(n)) https://oeis.org/A064491
		Smallest k such that d(k) = n https://oeis.org/A005179
			a(p) = 2^(p-1) for primes p
			相关题目 https://codeforces.com/problemset/problem/27/E https://www.luogu.com.cn/problem/P1128
			质数的情况 https://oeis.org/A061286
	    Number of divisors of n^2 less than n https://oeis.org/A063647 Also number of ways to write 1/n as a difference of exactly 2 unit fractions
	        a(n) = (d(n^2)-1)/2

	n 的因子之和 σ(n) = Π(pi^(ei+1)-1)/(pi-1) https://oeis.org/A000203 todo 相关题目 https://www.luogu.com.cn/problem/P1593
	    线性筛求法见后面
	    max(σ(i)), i=1..10^n https://oeis.org/A066410
	         对应的 n https://oeis.org/A066424
	    Smallest k such that sigma(k) = n http://oeis.org/A051444
		σ(n) 前缀和 = Σ{k=1..n} k*floor(n/k) https://oeis.org/A024916
		真因子之和 https://oeis.org/A001065 σ(n)-n
		完全数/完美数/完备数 https://oeis.org/A000396 Perfect numbers (σ(n) = 2n)
			https://en.wikipedia.org/wiki/Perfect_number
		过剩数/丰数/盈数 https://oeis.org/A005101 Abundant numbers (σ(n) > 2n)
			https://en.wikipedia.org/wiki/Abundant_number
		亏数/缺数/不足数 https://oeis.org/A005100 Deficient numbers (σ(n) < 2n)
			https://en.wikipedia.org/wiki/Deficient_number
			https://ac.nowcoder.com/acm/contest/10322/A O(nlogn) 可以先预处理因子

	n 的因子之积 μ(n) = n^(d(n)/2.0) https://oeis.org/A007955
	because we can form d(n)/2 pairs from the factors, each with product n
		若 n 是完全平方数，则 ei+1 全为奇数，此时可以计算 [n^(1/2)]^d(n)
		否则，ei+1 中必有偶数，将其除二
		相关题目 https://codeforces.com/problemset/problem/615/D

	n 的因子的差分表的最后一个数 https://oeis.org/A187202 https://oeis.org/A187203
	NOTE: a(2^k) = 1

		正数 https://oeis.org/A193671
		零   https://oeis.org/A187204
		负数 https://oeis.org/A193672

	d(d(...d(n))) 迭代至 2 所需要的迭代次数
	0,0,1,0,2,0,2,1,2,0,3,0,2,2,1,0,3,0,3,2,2,0,3,1,2,2,3

	高合成数/反质数 Highly Composite Numbers https://oeis.org/A002182
	https://oi-wiki.org/math/prime/#_7
	性质：一个高合成数一定是由另一个高合成数乘一个质数得到
	见进阶指南 p.140-141
	Number of divisors of n-th highly composite number https://oeis.org/A002183
	Number of highly composite numbers not divisible by n https://oeis.org/A199337
	求出不超过 n 的最大的反质数 https://www.luogu.com.cn/problem/P1463

	Largest divisor of n having the form 2^i*5^j http://oeis.org/A132741
	a(n) = A006519(n)*A060904(n) = 2^A007814(n)*5^A112765(n)

	Squarefree numbers https://oeis.org/A005117 (介绍了一种筛法)
	Numbers that are not divisible by a square greater than 1
	Lim_{n->infinity} a(n)/n = Pi^2/6

		Numbers that are not squarefree https://oeis.org/A013929
		Numbers that are divisible by a square greater than 1

	a(n) = Min {m>n | m has same prime factors as n ignoring multiplicity} https://oeis.org/A065642
		Numbers such that a(n)/n is not an integer are listed in https://oeis.org/A284342
	*/

	// 枚举一个数的全部因子
	divisors := func(n int64) (ds []int64) {
		for d := int64(1); d*d <= n; d++ {
			if n%d == 0 {
				ds = append(ds, d)
				if d*d < n {
					ds = append(ds, n/d)
				}
			}
		}
		//sort.Slice(ds, func(i, j int) bool { return ds[i] < ds[j] })
		return
	}

	// 不需要排序的写法
	divisors = func(n int64) (ds []int64) {
		ds2 := []int64{}
		for d := int64(1); d*d <= n; d++ {
			if n%d == 0 {
				ds = append(ds, d)
				if d*d < n {
					ds2 = append(ds2, n/d)
				}
			}
		}
		for i := len(ds2) - 1; i >= 0; i-- {
			ds = append(ds, ds2[i])
		}
		return
	}

	divisorPairs := func(n int64) (ds [][2]int64) {
		for d := int64(1); d*d <= n; d++ {
			if n%d == 0 {
				ds = append(ds, [2]int64{d, n / d})
			}
		}
		return
	}

	doDivisors := func(n int, do func(d int)) {
		for d := 1; d*d <= n; d++ {
			if n%d == 0 {
				do(d)
				if d*d < n {
					do(n / d)
				}
			}
		}
	}
	doDivisors2 := func(n int, do func(d1, d2 int)) {
		for d := 1; d*d <= n; d++ {
			if n%d == 0 {
				do(d, n/d)
			}
		}
	}

	// Number of odd divisors of n https://oeis.org/A001227
	// a(n) = d(2*n) - d(n)
	// 亦为整数 n 分拆成若干连续整数的方法数
	// Number of partitions of n into consecutive positive integers including the trivial partition of length 1
	// e.g. 9 = 2+3+4 or 4+5 or 9 so a(9)=3
	oddDivisorsNum := func(n int) (ans int) {
		for i := 1; i*i <= n; i++ {
			if n%i == 0 {
				if i&1 == 1 {
					ans++
				}
				if i*i < n && n/i&1 == 1 {
					ans++
				}
			}
		}
		return
	}

	// 因子的中位数（偶数个因子时取小的那个）
	// Lower central (median) divisor of n https://oeis.org/A060775
	// EXTRA: Largest divisor of n <= sqrt(n) https://oeis.org/A033676
	maxSqrtDivisor := func(n int) int {
		for d := int(math.Sqrt(float64(n))); ; d-- {
			if n%d == 0 {
				return d
			}
		}
	}

	// 预处理: [1,mx] 范围内数的所有因子
	// 复杂度 O(nlogn)
	// NOTE: 1~n 的因子个数总和大约为 nlogn
	// NOTE: divisors[x] 为奇数 => x 为完全平方数 https://oeis.org/A000290
	// NOTE: halfDivisors(x) 为 ≤√x 的因数集合 https://oeis.org/A161906
	divisorsAll := func() {
		const mx int = 1e6
		divisors := [mx + 1][]int{}
		for i := 1; i <= mx; i++ {
			for j := i; j <= mx; j += i {
				divisors[j] = append(divisors[j], i)
			}
		}

		{
			// 统计因子个数 d(n)
			// NOTE: 复杂度可以做到线性 https://codeforces.com/contest/920/submission/76859782
			// 相关 OEIS 见上面的注释块
			const mx int = 1e6
			d := [mx + 1]int{}
			for i := 1; i <= mx; i++ {
				for j := i; j <= mx; j += i {
					d[j]++
				}
			}
		}

		{
			// 去掉 1 作为因子
			const mx = 1e6
			divisors := [mx + 1][]int{1: {1}} // 仅保留 1 的因子 1
			for i := 2; i <= mx; i++ {
				for j := i; j <= mx; j += i {
					divisors[j] = append(divisors[j], i)
				}
			}
		}

		{
			// 线性筛求 n 的因子之和 σ(n)
			// https://codeforces.com/contest/1512/problem/G
			const mx int = 1e7
			d := make([]int, mx+1)
			d[1] = 1
			s := make([]int, mx+1)
			primes := []int{}
			for i := 2; i <= mx; i++ {
				if d[i] == 0 {
					s[i] = 1 + i
					d[i] = s[i]
					primes = append(primes, i)
				}
				for _, p := range primes {
					if p*i > mx {
						break
					}
					if i%p == 0 {
						s[p*i] = s[i]*p + 1
						d[p*i] = d[i] / s[i] * s[p*i]
						break
					}
					s[p*i] = 1 + p
					d[p*i] = d[i] * s[p*i]
				}
			}
		}

		isSquareNumber := func(x int) bool { return len(divisors[x])&1 == 1 }
		halfDivisors := func(x int) []int { d := divisors[x]; return d[:(len(d)-1)/2+1] }

		_, _ = isSquareNumber, halfDivisors
	}

	// LPF(n): least prime dividing n (when n > 1); a(1) = 1 https://oeis.org/A020639
	// 有时候数据范围比较大，用 primeFactorsAll 预处理会 MLE，这时候就要用 LPF 了（同样是预处理但是内存占用低）
	// 先预处理出 LPF，然后对要处理的数 v 不断地除 LPF(v) 直到等于 1
	// 		LPF 前缀和 https://oeis.org/A046669 前缀积 https://oeis.org/A072486
	//		n+LPF(n) https://oeis.org/A061228 the smallest number greater than n which is not coprime to n
	// 		n-LPF(n) https://oeis.org/A046666
	// 			迭代至 0 的次数 https://oeis.org/A175126 相关题目 https://codeforces.com/contest/1076/problem/B
	//		n*LPF(n) https://oeis.org/A285109
	// 		n/LPF(n) https://oeis.org/A032742 即 n 的最大因子 = Max{gcd(n,j); j=n+1..2n-1}
	//
	//		只考虑奇质数 http://oeis.org/A078701
	//
	// GPF(n): greatest prime dividing n, for n >= 2; a(1)=1 https://oeis.org/A006530
	//		GPF(p-1) http://oeis.org/A023503
	//		GPF(p+1) http://oeis.org/A023509
	// 		GPF 前缀和 https://oeis.org/A046670 前缀积 https://oeis.org/A104350
	//		n+GPF(n) https://oeis.org/A070229 the next m>n such that GPF(n)|m
	// 		n-GPF(n) https://oeis.org/A076563
	// 			迭代至 0 的次数 https://oeis.org/A309892
	// 		n*GPF(n) https://oeis.org/A253560
	// 		n/GPF(n) https://oeis.org/A052126
	//      a(1)=1, a(n+1)=a(n)+GPF(a(n)) https://oeis.org/A076271
	//
	// 		n/LPF(n)*GPF(n) https://oeis.org/A130064
	// 		n/GPF(n)*LPF(n) https://oeis.org/A130065
	//
	lpfAll := func() {
		const mx int = 1e6
		lpf := [mx + 1]int{1: 1}
		for i := 2; i <= mx; i++ {
			if lpf[i] == 0 {
				for j := i; j <= mx; j += i {
					if lpf[j] == 0 { // 去掉这个判断就变成求 GPF，也可以用来（从大到小地）分解质因数
						lpf[j] = i
					}
				}
			}
		}

		// EXTRA: 分解 x
		// lpf[x]==p 也可以写成 x%p==0
		factorize := func(x int) {
			for x > 1 {
				p := lpf[x]
				e := 1
				for x /= p; lpf[x] == p; x /= p {
					e++
				}
				// do(p,e) ...

			}
		}

		// EXTRA: https://oeis.org/A007913 Squarefree part of n (also called core(n))
		// a(n) is the smallest positive number m such that n/m is a square
		// https://oeis.org/A013928 Number of (positive) squarefree numbers < n
		// https://oeis.org/A055204 core(n!)
		//     log a(n) ~ n log 2
		//     Square root of largest square dividing n! https://oeis.org/A055772
		// https://oeis.org/A008833 n/core(n)   Largest square dividing n
		// https://oeis.org/A055071 n!/core(n!) Largest square dividing n!
		// https://codeforces.com/contest/1470/problem/B
		// https://codeforces.com/contest/1497/problem/E2
		core := func(x int) int {
			c := 1
			for x > 1 {
				p := lpf[x]
				e := 1
				for x /= p; lpf[x] == p; x /= p {
					e ^= 1
				}
				if e > 0 {
					c *= p
				}
			}
			return c
		}

		// EXTRA: https://oeis.org/A007947 Largest squarefree number dividing n: the squarefree kernel of n, rad(n), radical of n
		// https://oeis.org/A034386 rad(n!) Primorial numbers (second definition): n# = product of primes <= n
		//                                  = rad(LCM(1,...,n))
		//                                  = LCM(core(1), core(2), core(3), ..., core(n))
		// https://oeis.org/A003557 n/rad(n)  n divided by largest squarefree divisor of n
		// https://oeis.org/A049614 n!/rad(n!)
		rad := func(x int) int {
			r := 1
			for x > 1 {
				p := lpf[x]
				r *= p
				for x /= p; lpf[x] == p; x /= p {
				}
			}
			return r
		}

		_, _, _ = factorize, core, rad
	}

	// 预处理: [2,mx] 范围内数的不同质因子，例如 factors[12] = [2,3]
	// for i>=2, factors[i][0] == i means i is prime
	primeFactorsAll := func() {
		const mx int = 1e6
		factors := [mx + 1][]int{}
		for i := 2; i <= mx; i++ {
			if len(factors[i]) == 0 {
				for j := i; j <= mx; j += i {
					factors[j] = append(factors[j], i)
				}
			}
		}
	}

	// 预处理: [2,mx] 的不同的质因子个数 omega(n) https://oeis.org/A001221
	// https://en.wikipedia.org/wiki/Prime_omega_function
	// 莫比乌斯反演 https://oeis.org/A062799 = Sum_{d|n} omega(d)
	distinctPrimesCountAll := func() {
		const mx int = 1e6
		omega := make([]int, mx+1) // int8
		for i := 2; i <= mx; i++ {
			if omega[i] == 0 {
				for j := i; j <= mx; j += i {
					omega[j]++
				}
			}
		}

		{
			// 线性筛
			omega := make([]int, mx+1) // int8
			primes := []int{}
			for i := 2; i <= mx; i++ {
				if omega[i] == 0 {
					omega[i] = 1
					primes = append(primes, i)
				}
				for _, p := range primes {
					if p*i > mx {
						break
					}
					if i%p == 0 {
						omega[p*i] = omega[i]
						break
					}
					omega[p*i] = omega[i] + 1
				}
			}
		}

		// EXTRA: 前缀和，即 omega(n!) https://oeis.org/A013939
		for i := 3; i <= mx; i++ {
			omega[i] += omega[i-1]
		}

		// EXTRA:
		// http://oeis.org/A034444 Number of unitary divisors of n (d such that d divides n, gcd(d, n/d) = 1)
		// a(n) = 1<<omega[n]
		// 另一种说法是，已知 LCM(x,y) 和 GCD(x,y)，求 (x,y) 的数量
		// 由于 (x/GCD) * (y/GCD) = LCM/GCD，且 x/GCD 和 y/GCD 互质
		// 所以相当于是在求 a(LCM/GCD) = 1<<omega[LCM/GCD]
		// 相关题目 https://codeforces.com/problemset/problem/1499/D

		// EXTRA:
		// http://oeis.org/A007875 Number of ways of writing n as p*q, with p <= q, gcd(p, q) = 1
		// a(n) = 1<<(omega[n]-1)
	}

	// 欧拉函数（互质的数的个数）Euler totient function https://oeis.org/A000010
	// https://en.wikipedia.org/wiki/Euler%27s_totient_function
	// https://oi-wiki.org/math/euler/
	// 前缀和见下面的「phi 求和相关」
	// NOTE: φ(φ...(n)) 收敛到 1 的迭代次数是 log 级别的：奇数减一，偶数减半 https://oeis.org/A003434
	// GCD(n, φ(n)) https://oeis.org/A009195
	// n+φ(n) http://oeis.org/A121048
	// n-φ(n) https://oeis.org/A051953 called Cototient
	// n*φ(n) https://oeis.org/A002618 = φ(n^2)
	// φ(n)|n https://oeis.org/A007694 iff n = 1 or n = 2^w * 3^u for w > 0 and u >= 0
	// 		n/φ(n) https://oeis.org/A049237
	//      the quotients can take only 3 distinct values:
	//			n/φ(n) = 1 iff n = 1
	//			n/φ(n) = 2 iff n = 2^w, w >= 1
	//			n/φ(n) = 3 iff n = 2^w * 3^u, w >= 1, u >= 1
	// n%φ(n) https://oeis.org/A068494
	// a(1)=1, a(n+1)=a(n)+φ(a(n)) https://oeis.org/A074693
	// Least number k such that phi(k) = n https://oeis.org/A002181    Inverse of Euler totient function
	// Number of values of k such that phi(k) = n https://oeis.org/A058277
	// φ集合 https://oeis.org/A002202
	// φ补集 https://oeis.org/A007617
	// https://oeis.org/A023900 Dirichlet inverse of Euler totient function
	// https://oeis.org/A070194 小于 n 且与 n 互质的数的最大间隔
	// https://oeis.org/A023896 小于 n 且与 n 互质的数之和 a(n) = phi(n^2)/2 = n*phi(n)/2
	// https://oeis.org/A053818 小于 n 且与 n 互质的数的平方之和 If n = p_1^e_1 * ... *p_r^e_r then a(n) = n^2*phi(n)/3 + (-1)^r*p_1*..._p_r*phi(n)/6 = n^2*A000010(n)/3 + n*A023900(n)/6, n>1
	// https://oeis.org/A053819 小于 n 且与 n 互质的数的立方之和 a(n) = n^2/4*(n*A000010(n) + A023900(n)), n>1

	// phi 求和相关
	// ∑φ(i) https://oeis.org/A002088 #{(x,y): 1<=x<=y<=n, gcd(x,y)=1}
	// a(n) = (3*n^2)/(Pi^2) + O(nlogn)，近似看成 n^2 / 3
	//      = (1/2)*Sum_{k>=1} mu(k)*floor(n/k)*floor(1+n/k)
	//     相关题目 https://codeforces.com/problemset/problem/1009/D
	// 1, 2, 4, 6, 10, 12, 18, 22, 28, 32, 42, 46, 58, 64, 72, 80, 96, 102
	// ∑φ(i)-1 http://oeis.org/A015614 #{(x,y): 1<=x<y<=n, gcd(x,y)=1}
	// 0, 1, 3, 5, 9, 11, 17, 21, 27, 31, 41, 45, 57, 63, 71, 79, 95, 101

	// Number of numbers "unrelated to n" http://oeis.org/A045763
	// m < n such that m is neither a divisor of n nor relatively prime to n
	// a(n) = n + 1 - d(n) - phi(n); where d(n) is the number of divisors of n

	// Unitary totient (or unitary phi) function uphi(n) http://oeis.org/A047994

	// 计算单个数 n 的欧拉函数（互质的数的个数）Euler totient function
	calcPhi := func(n int) int {
		ans := n
		for i := 2; i*i <= n; i++ {
			if n%i == 0 {
				ans = ans / i * (i - 1)
				for n /= i; n%i == 0; n /= i {
				}
			}
		}
		if n > 1 {
			ans = ans / n * (n - 1)
		}
		return ans
	}

	// 预处理 [1,mx] 欧拉函数
	initPhi := func() {
		const mx int = 1e6
		phi := [mx + 1]int{}
		for i := 1; i <= mx; i++ {
			phi[i] = i
		}
		for i := 2; i <= mx; i++ {
			if phi[i] == i {
				for j := i; j <= mx; j += i {
					phi[j] = phi[j] / i * (i - 1)
				}
			}
		}
	}

	// 线性筛 https://oi-wiki.org/math/sieve/#_8
	// https://www.luogu.com.cn/discuss/show/213297
	sievePhi := func() {
		const mx int = 1e6
		phi := [mx + 1]int{1: 1}
		primes := []int{}
		vis := [mx + 1]bool{}
		for i := 2; i <= mx; i++ {
			if !vis[i] {
				phi[i] = i - 1
				primes = append(primes, i)
			}
			for _, p := range primes {
				v := p * i
				if v > mx {
					break
				}
				vis[v] = true
				if i%p == 0 {
					phi[v] = phi[i] * p
					break
				}
				phi[v] = phi[i] * (p - 1)
			}
		}
	}

	// 扩展欧拉定理（降幂公式）
	// https://oi-wiki.org/math/fermat/#_5
	// https://zhuanlan.zhihu.com/p/42632291
	// https://blog.csdn.net/synapse7/article/details/19610361
	// a^b ≡ a^(b mod φ(m)) (mod m), gcd(a,m)=1
	// a^b ≡ a^(b mod φ(m) + φ(m)) (mod m), gcd(a,m)!=1 且 b>φ(m)
	// 模板题 https://www.luogu.com.cn/problem/P5091
	// 例题 https://codeforces.com/problemset/problem/615/D
	// https://cses.fi/problemset/task/1712
	exPhi := func(a, b int64, m int) int64 {
		phi := int64(calcPhi(m))
		if b >= phi {
			b = b%phi + phi
		}
		return pow(a, b, int64(m))
	}

	// 原根
	// https://en.wikipedia.org/wiki/Primitive_root_modulo_n
	// https://oi-wiki.org/math/primitive-root/
	// https://cp-algorithms.com/algebra/primitive-root.html
	// todo 原根&离散对数相关 https://www.luogu.com.cn/blog/command-block/yuan-gen-li-san-dui-shuo-xiang-guan
	// https://oeis.org/A033948 Numbers that have a primitive root (the multiplicative group modulo n is cyclic)
	//     The sequence consists of 1, 2, 4 and numbers of the form p^i and 2p^i, where p is an odd prime and i > 0
	// https://oeis.org/A046144 Number of primitive roots modulo n
	//    https://oeis.org/A010554 a(n) = phi(phi(n))
	// https://oeis.org/A008330 Number of primitive roots of n-th prime = phi(p-1)
	// https://oeis.org/A046145 Smallest primitive root modulo n, or 0 if no root exists
	// https://oeis.org/A001918 Smallest primitive root of n-th prime
	// https://oeis.org/A046146 Largest primitive root (<n) modulo n, or 0 if no root exists
	// https://oeis.org/A071894 Largest primitive root (<p) of n-th prime p
	// https://oeis.org/A056619 Smallest prime with primitive root of n or 0 if no such prime exists
	// https://oeis.org/A023049 Smallest prime > n having primitive root n, or 0 if no such prime exists
	// 模板题 https://www.luogu.com.cn/problem/U125141
	//
	// 返回 n 的最小的原根, n >= 2
	// 不存在时返回 -1
	// 由于有 phi(phi(n)) 个原根，密度足够大，最小原根可以很快找到，复杂度约为 O(n^0.25logn)
	primitiveRoot := func(n int) int {
		var gcd func(_, _ int) int
		var pow func(_, _, _ int) int

		if n != 2 && n != 4 {
			x := n
			if x&1 == 0 {
				x /= 2
			}
			if x&1 == 0 || len(primeDivisors(x)) > 1 {
				return -1
			}
		}

		pn := calcPhi(n)
		ps := primeDivisors(pn)
	o:
		for g := 1; ; g++ {
			if gcd(g, n) > 1 {
				continue
			}
			for _, p := range ps {
				if pow(g, pn/p, n) == 1 {
					continue o
				}
			}
			return g
		}
	}

	// 返回 n 的所有原根
	// n 没有原根时返回空切片
	// 模板题 https://www.luogu.com.cn/problem/P6091
	primitiveRootsAll := func(n int) []int {
		var gcd func(_, _ int) int

		rt0 := primitiveRoot(n)
		if rt0 < 0 {
			return nil
		}
		pn := calcPhi(n)
		ppn := calcPhi(pn)
		roots := make([]int, 0, ppn)
		for i, rtPow := 1, 1; len(roots) < ppn; i++ {
			rtPow = rtPow * rt0 % n
			if gcd(i, pn) == 1 {
				roots = append(roots, rtPow)
			}
		}
		sort.Ints(roots)
		return roots
	}

	/* 同余 逆元
	http://oeis.org/A006254 a(n) = 2^-1 mod p         Numbers n such that 2n-1 is prime
	http://oeis.org/A283419 3^-1 mod p
	http://oeis.org/A006093 (p-1)^-1 mod p     prime(n) - 1
	http://oeis.org/A040976 (p/2)^-1 mod p     prime(n) - 2
	http://oeis.org/A090938 Least multiple of n == 1 (mod prime(n))
	http://oeis.org/A090939 Least multiple of n == -1 (mod prime(n))
	http://oeis.org/A091185 a(n) = A090938(n)/n      n^-1 mod prime(n)
	*/

	// 二元一次不定方程（线性丢番图方程中的一种）https://en.wikipedia.org/wiki/Diophantine_equation
	// exgcd solve equation ax+by=gcd(a,b)
	// we have |x|<=b and |y|<=a in result (x,y)
	// https://cp-algorithms.com/algebra/extended-euclid-algorithm.html
	var exgcd func(a, b int64) (gcd, x, y int64)
	exgcd = func(a, b int64) (gcd, x, y int64) {
		if b == 0 {
			return a, 1, 0
		}
		gcd, y, x = exgcd(b, a%b)
		y -= a / b * x
		return
	}

	// a*x + b*y = c 的通解为
	// x = (c/g)*x0 + (b/g)*k
	// y = (c/g)*y0 - (a/g)*k
	// 其中 g = gcd(a,b)
	// x0 和 y0 是 ax+by=g 的一组特解（即 exgcd(a,b) 的返回值）
	//
	// 为方便讨论，这里要求输入的 a b c 必须为正整数
	// 注意若输入超过 1e9 可能要用高精
	// 返回：正整数解的个数（无解时为 -1，无正整数解时为 0）
	//      x 取最小正整数时的解 x1 y1，此时 y1 是最大正整数解
	//      y 取最小正整数时的解 x2 y2，此时 x2 是最大正整数解
	// 相关论文 THE NUMBER OF SOLUTIONS TO ax + by = n http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.376.403
	// 相关题目 https://www.luogu.com.cn/problem/P5656
	// 使非负解 x+y 尽量小 https://codeforces.com/problemset/problem/1244/C
	//    最简单的做法就是 min(x1+y1, x2+y2)
	// 需要转换一下符号 https://atcoder.jp/contests/abc186/tasks/abc186_e
	solveLinearDiophantineEquations := func(a, b, c int64) (n, x1, y1, x2, y2 int64) {
		g, x0, y0 := exgcd(a, b)

		// 无解
		if c%g != 0 {
			n = -1
			return
		}

		a /= g
		b /= g
		c /= g
		x0 *= c
		y0 *= c

		x1 = x0 % b
		if x1 <= 0 { // 若要求的是非负整数解，去掉等号
			x1 += b
		}
		k1 := (x1 - x0) / b
		y1 = y0 - k1*a

		y2 = y0 % a
		if y2 <= 0 { // 若要求的是非负整数解，去掉等号
			y2 += a
		}
		k2 := (y0 - y2) / a
		x2 = x0 + k2*b

		// 无正整数解
		if y1 <= 0 {
			return
		}

		// k 越大 x 越大
		n = k2 - k1 + 1
		return
	}

	// 关于 ax+by+cz=n 的解的个数（或者是三币种找零）
	// On the number of solutions of the Diophantine equation of Frobenius – General case https://core.ac.uk/download/pdf/14375587.pdf
	// The Number of Solutions to ax + by + cz = n and its Relation to Quadratic Residues https://cs.uwaterloo.ca/journals/JIS/VOL23/Binner/binner4.pdf
	// 上面这篇提出了一个 O(log max(a,b,c)) 的算法来求 N(a,b,c;n)
	// https://oeis.org/A001399 N(1,2,3;n) = round((n+3)^2/12)
	// https://oeis.org/A000115 N(1,2,5;n) = round((n+4)^2/20)
	// https://oeis.org/A008672 N(1,3,5;n) = round((n+3)*(n+6)/30)  =  floor((n^3+9n+30)/30)
	// https://oeis.org/A005044 N(2,3,4;n) = round(n^2/12)-floor(n/4)*floor((n+2)/4)      a(n) = round(n^2/48) if n is even; a(n) = round((n+3)^2/48) if n is odd
	// https://oeis.org/A025795 N(2,3,5;n)
	// https://oeis.org/A008680 N(3,4,5;n)

	// 任意非零模数逆元 ax ≡ 1 (mod m)
	// 返回最小正整数解
	// 模板题 https://www.luogu.com.cn/problem/P1082
	invM := func(a, m int64) int64 { _, x, _ := exgcd(a, m); return (x%m + m) % m }

	// 费马小定理求质数逆元
	// ax ≡ 1 (mod p)
	// x^-1 ≡ a^(p-2) (mod p)
	invP := func(a, p int64) int64 { return pow(a, p-2, p) }

	// 有理数取模
	// 模板题 https://www.luogu.com.cn/problem/P2613
	divM := func(a, b, m int64) int64 { return a * invM(b, m) % m }
	divP := func(a, b, p int64) int64 { return a * invP(b, p) % p }

	// 线性求逆元·其一
	// 求 1, 2, ..., p−1 mod p 的逆元
	// http://blog.miskcoo.com/2014/09/linear-find-all-invert
	// https://www.zhihu.com/question/59033693
	// 模板题 https://www.luogu.com.cn/problem/P3811
	initAllInv := func(p int) []int64 {
		inv := make([]int64, p)
		inv[1] = 1
		for i := 2; i < p; i++ {
			inv[i] = int64(p-p/i) * inv[p%i] % int64(p)
		}
		return inv
	}

	// 线性求逆元·其二（离线逆元）
	// 求 a1, a2, ..., an mod p 的逆元
	// 根据 ai^-1 ≡ Πai/ai * (Πai)^-1 (mod p)，求出 Πai 的前缀积和后缀积可以得到 Πai/ai，从而求出 ai^-1 mod p
	// https://zhuanlan.zhihu.com/p/86561431
	// 模板题 https://www.luogu.com.cn/problem/P5431
	calcAllInv := func(a []int64, p int64) []int64 {
		n := len(a)
		pre := make([]int64, n+1)
		pre[0] = 1
		for i, v := range a {
			pre[i+1] = pre[i] * v % p
		}
		invMulAll := invP(pre[n], p)
		suf := make([]int64, n+1)
		suf[n] = 1
		for i := len(a) - 1; i > 0; i-- { // i=0 不用求
			suf[i] = suf[i+1] * a[i] % p
		}
		inv := make([]int64, n)
		for i, pm := range pre[:n] {
			inv[i] = pm * suf[i+1] % p * invMulAll % p
		}
		return inv
	}

	// 模数两两互质的线性同余方程组 - 中国剩余定理 (CRT)
	// x ≡ bi (mod mi), bi 与 mi 互质且 Πmi <= 1e18
	// https://blog.csdn.net/synapse7/article/details/9946013
	// https://codeforces.com/blog/entry/61290
	// 模板题 https://www.luogu.com.cn/problem/P1495
	crt := func(B, M []int64) (x int64) {
		m := int64(1)
		for _, mi := range M {
			m *= mi
		}
		for i, mi := range M {
			Mi := m / mi
			_, inv, _ := exgcd(Mi, mi)
			x = (x + B[i]*Mi*inv) % m
		}
		x = (x + m) % m
		return
	}

	// 扩展中国剩余定理 (EXCRT)
	// ai * x ≡ bi (mod mi)
	// 解为 x ≡ b (mod m)
	// 有解时返回 (b, m)，无解时返回 (0, -1)
	// 推导过程见《挑战程序设计竞赛》P292
	// 注意乘法溢出的可能
	// 推荐 https://blog.csdn.net/niiick/article/details/80229217
	// todo 模板题 https://www.luogu.com.cn/problem/P4777 https://www.luogu.com.cn/problem/P4774
	// https://codeforces.com/contest/1500/problem/B
	excrt := func(A, B, M []int64) (x, m int64) {
		m = 1
		for i, mi := range M {
			a, b := A[i]*m, B[i]-A[i]*x
			d := gcd(a, mi)
			if b%d != 0 {
				return 0, -1
			}
			t := divM(b/d, a/d, mi/d)
			x += m * t
			m *= mi / d
		}
		x = (x%m + m) % m
		return
	}

	// 另一种写法，参考进阶指南 p.155
	// todo 待整理
	// excrt := func(a, m []int) (x int) {
	//	x = a[0]
	//	M := m[0]
	//	for i := 1; i < len(a); i++ {
	//		mi := m[i]
	//		c := (a[i] - x%mi + mi) % mi
	//		gcd, inv, _ := exgcd(M, mi)
	//		if c%gcd != 0 {
	//			return -1
	//		}
	//		c /= gcd
	//		mi /= gcd
	//		inv = inv * c % mi
	//		x += inv * M
	//		M *= mi
	//		x %= M
	//	}
	//	x = (x + M) % M
	//	return
	//}

	// 离散对数 - 小步大步算法 (BSGS)
	// a^x ≡ b (mod p)，a 和 p 互质
	// 有解时返回 x，无解时返回 -99，这样可以让 exBSGS 中的 +1 操作不影响无解的判断
	// 时间复杂度 O(√p)
	// 见进阶指南 p.155
	// todo https://www.luogu.com.cn/blog/command-block/yuan-gen-li-san-dui-shuo-xiang-guan
	//      http://blog.miskcoo.com/2015/05/discrete-logarithm-problem
	//      https://www.luogu.com.cn/blog/hzoiliuchang/shuo-lun-zhi-bsgs-suan-fa
	// 模板题 https://www.luogu.com.cn/problem/P3846
	babyStepGiantStep := func(a, b, p, k int64) int64 { // 非 exBSGS 下 k=1
		b %= p
		t := int64(math.Sqrt(float64(p))) + 1
		mp := map[int64]int64{}
		for j, x := int64(0), b; j < t; j++ {
			mp[b] = j
			x = x * a % p
		}
		a = pow(a, t, p)
		if a == 0 {
			if b == 0 {
				return 1
			}
			return -99
		}
		for i, x := int64(0), k; i < t; i++ {
			if j, ok := mp[x]; ok && i*t >= j {
				return i*t - j
			}
			x = x * a % p
		}
		return -99
	}

	// 拓展大步小步算法
	// a^x ≡ b (mod m)，a 和 m 不一定互质
	// https://zhuanlan.zhihu.com/p/132603308
	// 模板题 https://www.luogu.com.cn/problem/P4195
	var _exBSGS func(a, b, m, k int64) int64
	_exBSGS = func(a, b, m, k int64) int64 {
		if b == 1 {
			return 0
		}
		if a == 0 && b == 0 {
			return 1
		}
		g := gcd(a, m)
		if b%g > 0 {
			return -99
		}
		if g == 1 {
			return babyStepGiantStep(a, b, m, k%m)
		}
		return 1 + _exBSGS(a, b/g, m/g, k*a/g%m)
	}

	exBSGS := func(a, b, m int64) int64 {
		x := _exBSGS(a%m, b%m, m, 1)
		phiM := int64(calcPhi(int(m)))
		if x > phiM {
			x = x%phiM + phiM
		}
		return x
	}

	// 二次剩余 x^2 ≡ a (mod p)
	// 一个数 a，如果不是 p 的倍数且模 p 同余于某个数的平方，则称 a 为模 p 的二次剩余
	// https://en.wikipedia.org/wiki/Quadratic_residue
	// https://en.wikipedia.org/wiki/Cipolla%27s_algorithm
	// https://oi-wiki.org/math/quad-residue/
	// https://blog.csdn.net/doyouseeman/article/details/52033204
	// Tonelli-Shanks https://www.luogu.com.cn/blog/242973/solution-p5491
	// 模板题 https://www.luogu.com.cn/problem/P5491
	modSqrt := func(x, p int64) []int64 { // p 必须是奇素数
		if x == 0 {
			return []int64{0}
		}
		x0 := new(big.Int).ModSqrt(big.NewInt(x), big.NewInt(p))
		if x0 == nil {
			return nil
		}
		// 如果要求小的在前，注意判断
		return []int64{x0.Int64(), p - x0.Int64()}
	}

	// 判断 a 是否为模 p 的二次剩余，p 必须是奇素数
	// Jacobi 符号为 -1
	isQuadraticResidue := func(a, p int64) bool {
		return big.Jacobi(big.NewInt(a), big.NewInt(p)) < 0
	}

	// todo N 次剩余 / 高次同余方程 x^a ≡ b (mod p)
	// todo 模板题 https://www.luogu.com.cn/problem/P5668

	// https://oeis.org/A072994 Number of solutions to x^n ≡ 1 (mod n), 1<=x<=n
	// Least k > 0 such that the number of solutions to x^k == 1 (mod k) 1 <= x <= k is equal to n, or 0 if no such k exists https://oeis.org/A072995

	// http://oeis.org/A182865 Minimal number of quadratic residues
	// a(n) is the least integer m such that any nonzero square is congruent (mod n) to one of the squares from 1 to m^2
	// 把这题的 1000 改成 i，则至多需要枚举到 a(i) https://ac.nowcoder.com/acm/contest/6489/A

	// https://oeis.org/A000224 Number of distinct squares residues mod n
	// Multiplicative with a(p^e) = floor(p^e/6) + 2 IF p = 2 ELSE floor(p^(e+1)/(2p + 2)) + 1
	// https://oeis.org/A046530 Number of distinct cubic residues mod n

	/* 阶乘 组合数/二项式系数 */

	// https://en.wikipedia.org/wiki/Factorial
	// https://oeis.org/A000142
	// https://en.wikipedia.org/wiki/Stirling%27s_approximation
	// n! ~ √(2πn)*(n/e)^n
	factorial := []int{1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800 /*10!*/, 39916800, 479001600}

	// 等差数列的乘积转换成阶乘 https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e

	// https://oeis.org/A003070 a(n) = ceiling(log_2(n!))
	// https://oeis.org/A067850 Highest power of 2 not exceeding n!
	// https://oeis.org/A049606 Largest odd divisor of n!
	// https://oeis.org/A240533 a(n) = numerators of n!/10^n

	calcFactorial := func(n int) int64 {
		res := int64(1) % mod
		for i := 2; i <= n; i++ {
			res = res * int64(i) % mod
		}
		return res
	}
	// n 小于 1 时返回 1
	calcFactorialBig := func(n int) *big.Int {
		return new(big.Int).MulRange(1, int64(n))
	}

	initFactorial := func() {
		const mx int = 1e6
		F := [mx + 1]int64{1}
		for i := 1; i <= mx; i++ {
			F[i] = F[i-1] * int64(i) % mod
		}
	}

	// 阶乘模质数（质数较小）
	// 时间复杂度 O(plogn)
	// todo 待整理 https://cp-algorithms.com/algebra/factorial-modulo.html
	_factorial := func(n, p int) int {
		res := 1
		for ; n > 1; n /= p {
			if n/p&1 == 1 {
				res = res * (p - 1) % p
			}
			for i := 2; i <= n%p; i++ {
				res = res * i % p
			}
		}
		return res
	}

	// 双阶乘
	// https://en.wikipedia.org/wiki/Double_factorial

	// 偶阶乘
	// https://oeis.org/A000165 Double factorial of even numbers: (2n)!! = 2^n*n!
	calcEvenFactorialBig := func(n int) *big.Int {
		return new(big.Int).Lsh(new(big.Int).MulRange(1, int64(n)), uint(n))
	}

	// 奇阶乘
	// https://oeis.org/A001147 Double factorial of odd numbers: (2*n-1)!! = 1*3*5*...*(2*n-1) = A(2*n,n) / 2^n
	// 1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, 654729075, 13749310575, 316234143225, 7905853580625, ...
	// Number of ways to choose n disjoint pairs of items from 2*n items
	// Number of perfect matchings in the complete graph K(2n)
	// 相关题目 LC1359/双周赛20D 有效的快递序列数目 https://leetcode-cn.com/contest/biweekly-contest-20/problems/count-all-valid-pickup-and-delivery-options/
	// 奇阶乘模 2^64 http://acm.hdu.edu.cn/showproblem.php?pid=6481 https://www.90yang.com/hdu6481-a-math-problem/
	calcOddFactorialBig := func(n int) *big.Int {
		return new(big.Int).Rsh(new(big.Int).MulRange(int64(n+1), int64(2*n)), uint(n))
	}

	// https://oeis.org/A010786 Floor-factorial numbers: a(n) = Product_{k=1..n} floor(n/k)
	// 1, 2, 3, 8, 10, 36, 42, 128, 216, 600, 660, 3456, 3744, 9408, 18900, 61440, 65280, 279936, 295488, 1152000, 2116800, 4878720, 5100480, 31850496, 41472000, 93450240, 163762560, 568995840, 589317120, 3265920000, 3374784000
	// https://oeis.org/A309912 a(n) = Product_{p prime, p <= n} floor(n/p)

	// binomial(n, floor(n/2)) https://oeis.org/A001405
	// a(n) ~ 2^n / sqrt(π * n/2)
	// 从一个大小为 n 的集合的子集中随机选一个，选到 n/2 大小的子集的概率是 1 / sqrt(π * n/2)
	// Sperner's theorem says that this is the maximal number of subsets of an n-set such that no one contains another
	// EXTRA: https://oeis.org/A000984 Central binomial coefficients: binomial(2*n,n) = (2*n)!/(n!)^2
	// EXTRA: https://oeis.org/A100071 a(n) = n * A001405(n-1) = 1, 2, 6, 12, 30, 60, 140, 280, 630, 1260, ...
	//                                 a(n) = a(n-1) * n / floor(n/2)
	// EXTRA: https://oeis.org/A107373 a(n) = (n/2) * A001405(n-1) - 2^(n-2)
	combHalf := []int64{
		1, 1, 2, 3, 6, 10, 20, 35, 70, 126, // C(9,4)
		252, 462, 924, 1716, 3432, 6435, 12870, 24310, 48620, 92378, // C(19,9)
		184756, 352716, 705432, 1352078, 2704156, 5200300, 10400600, 20058300, 40116600, 77558760, // C(29,14)
		155117520, 300540195, 601080390, 1166803110, 2333606220, 4537567650, 9075135300, 17672631900, 35345263800, 68923264410, // C(39,19)
		137846528820, 269128937220, 538257874440, 1052049481860, 2104098963720, 4116715363800, 8233430727600, 16123801841550, 32247603683100, 63205303218876, // C(49,24)
		126410606437752, 247959266474052, 495918532948104, 973469712824056, 1946939425648112, 3824345300380220, 7648690600760440, 15033633249770520, 30067266499541040, 59132290782430712, // C(59,29)
		118264581564861424, 232714176627630544, 465428353255261088, 916312070471295267, 1832624140942590534, 3609714217008132870, 7219428434016265740, // C(66,33)
	}

	// 组合数/二项式系数
	// 不取模，仅适用于小范围的 n 和 k
	initComb := func() {
		const mx = 60
		C := [mx + 1][mx + 1]int64{}
		for i := 0; i <= mx; i++ {
			C[i][0], C[i][i] = 1, 1
			for j := 1; j < i; j++ {
				C[i][j] = C[i-1][j-1] + C[i-1][j]
			}
		}
	}

	// O(n) 预处理阶乘及其逆元，O(1) 求组合数
	{
		const mod int64 = 1e9 + 7
		const mx int = 2e6
		F := [mx + 1]int64{1}
		for i := 1; i <= mx; i++ {
			F[i] = F[i-1] * int64(i) % mod
		}
		pow := func(x, n int64) int64 {
			//x %= mod
			res := int64(1)
			for ; n > 0; n >>= 1 {
				if n&1 == 1 {
					res = res * x % mod
				}
				x = x * x % mod
			}
			return res
		}
		invF := [...]int64{mx: pow(F[mx], mod-2)}
		for i := mx; i > 0; i-- {
			invF[i-1] = invF[i] * int64(i) % mod
		}
		C := func(n, k int) int64 {
			if k < 0 || k > n {
				return 0
			}
			return F[n] * invF[k] % mod * invF[n-k] % mod
		}
		P := func(n, k int) int64 {
			if k < 0 || k > n {
				return 0
			}
			return F[n] * invF[n-k] % mod
		}

		// EXTRA: 卢卡斯定理
		// 注意初始化 F 和 invF 时 mx 取 mod-1
		// https://www.luogu.com.cn/problem/P3807
		// https://www.luogu.com.cn/problem/P7386
		var lucas func(n, k int64) int64
		lucas = func(n, k int64) int64 {
			if k == 0 {
				return 1
			}
			return C(int(n%mod), int(k%mod)) * lucas(n/mod, k/mod) % mod
		}

		// 可重组合 https://en.wikipedia.org/wiki/Combination#Number_of_combinations_with_repetition
		// 方案数 H(n,k)=C(n+k-1,k) https://oeis.org/A059481
		// 相当于把 k 个无区别的球放入 n 个有区别的盒子中，且允许空盒的方案数
		//		隔板法：把 n 个盒子当做 n-1 个隔板，这样相当于总共有 k+n-1个位置，从中选择 k 个位置放球，剩下的位置放隔板。这样就把 k 个球划分成了 n 份，放入对应的盒子中
		// NOTE: mx 要开两倍空间！
		H := func(n, k int) int64 { return C(n+k-1, k) }
		// 也相当于，给出长度和元素范围，求有多少种非降序列
		// 也可以理解成在长度和取值范围-1的格点上走单调路径
		H = func(range_, length int) int64 { return C(range_+length-1, length) }

		// 卡特兰数 Cn = C(2n,n)/(n+1) = C(2n,n)-C(2n,n-1)
		// https://en.wikipedia.org/wiki/Catalan_number
		// https://oeis.org/A000108
		// 从 n=0 开始：1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012, 742900, 2674440, 9694845, 35357670, 129644790
		// 所有在 n×n 格点中不越过对角线的单调路径的个数
		// Number of noncrossing partitions of the n-set (不相交握手问题) LC1259/双周赛13D https://leetcode-cn.com/contest/biweekly-contest-13/problems/handshakes-that-dont-cross/
		//
		// 将全部偶数提取一个 2，可得 (2n)! = 1*3*5*...*(2n-1) * (2^n) * (n!)
		// 故 C(2*n,n)/(n+1) = (2*n)!/(n!)/(n+1)! = 1*3*5*...*(2n-1)*(2^n)/(n+1)!
		// 又由于 n! 的 2 的因子个数 = n/2 + n/4 + ... + n/2^k <= n-1 当且仅当 n 为 2^k 时取到等号
		// 对比分子分母的 2 的因子个数，可以得出如下结论：
		//     当且仅当 n+1 为 2^k 时，卡特兰数为奇数
		//
		// EXTRA: 高维的情况 https://loj.ac/p/6051
		Catalan := func(n int) int64 { return F[2*n] * invF[n+1] % mod * invF[n] % mod }
		Catalan = func(n int) int64 { return new(big.Int).Rem(new(big.Int).Div(new(big.Int).Binomial(int64(2*n), int64(n)), big.NewInt(int64(n+1))), big.NewInt(mod)).Int64() }

		// 某些组合题可能用到
		pow2 := [mx + 1]int64{1}
		for i := 1; i <= mx; i++ {
			pow2[i] = pow2[i-1] << 1 % mod
		}

		_ = []interface{}{C, P, H, Catalan}
	}

	// 适用于 n 巨大但 k 或 n-k 较小的情况
	comb := func(n, k int) int64 {
		if k > n-k {
			k = n - k
		}
		var a, b int64 = 1, 1
		for i := 1; i <= k; i++ {
			a = a * int64(n) % mod
			n--
			b = b * int64(i) % mod
		}
		return divP(a, b, mod)
	}

	// 另类组合数求法
	{
		const mod int64 = 1e9 + 7
		var n, k int64
		// 注意当 n 为负数时，可能会算出非 0 的结果，这种情况要特判
		// 当 0 <= n < k 时结果为 0
		_ = new(big.Int).Binomial(n, k).Int64() // small
		_ = new(big.Int).Rem(new(big.Int).Binomial(n, k), big.NewInt(mod)).Int64()
		_ = int64(math.Round(math.Gamma(float64(n+1)) / math.Gamma(float64(k+1)) / math.Gamma(float64(n-k+1))))
	}

	// 扩展卢卡斯
	// todo https://blog.csdn.net/niiick/article/details/81064156
	// https://blog.csdn.net/skywalkert/article/details/52553048
	// https://blog.csdn.net/skywalkert/article/details/104681101
	// https://cp-algorithms.com/combinatorics/binomial-coefficients.html
	// 模板题 https://www.luogu.com.cn/problem/P4720
	// 古代猪文 https://www.luogu.com.cn/problem/P2480

	// 斯特林数（斯特林轮换数，斯特林子集数）
	// https://en.wikipedia.org/wiki/Stirling_number
	// https://oi-wiki.org/math/stirling/
	// todo 斯特林数的四种求法 https://www.luogu.com.cn/blog/command-block/si-te-lin-shuo-zong-jie
	// todo https://www.luogu.com.cn/blog/xzc/zu-ge-shuo-xue-hu-si-te-lin-shuo
	// https://blog.csdn.net/ACdreamers/article/details/8521134
	// Stirling numbers of the first kind, s(n,k) https://oeis.org/A008275
	//    将 n 个元素排成 k 个非空循环排列的方法数
	//    s(n,k) 的递推公式： s(n,k)=(n-1)*s(n-1,k)+s(n-1,k-1), 1<=k<=n-1
	//    边界条件：s(n,0)=0, n>=1    s(n,n)=1, n>=0
	//    todo https://www.luogu.com.cn/problem/P5408
	//         https://www.luogu.com.cn/problem/P5409
	// Stirling numbers of the second kind, S2(n,k) https://oeis.org/A008277
	//    将 n 个元素拆分为 k 个非空集的方法数
	//    S2(n, k) = (1/k!) * Σ{i=0..k} (-1)^(k-i)*binomial(k, i)*i^n.
	//    S2(n,k) 的递推公式：S2(n,k)=k*S2(n-1,k)+S2(n-1,k-1), 1<=k<=n-1
	//    边界条件：S(n,0)=0, n>=1    S(n,n)=1, n>=0
	//    https://www.luogu.com.cn/problem/P5395
	//    todo https://www.luogu.com.cn/problem/P5396
	//    https://oeis.org/A019538 n 个位置，每个位置填 [1,k] 之间的数，要求每个数字至少出现一次 => k!*S2(n,k)
	// Generalized Stirling numbers: a(n) = n! * Sum_{k=0..n-1} (k+1)/(n-k) https://oeis.org/A001705
	// Unsigned Stirling numbers of first kind: s(n+1,2): a(n+1) = (n+1)*a(n) + n! https://oeis.org/A000254
	// todo 斯特林数，斯特林反演初探 https://www.yijan.co/si-te-lin-shu-si-te-lin-fan-yan-chu-tan/
	// todo https://codeforces.com/contest/1278/problem/F 洛谷有艹标算的题解
	stirling2 := func(n, k int) int64 {
		s2 := make([][]int64, n+1)
		for i := range s2 {
			s2[i] = make([]int64, n+1)
		}
		s2[0][0] = 1
		for i := 1; i <= n; i++ {
			for j := 1; j <= i; j++ {
				s2[i][j] = (s2[i-1][j-1] + int64(j)*s2[i-1][j]) % mod
			}
		}
		return s2[n][k]
	}

	// 第二类斯特林数·行
	// https://www.luogu.com.cn/problem/P5395
	stirling2RowPoly := func(n int) poly {
		F := make([]int64, n+1)
		F[0] = 1
		for i := 1; i <= n; i++ {
			F[i] = F[i-1] * int64(i) % P
		}
		invF := make(poly, n+1)
		invF[n] = _pow(F[n], P-2)
		for i := n; i > 0; i-- {
			invF[i-1] = invF[i] * int64(i) % P
		}
		a := make(poly, n+1)
		b := make(poly, n+1)
		for i, v := range invF {
			if i&1 == 0 {
				a[i] = v
			} else {
				a[i] = P - v
			}
			b[i] = _pow(int64(i), n) * v % P
		}
		return a.conv(b)[:n+1]
	}

	// 贝尔数：基数为 n 的集合的划分方法数 https://oeis.org/A000110
	// https://en.wikipedia.org/wiki/Bell_number
	// 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 4213597, 27644437, 190899322, 1382958545, ...
	// B(n+1) = Sum_{k=0..n} C(n,k)*B(k)
	// B(n) = Sum_{k=1..n} Stirling2(n,k)
	bell := func(n int) (res int64) {
		s2 := make([][]int64, n+1)
		for i := range s2 {
			s2[i] = make([]int64, n+1)
		}
		s2[0][0] = 1
		for i := 1; i <= n; i++ {
			for j := 1; j <= i; j++ {
				s2[i][j] = (s2[i-1][j-1] + int64(j)*s2[i-1][j]) % mod
			}
		}
		for _, v := range s2[n][1:] {
			res += v
		}
		return res % mod
	}

	// 贝尔数的多项式求法
	// https://blog.csdn.net/a_forever_dream/article/details/106489066
	// https://www.luogu.com.cn/problem/P5748
	bellPoly := func(n int) poly {
		F := make([]int64, n+1)
		F[0] = 1
		for i := 1; i <= n; i++ {
			F[i] = F[i-1] * int64(i) % P
		}
		invF := make(poly, n+1)
		invF[n] = _pow(F[n], P-2)
		for i := n; i > 1; i-- { // 注意为了计算下面的 exp，invF[0] = 0
			invF[i-1] = invF[i] * int64(i) % P
		}

		b := invF.exp()
		for i, v := range b {
			b[i] = v * F[i] % P
		}
		return b
	}

	// 欧拉数 https://oeis.org/A008292
	// https://en.wikipedia.org/wiki/Eulerian_number
	// T(n, k) = k * T(n-1, k) + (n-k+1) * T(n-1, k-1), T(1, 1) = 1
	// T(n, k) = Sum_{j=0..k} (-1)^j * (k-j)^n * C(n+1, j)
	// todo 浅谈欧拉数 https://www.luogu.com.cn/blog/Karry5307/eulerian-numbers

	// 莫比乌斯函数 mu https://oeis.org/A008683
	// https://oi-wiki.org/math/mobius/#_11
	// 线性筛 https://oi-wiki.org/math/sieve/#_9
	// 前缀和 https://oeis.org/A002321 Mertens's function 梅滕斯函数
	//    https://en.wikipedia.org/wiki/Mertens_function
	//    |a(n)| = O(x^(1/2 + eps))
	//    零点 https://oeis.org/A028442
	sieveMu := func() {
		const mx int = 1e6
		mu := [mx + 1]int{1: 1}
		primes := []int{}
		vis := [mx + 1]bool{}
		for i := 2; i <= mx; i++ {
			if !vis[i] {
				mu[i] = -1
				primes = append(primes, i)
			}
			for _, p := range primes {
				v := p * i
				if v > mx {
					break
				}
				vis[v] = true
				if i%p == 0 {
					mu[v] = 0
					break
				}
				mu[v] = -mu[i]
			}
		}
	}

	// 莫比乌斯反演（岛娘推荐！https://zhuanlan.zhihu.com/p/133761303）
	// todo https://oi-wiki.org/math/mobius/

	//

	// 数论分块/除法分块/整除分块
	// https://oi-wiki.org/math/mobius/#_3
	//     https://oeis.org/A006218
	//     a(n) = Σ{k=1..n} floor(n/k)
	//          = n * (log(n) + 2*gamma - 1) + O(sqrt(n))
	//          也就是说平均每项的贡献约为 log(n)
	//     a(n) = Σ{k=1..n} floor(n/k)
	//          = 2*(Σ{i=1..floor(sqrt(n))} floor(n/i)) - floor(sqrt(n))^2
	//          因此 a(n) % 2 == floor(sqrt(n)) % 2
	//     a(n) 前缀和 = Sum_{k=1..n-1} Sum_{i=1..n-1} floor(k/i) https://oeis.org/A078567
	// 恒等式 n%i = n-(n/i)*i
	// ∑n/i https://www.luogu.com.cn/problem/P1403 n=1e18 的做法见 https://www.luogu.com.cn/problem/SP26073
	// ∑k%i 代码见下面的 floorLoopK https://www.luogu.com.cn/problem/P2261
	// ∑(n/i)*(n%i) https://ac.nowcoder.com/acm/contest/9005/C
	// todo https://codeforces.com/contest/1202/problem/F
	// ∑∑(n%i)*(m%j) 代码见下面的 floorLoop2 https://www.luogu.com.cn/problem/P2260
	//
	// EXTRA: 一些另类的求和
	// https://oeis.org/A116477 a(n) = Sum_{1<=k<=n, gcd(k,n)=1} floor(n/k)
	//                          sum{k|n} a(k) = sum{k=1 to n} d(k) = https://oeis.org/A006218
	// https://oeis.org/A013939 a(n) = Sum_{k = 1..n} floor(n/prime(k)) = omega(n!)
	//
	// EXTRA: n/k (k=1..n) 的不同数字的个数 https://oeis.org/A055086
	//        = floor(sqrt(4*n+1)) - 1 (注意用 int64)
	floorLoop := func(n int64) (sum int64) {
		for l, r := int64(1), int64(0); l <= n; l = r + 1 {
			h := n / l
			r = n / h
			w := r - l + 1
			sum += h * w
		}
		return
	}

	// ∑x/i, i in [low,up]
	// https://codeforces.com/contest/1485/problem/C
	floorLoopRange := func(low, up, x int64) (sum int64) {
		min := func(a, b int64) int64 {
			if a < b {
				return a
			}
			return b
		}
		for l, r := low, int64(0); l <= up; l = r + 1 {
			h := x / l
			if h == 0 {
				break
			}
			r = min(x/h, up)
			w := r - l + 1
			sum += h * w
		}
		return
	}

	// 余数求和
	// ∑k%i (when k=n its ∑n%i)
	// = ∑k-(k/i)*i
	// = n*k-∑(k/i)*i
	// 对于 [l,r] 范围内的 i，k/i 不变，此时 ∑(k/i)*i = (k/i)*∑i = (k/i)*(l+r)*(r-l+1)/2
	// https://www.luogu.com.cn/problem/P2261
	floorLoopK := func(n, k int64) int64 {
		min := func(a, b int64) int64 {
			if a < b {
				return a
			}
			return b
		}
		sum := n * k
		for l, r := int64(1), int64(0); l <= n; l = r + 1 {
			h := k / l
			if h > 0 {
				r = min(k/h, n)
			} else {
				r = n
			}
			w := r - l + 1
			s := (l + r) * w / 2
			sum -= h * s
		}
		return sum
	}

	// 二维整除分块
	// Σ{i=1..min(n,m)} floor(n/i)*floor(m/i)
	// https://www.luogu.com.cn/blog/command-block/zheng-chu-fen-kuai-ru-men-xiao-ji
	// todo ∑∑(n%i)*(m%j) 模积和 https://www.luogu.com.cn/problem/P2260
	floorLoop2 := func(n, m int64) (sum int64) {
		min := func(a, b int64) int64 {
			if a < b {
				return a
			}
			return b
		}
		for l, r := int64(1), int64(0); l <= min(n, m); l = r + 1 {
			hn, hm := n/l, m/l
			r = min(n/hn, m/hm)
			w := r - l + 1
			sum += hn * hm * w
		}
		return
	}

	// 杜教筛 - 积性函数前缀和
	// 复杂度 O(n^(2/3))
	// 【推荐】https://blog.csdn.net/weixin_43914593/article/details/104229700 算法竞赛专题解析（4）：杜教筛--以及积性函数的前世今生
	// https://www.luogu.com.cn/blog/command-block/du-jiao-shai
	// http://baihacker.github.io/main/
	// The prefix-sum of multiplicative function: the black algorithm http://baihacker.github.io/main/2020/The_prefix-sum_of_multiplicative_function_the_black_algorithm.html
	// The prefix-sum of multiplicative function: Dirichlet convolution http://baihacker.github.io/main/2020/The_prefix-sum_of_multiplicative_function_dirichlet_convolution.html
	// The prefix-sum of multiplicative function: powerful number sieve http://baihacker.github.io/main/2020/The_prefix-sum_of_multiplicative_function_powerful_number_sieve.html
	// todo 浅谈一类积性函数的前缀和 + 套题 https://blog.csdn.net/skywalkert/article/details/50500009
	// 模板题 https://www.luogu.com.cn/problem/P4213
	// todo ∑∑i*j*gcd(i,j) https://www.luogu.com.cn/problem/solution/P3768
	sieveDu := func() {
		const mx int = 1e6
		phi := [mx + 1]int64{1: 1}
		mu := [mx + 1]int{1: 1}
		primes := []int{}
		vis := [mx + 1]bool{}
		for i := 2; i <= mx; i++ {
			if !vis[i] {
				phi[i] = int64(i - 1)
				mu[i] = -1
				primes = append(primes, i)
			}
			for _, p := range primes {
				v := p * i
				if v > mx {
					break
				}
				vis[v] = true
				if i%p == 0 {
					phi[v] = phi[i] * int64(p)
					mu[v] = 0
					break
				}
				phi[v] = phi[i] * int64(p-1)
				mu[v] = -mu[i]
			}
		}
		for i := 0; i < mx; i++ {
			phi[i+1] += phi[i]
			mu[i+1] += mu[i]
		}

		cachePhi := map[int]int64{}
		var sumPhi func(int) int64
		sumPhi = func(n int) int64 {
			if n <= mx {
				return phi[n]
			}
			if s := cachePhi[n]; s > 0 {
				return s
			}
			m := int64(n)
			res := m * (m + 1) / 2
			for l, r := int64(2), int64(0); l <= m; l = r + 1 {
				h := m / l
				r = m / h
				res -= (r - l + 1) * sumPhi(int(h))
			}
			cachePhi[n] = res
			return res
		}

		cacheMu := map[int]int{}
		var sumMu func(int) int
		sumMu = func(n int) int {
			if n <= mx {
				return mu[n]
			}
			if s, has := cacheMu[n]; has {
				return s
			}
			res := 1
			for l, r := 2, 0; l <= n; l = r + 1 {
				h := n / l
				r = n / h
				res -= (r - l + 1) * sumMu(h)
			}
			cacheMu[n] = res
			return res
		}
	}

	// Min25 筛 - 积性函数前缀和
	// https://zhuanlan.zhihu.com/p/60378354
	// https://oi-wiki.org/math/min-25/
	// todo 模板题 https://www.luogu.com.cn/problem/P5325

	// 一篇新论文，复杂度为 O((nlogn)^(3/5))
	// Summing μ(n): a faster elementary algorithm
	// https://arxiv.org/pdf/2101.08773.pdf

	//

	// 埃及分数 - 不同的单位分数的和 (IDA*)
	// https://www.luogu.com.cn/problem/UVA12558
	// 贪婪算法：将一项分数分解成若干项单分子分数后的项数最少，称为第一种好算法；最大的分母数值最小，称为第二种好算法
	// https://en.wikipedia.org/wiki/Egyptian_fraction
	// https://oeis.org/A006585 number of solutions
	// https://oeis.org/A247765 Table of denominators in the Egyptian fraction representation of n/(n+1) by the greedy algorithm
	// https://oeis.org/A100678 Number of Egyptian fractions in the representation of n/(n+1) via the greedy algorithm
	// https://oeis.org/A100695	Largest denominator used in the Egyptian fraction representation of n/(n+1) by the greedy algorithm
	//
	// 		埃尔德什-施特劳斯猜想（Erdős–Straus conjecture）https://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93Straus_conjecture

	/* 斐波那契数列 F(n) https://oeis.org/A000045
	性质 https://oi-wiki.org/math/fibonacci/
	https://oeis.org/A000071 F(n) 前缀和 = F(n+2)-1
	∑k=[0,n]C(n,k)F(k) = F(2n)
	https://oeis.org/A007598 F^2(n)    a(n) = 2*a(n-1) + 2*a(n-2) - a(n-3), n > 2. a(0)=0, a(1)=1, a(2)=1
	                                  a(n) = (F(n)*F(n+4)-3*F(n)*F(n+1))/2
	https://oeis.org/A001690 补集
	https://oeis.org/A022307 F(n) 的不同的质因子个数
	https://oeis.org/A001175 N(m) = {F}%m 的周期    Pisano periods / Pisano numbers https://en.wikipedia.org/wiki/Pisano_period
	                         N(m) = LCM(N(p1^e1), ..., N(pk^ek))
	https://oeis.org/A060305 N(p) = {F}%p 的周期 https://blog.csdn.net/acdreamers/article/details/10983813
	https://oeis.org/A003893 F(n)%10
	https://oeis.org/A001605 使 F(n) 为质数的 n
	https://oeis.org/A191797 C(F(n), 2)
	https://oeis.org/A001177 Fibonacci entry points: a(n) = least k >= 1 such that n divides Fibonacci number
	异或和 F(n) 1,0,2,1,4,12,1,20,54,1,88,200,33,344,826,225,1756,3268,7313,1788
	定义 f(m) 为最小的满足 F(i)+F(j) ≡ 0 (mod m) 的 i (j<=i)，f(m) 大概是 O(√m) 的

	其他相关序列
	https://oeis.org/A000213 Tribonacci numbers: a(n)=a(n-1)+a(n-2)+a(n-3) with a(0)=a(1)=a(2)=1
	https://oeis.org/A000931 Padovan sequence (or Padovan numbers): a(n)=a(n-2)+a(n-3) with a(0)=1, a(1)=a(2)=0
	*/

	// https://oeis.org/A195264 Iterate x -> A080670(x) (replace x with the concatenation of the primes and exponents in its prime factorization)
	// starting at n until reach 1 or a prime; or -1 if a prime is never reached
	// https://www.zhihu.com/question/48612677/answer/487252829

	_ = []interface{}{
		primes, primes10, primes10_,
		sqCheck, cubeCheck, sqrt, cbrt, bottomDiff,
		gcd, gcdPrefix, gcdSuffix, lcm, frac, countDifferentSubsequenceGCDs, floorSum,
		isPrime, sieve, sieveEuler, sieveEulerTemplate, factorize, primeDivisors, powerOfFactorialPrimeDivisor, primeExponentsCountAll,
		divisors, divisorPairs, doDivisors, doDivisors2, oddDivisorsNum, maxSqrtDivisor, divisorsAll, primeFactorsAll, lpfAll, distinctPrimesCountAll,
		calcPhi, initPhi, sievePhi, exPhi,
		primitiveRoot, primitiveRootsAll,
		exgcd, solveLinearDiophantineEquations, invM, invP, divM, divP, initAllInv, calcAllInv,
		crt, excrt,
		babyStepGiantStep, exBSGS,
		modSqrt, isQuadraticResidue,
		factorial, calcFactorial, calcFactorialBig, initFactorial, _factorial, calcEvenFactorialBig, calcOddFactorialBig, combHalf, initComb, comb,
		stirling2, stirling2RowPoly,
		bell, bellPoly,
		sieveMu,
		floorLoop, floorLoopRange, floorLoopK, floorLoop2,
		sieveDu,
	}
}

/* 组合数学
https://en.wikipedia.org/wiki/Combination
https://en.wikipedia.org/wiki/Enumerative_combinatorics
二项式定理 https://en.wikipedia.org/wiki/Binomial_theorem

todo 组合数性质 | 二项式推论 https://oi-wiki.org/math/combination/#_13
一些常用组合恒等式的解释 https://www.zhihu.com/question/26094736
递推式 C(n-1, k-1) + C(n-1, k) = C(n, k)
上项求和 C(r, r) + C(r+1, r) + ... + C(n, r) = C(n+1, r+1)   相关题目 https://www.luogu.com.cn/problem/P7386
上式亦为 C(n, 0) + C(n+1, 1) + ... + C(n+m, m) = C(n+m+1, m) 相关题目 https://atcoder.jp/contests/abc154/tasks/abc154_f
范德蒙德恒等式 https://zh.wikipedia.org/wiki/%E8%8C%83%E5%BE%B7%E8%92%99%E6%81%92%E7%AD%89%E5%BC%8F
∑i=[0,k] C(n,i)*C(m,k-i) = C(n+m,k)
∑i>=n and k-i>=m C(i,n)*C(k-i,m) = C(k+1,n+m+1)    https://www.luogu.com.cn/blog/hanzhongtlx/ti-xie-0-1-trie
组合恒等式之万金油方法 https://zhuanlan.zhihu.com/p/25195967
∑i*C(n,i) = n*2^(n-1)

NOTE: 涉及到相邻的组合问题：可以考虑当前位置和左侧位置所满足的性质（例题 https://atcoder.jp/contests/abc167/tasks/abc167_e）

隔三组合数 https://oeis.org/A024493 https://oeis.org/A024494 https://oeis.org/A024495 C(n,0) + C(n,3) + ... + C(n,3[n/3])
隔四组合数 https://oeis.org/A038503 https://oeis.org/A038504 https://oeis.org/A038505 https://oeis.org/A000749

Tetrahedral (or triangular pyramidal) numbers: a(n) = C(n+2,3) = n*(n+1)*(n+2)/6 http://oeis.org/A000292
a(n) = Sum_{1<=i<=j<=n} j-i
a(n) = sum of all the possible products p*q where (p,q) are ordered pairs and p + q = n + 1
a(n) = 长度为 n 的字符串的所有子串长度之和

隔板法 https://zh.wikipedia.org/wiki/%E9%9A%94%E6%9D%BF%E6%B3%95
放球问题（总结得不错）https://baike.baidu.com/item/%E6%94%BE%E7%90%83%E9%97%AE%E9%A2%98
	扩展例题 https://codeforces.com/problemset/problem/893/E
todo 可重集排列组合 https://oi-wiki.org/math/combination/
todo https://codeforces.com/problemset/problem/451/E
不相邻的排列 https://oi-wiki.org/math/combination/#_10
错排 https://oeis.org/A000166 subfactorial numbers  a[n]=(n-1)*(a[n-1]+a[n-2])  https://zh.wikipedia.org/wiki/%E9%94%99%E6%8E%92%E9%97%AE%E9%A2%98
	https://oeis.org/A082491 n! * A000166(n)   a(n+2) = (n+2)*(n+1)*(a(n+1)+(n+1)*a(n))
	https://oeis.org/A000255 错排的比较对象的范围是 [1,n+1]  a(n) = n*a(n-1) + (n-1)*a(n-2), a(0) = a(1) = 1
	https://oeis.org/A000153 错排的比较对象的范围是 [1,n+2]  a(n) = n*a(n-1) + (n-2)*a(n-2), a(0) = 0, a(1) = 1
	https://oeis.org/A000261 错排的比较对象的范围是 [1,n+3]  a(n) = n*a(n-1) + (n-3)*a(n-2), a(1) = 0, a(2) = 1
	https://oeis.org/A001909 错排的比较对象的范围是 [1,n+4]  a(n) = n*a(n-1) + (n-4)*a(n-2), a(2) = 0, a(3) = 1
		https://atcoder.jp/contests/abc172/tasks/abc172_e
圆排列 https://zh.wikipedia.org/wiki/%E5%9C%86%E6%8E%92%E5%88%97
    Q(n,n) = (n-1)!

二阶递推数列通项 https://zhuanlan.zhihu.com/p/75096951
凯莱公式 Cayley’s formula: the number of trees on n labeled vertices is n^(n-2).
普吕弗序列 Prüfer sequence: 由树唯一地产生的序列
约瑟夫问题 Josephus Problem https://cp-algorithms.com/others/josephus_problem.html https://en.wikipedia.org/wiki/Josephus_problem
Stern-Brocot 树与 Farey 序列 https://oi-wiki.org/misc/stern-brocot/ https://cp-algorithms.com/others/stern_brocot_tree_farey_sequences.html

* 生成函数/母函数 *
https://en.wikipedia.org/wiki/Generating_function
https://oi-wiki.org/math/gen-func/intro/
todo 一些常见数列的生成函数推导 https://www.luogu.com.cn/blog/nederland/girl-friend
整数分拆 https://oeis.org/A000041 https://en.wikipedia.org/wiki/Partition_(number_theory)
https://oeis.org/A104513 The number of consecutive integers > 1 beginning with A104512(n), the sum of which equals n, or 0 if impossible.
						a(n)=0 iff n=2^k
https://oeis.org/A069283 将 n 分拆成至少两个连续整数的方法数 = n 的奇因子数 - 1
						见上面的 oddDivisorsNum 函数
https://oeis.org/A018819 Binary partition function: number of partitions of n into powers of 2
	相关题目 https://www.luogu.com.cn/problem/P6065 http://poj.org/problem?id=2229

	质数分拆
	https://oeis.org/A061358 Number of ways of writing n=p+q with p, q primes and p>=q
	https://oeis.org/A067187 Numbers that can be expressed as the sum of two primes in exactly one way
	https://oeis.org/A068307 number of partitions of n into a sum of three primes
	https://oeis.org/A071335 Number of partitions of n into a sum of at most three primes
	https://oeis.org/A023022 Number of partitions of n into two relatively prime parts
		a(n) = phi(n)/2 for n >= 3

https://oeis.org/A000404 Numbers that are the sum of 2 nonzero squares
https://oeis.org/A003325 Numbers that are the sum of 2 positive cubes

Maximum product of two integers whose sum is n https://oeis.org/A002620
Quarter-squares: floor(n/2)*ceiling(n/2). Equivalently, floor(n^2/4)
https://oeis.org/A024206 = A002620(n+1)-1 = floor((n-1)(n+3)/4)

	Maximal product of three numbers with sum n: a(n) = max(r*s*t), n = r+s+t https://oeis.org/A006501
	a(n) = floor(n/3)*floor((n+1)/3)*floor((n+2)/3)
	Expansion of (1+x^2) / ( (1-x)^2 * (1-x^3)^2 )

	Maximal product of four nonnegative integers whose sum is n https://oeis.org/A008233
	a(n) = floor(n/4)*floor((n+1)/4)*floor((n+2)/4)*floor((n+3)/4)

	...

	相关题目 https://codeforces.com/problemset/problem/1368/B

记 A = [1,2,...,n]，A 的全排列中与 A 的最大差值为 n^2/2 https://oeis.org/A007590
Maximum sum of displacements of elements in a permutation of (1..n)
For example, with n = 9, permutation (5,6,7,8,9,1,2,3,4) has displacements (4,4,4,4,4,5,5,5,5) with maximal sum = 40

https://oeis.org/A176127 The number of permutations of {1,2,...,n,1,2,...,n} with the property that there are k numbers between the two k's in the set for k=1,...,n
相关题目：《程序员的算法趣题》Q53 同数包夹

n married couples are seated in a row so that every wife is to the left of her husband
若不考虑顺序，则所有排列的个数为 (2n)!
考虑顺序可以发现，对于每一对夫妻来说，妻子在丈夫左侧的情况和在右侧的情况相同且不同对夫妻之间是独立的
因此每有一对夫妻，符合条件的排列个数就减半
所以结果为 a(n) = (2n)!/2^n https://oeis.org/A000680
或者见这道题目的背景 LC1359 https://leetcode-cn.com/problems/count-all-valid-pickup-and-delivery-options/

NxN 大小的对称置换矩阵的个数 http://oeis.org/A000085
这里的对称指仅关于主对角线对称
a[i] = (a[i-1] + (i-1)*a[i-2]) % mod
The number of n X n symmetric permutation matrices
Number of self-inverse permutations on n letters, also known as involutions; number of standard Young tableaux with n cells
Proof of the recurrence relation a(n) = a(n-1) + (n-1)*a(n-2):
	number of involutions of [n] containing n as a fixed point is a(n-1);
	number of involutions of [n] containing n in some cycle (j, n),
	where 1 <= j <= n-1, is (n-1) times the number of involutions of [n] containing the cycle (n-1 n) = (n-1)*a(n-2)
相关题目 https://ac.nowcoder.com/acm/contest/5389/C

The number of 3 X n matrices of integers for which the upper-left hand corner is a 1,
the rows and columns are weakly increasing, and two adjacent entries differ by at most 1
a(n+2) = 5*a(n+1) - 2*a(n), with a(0) = 1, a(1) = 4
https://oeis.org/A052913
相关题目 LC1411/周赛184D https://leetcode-cn.com/problems/number-of-ways-to-paint-n-x-3-grid/ https://leetcode-cn.com/contest/weekly-contest-184/

男厕问题 / 电话问题 https://oeis.org/A185456
Assume that the first person to use a bank of payphones selects one at the end,
and all subsequent users select the phone which puts them farthest from the current phone users.
U(n) is the smallest number of phones such that n may be used without any two adjacent phones being used
https://www.zhihu.com/question/278361000/answer/1004606685

https://oeis.org/A089934 Table T(n,k) of the number of n X k matrices on {0,1} without adjacent 0's in any row or column
https://oeis.org/A006506 上面这个 table 的对角线
    Number of configurations of non-attacking princes on an n X n board, where a "prince" attacks the four adjacent (non-diagonal) squares
    Also number of independent vertex sets in an n X n grid

https://oeis.org/A001224 The number of inequivalent ways to pack a 2 X n rectangle with dominoes
    If F(n) is the n-th Fibonacci number, then a(2n) = (F(2n+1) + F(n+2))/2 and a(2n+1) = (F(2n+2) + F(n+1))/2
    https://oeis.org/A060312

https://oeis.org/A001187 n 个节点的无向连通图的个数  Number of connected labeled graphs with n nodes
相关题目：https://www.acwing.com/problem/content/309/

韦德伯恩-埃瑟林顿数
https://oeis.org/A001190 Wedderburn-Etherington numbers: unlabeled binary rooted trees (every node has out-degree 0 or 2) with n endpoints (and 2n-1 nodes in all)
https://en.wikipedia.org/wiki/Wedderburn%E2%80%93Etherington_number
https://mathworld.wolfram.com/WeaklyBinaryTree.html 给出了如下公式：
    a(2n-1) = a(1)a(2n-2) + a(2)a(2n-3) + ... + a(n-1)a(n)
    a(2n)   = a(1)a(2n-1) + a(2)a(2n-2) + ... + a(n-1)a(n+1) + a(n)(a(n)+1)/2
https://oeis.org/A000598 Number of rooted ternary trees with n nodes
https://oeis.org/A268172 Binary-ternary Wedderburn-Etherington numbers
相关题目：《程序员的算法趣题》Q30 用插线板制作章鱼脚状线路

https://oeis.org/A003991 Multiplication table read by antidiagonals = (n-k+1)*k
https://oeis.org/A059036 = A003991(n, k) - 1

一些二进制的计数问题见 bits.go

CF 上的一些组合计数问题 http://blog.miskcoo.com/2015/06/codeforces-combinatorics-and-probabilities-problem

找出 50% 作弊者 https://codingcompetitions.withgoogle.com/codejam/round/000000000043580a/00000000006d1155
    讨论 https://codeforces.com/blog/entry/84822
*/
func combinatoricsCollection() {
	// 容斥原理 (PIE, the principle of inclusion and exclusion)
	// 参考《挑战程序设计竞赛》P296
	// https://codeforces.com/blog/entry/64625
	// https://ac.nowcoder.com/acm/contest/6219/C
	solveInclusionExclusion := func(a []int) (ans int64) {
		n := len(a)
		const mod int64 = 1e9 + 7 // 998244353
		for sub := uint(0); sub < 1<<n; sub++ {
			res := int64(0)
			for i, v := range a {
				if sub>>i&1 == 1 {
					// 视情况而定，有时候包含元素 i 表示考虑这种情况，有时候表示不考虑这种情况
					_ = v // do v...

				}
			}
			if bits.OnesCount(sub)&1 == 1 { // 某些题目是 0
				res = -res
			}
			ans += res // mod
		}
		ans = (ans%mod + mod) % mod // 注意最后变为非负
		return
	}

	_ = []interface{}{solveInclusionExclusion}
}

// 数值分析
// https://en.wikipedia.org/wiki/Numerical_analysis
func numericalAnalysisCollection() {
	type mathF func(x float64) float64

	// Simpson's 1/3 rule
	// https://en.wikipedia.org/wiki/Simpson%27s_rule
	// 证明过程 https://phqghume.github.io/2018/05/19/%E8%87%AA%E9%80%82%E5%BA%94%E8%BE%9B%E6%99%AE%E6%A3%AE%E6%B3%95/
	simpson := func(l, r float64, f mathF) float64 { h := (r - l) / 2; return h * (f(l) + 4*f(l+h) + f(r)) / 3 }

	// 不放心的话还可以设置一个最大递归深度 maxDeep
	// 15eps 的证明过程 http://www2.math.umd.edu/~mariakc/teaching/adaptive.pdf
	var _asr func(l, r, eps, A float64, f mathF) float64
	_asr = func(l, r, eps, A float64, f mathF) float64 {
		mid := l + (r-l)/2
		L := simpson(l, mid, f)
		R := simpson(mid, r, f)
		if math.Abs(L+R-A) <= 15*eps {
			return L + R + (L+R-A)/15
		}
		return _asr(l, mid, eps/2, L, f) + _asr(mid, r, eps/2, R, f)
	}

	// 自适应辛普森积分 Adaptive Simpson's Rule
	// https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
	// https://oi-wiki.org/math/integral/
	// https://cp-algorithms.com/num_methods/simpson-integration.html
	// 模板题 https://www.luogu.com.cn/problem/P4525 https://www.luogu.com.cn/problem/P4526
	asr := func(a, b, eps float64, f mathF) float64 { return _asr(a, b, eps, simpson(a, b, f), f) }

	//

	// 多项式插值
	// https://en.wikipedia.org/wiki/Polynomial_interpolation

	// 拉格朗日插值
	// 给定多项式上的 n 个点 (xi,yi)，求 f(k)
	// https://en.wikipedia.org/wiki/Lagrange_polynomial
	// https://oi-wiki.org/math/poly/lagrange/
	// 浅谈几种插值方法 https://www.luogu.com.cn/blog/zhang-xu-jia/ji-zhong-cha-zhi-fang-fa-yang-xie
	//
	// 模板题 https://www.luogu.com.cn/problem/P4781
	// todo https://www.luogu.com.cn/problem/P5667
	// 等幂和 https://codeforces.com/problemset/problem/622/F
	lagrangePolynomialInterpolation := func(xs, ys []int64, k int64) int64 {
		const mod = 998244353

		pow := func(x, n int64) int64 {
			x %= mod
			res := int64(1)
			for ; n > 0; n >>= 1 {
				if n&1 == 1 {
					res = res * x % mod
				}
				x = x * x % mod
			}
			return res
		}
		inv := func(a int64) int64 { return pow(a, mod-2) }
		div := func(a, b int64) int64 { return a % mod * inv(b) % mod }

		ans := int64(0)
		for i, xi := range xs {
			a, b := ys[i]%mod, int64(1)
			for j, x := range xs {
				if j != i {
					a = a * (k - x) % mod
					b = b * (xi - x) % mod
				}
			}
			ans += div(a, b)
		}
		ans = (ans%mod + mod) % mod
		return ans
	}

	_ = []interface{}{
		asr,
		lagrangePolynomialInterpolation,
	}
}