数学

模乘

快速乘

可能比 __int128 快
|a|,|b|<m，且m在longlong范围
编译环境支持64个二进制位精度的long doube（如x86-64平台的GNU G++）

LL mul(LL a, LL b, LL m) {
    u64 r = (u64)a * b - (u64)((LD)a / m * b + 0.5L) * m;
    return r < m ? r : r + m;
}

龟速乘

template<typename T>
T pow_mod(T a,T b,T m){
    T res=(m!=1);
    while(b){
        if(b%2)res=mul_mod(res,a,m);
        a=mul_mod(a,a,m);
        b/=2;
    }
    return res;
}

快速幂

longlong范围用fpl

LL fp(LL a, LL b, LL Mod) {
    LL res = (Mod != 1);
    for (; b; b >>= 1, a = a * a % Mod)
        if (b & 1)
            res = res * a % Mod;
    return res;
}

LL fpl(LL a, LL b, LL Mod) {
    LL res = (Mod != 1);
    for (; b; b >>= 1, a = mul(a, a, Mod))
        if (b & 1)
            res = mul(res, a, Mod);
    return res;
}

光速幂

给定 $(a,m)=1$ ，要多次询问 $f(b)=a^b\bmod m$

对 $a$ 的幂分块，预处理出 $a^0,a^1,a^2,a^3,...a^t,a^{2t},a^{3t},...,a^{t^2}$ ，常见的取法是 $t=\lceil \sqrt{\phi(m)} \rceil$ ，或 $t=\lceil \sqrt m \rceil$ ，不小于阶即可

有 $f(b)=f(b \mod{\phi(m)})=f(kt)*f(r)$ ，$O(1)$ 回答询问

最大公约数（欧几里得算法）

实际上 <algorithm> 有函数__gcd 在 cpp17 的 numeric 中有函数 std::gcd

template <typename T>
T gcd(T a, T b) {
	while (b){T t = b;b = a % b;a = t;}
	return a;
}

template <typename T>
T lcm(T a, T b) { return a / gcd(a, b) * b; }

template<typename T>
void exgcd(T a,T b,T&x,T&y,T&d){
	if (b == 0){ x = 1, y = 0, d = a;}
	else{
		exgcd(b, a % b, y, x, d);
		y -= (a / b) * x;
	}
}

binary-gcd

常数极小，longlong版本更换对应的64位函数

int bgcd(int a, int b) {
    int az = __builtin_ctz(a);
    int bz = __builtin_ctz(b);
    int z = min(az, bz);
    b >>= bz;
    while (a) {
        a >>= az;
        int diff = a - b;
        az = __builtin_ctz(diff);
        b = min(a, b), a = abs(diff);
    }
    return b << z;
}

乘法逆元

利用欧拉定理(费马小定理)通过快速幂求逆元代码最短

exgcd求逆元

需要 c++11

int getinv(int a, int p){
    int m = 0, n = 1, x = 1, y = 0, b = p;
    while (b){
        int q = a / b;
        tie(x, m) = make_tuple(m, x - q * m);
        tie(y, n) = make_tuple(n, y - q * n);
        tie(a, b) = make_tuple(b, a - q * b);
    }
    (x < 0) && (x += p);
    return x;
}

线性递推求逆元

模 $p$ 时，设正整数 $i$ 和 $p%i$ 均存在逆元，有

$$ i^{-1}=(p-\lfloor p/i\rfloor )*(p\mod\ i)^{-1} $$

可以 $O(n)$ 求出 $1..n$ 模固定质数 $p$ 的逆元，但是空间局部性不佳，对于 $n=10^7$ 数量级不如用对阶乘求逆的 $O(n+log P)$ 做法

线性同余方程组的合并（ex中国剩余定理）

用于合并同余方程组

$$ x\equiv a_i(mod\ m_i) $$

其中 $a_i\geq 0$ ，$m_i>0$ ，$2*lcm(m_1,...)\in long long$

前置算法 快速乘 和 exgcd

用 pair<LL,LL> 存 ${a_i,m_i}$

有解返回 True，并将第二个方程合并进第一个（得到最小非负解）；否则返回 False

注意题目求最小非负解还是正数解！！！

bool crt(pll &a, pll b){
	LL g, u, d;
	exgcd(a.second, b.second, u, d, g);
	d = b.first - a.first;if (d % g) return 0;
	const LL t = b.second / g;
	a.first += mul(u, d / g, t) * a.second; 
	a.second *= t;
	a.first %= a.second; 
    (a.first < 0) && (a.first += a.second); //注意题目要求最小非负解还是正数解
	return 1;
}

大指数运算（欧拉降幂）

$$ a^x=\begin{cases} a^{x\ mod \phi(m)}&(a,m)=1\\ a^x&(a,m)!=1且x< c\\ a^{(x-c)mod\ c+c}&(a,m)!=1且x\geq c\\ \end{cases} $$

其中 $c=m的质因子最高次数 +1$ ，是一个大小为 $O(log\ n)$ 的常数，为计算方便常用 $\phi(m)$ 代替它

只需在快速幂中特判，若 $b^c≥\phi(m)$ ，则返回 $b^{c\ mod\ \phi(m)+\phi(m)}$ ，否则返回 $b^{c\ mod\ \phi(m)}$

LL exfp(LL a,LL b,LL Mod){
    /*注意检查Mod<=1e9 且res初值为1*/
    bool f=0;LL res=1;
    for(;b;b>>=1){
        if(a>Mod){f=1;a%=Mod;}
        if(b&1){
            res=res*a;
            if(res>Mod){f=1;res%=Mod;}
        }
        a=a*a;
    }return f?res+Mod:res;
}

Miller-Rabin 素性检测

loj143 luoguP4718

$n<=10^{18}$

需要 mul fpl

http://miller-rabin.appspot.com/

对于 $n\leq 2^{64}$ 保证正确的7个base: 2, 325, 9375, 28178, 450775, 9780504, 1795265022

对于$ n \leq 3 * 10^{24}$，可以选 a = 2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41 （前 13 个质数）

bool isPrime(LL n) {
    if (n <= 3) return n > 1;
    if (n % 2 == 0 || n % 3 == 0) return 0;
    LL d = n - 1;
    while (d % 2 == 0) d /= 2;
    for (LL a : {2, 325, 9375, 28178, 450775, 9780504, 1795265022}) {
        if (a % n == 0) continue;
        LL t = d;
        LL u = fpl(a, t, n);
        while (t != (n - 1) && u != 1 && u != (n - 1)) {
            u = mul(u, u, n);
            t *= 2;
        }
        if (t % 2 == 0 && u != (n - 1)) return 0;
    }
    return 1;
}

__int128版本 51nod 质数检测V2 （陈旧的垃圾代码

// a b \in [0,m)
template<typename T>
T mul_mod(T a,T b,T m){
    T res=0;
    while(b){
        if(b%2){
            res+=a;
            if(res>m)res-=m;
        }
        a*=2;
        if(a>m)a-=m;
        b/=2;
    }
    return res;
}

template<typename T>
T pow_mod(T a,T b,T m){
    T res=(m!=1);
    while(b){
        if(b%2)res=mul_mod(res,a,m);
        a=mul_mod(a,a,m);
        b/=2;
    }
    return res;
}

template<typename T>
bool M_R(T a,const T&d,const int&t,const T&n){//注意a%n==0时应判定通过测试
    a%=n;if(a==0)return 1;
    a=pow_mod(a,d,n);
    if(a==1||a==(n-1))return 1;
    for(int i=1;i<t;++i){
        a=mul_mod(a,a,n);
        if(a==(n-1))return 1;
        if(a==1)return 0;
    }
    return 0;
}

template<typename T>
bool isPrime(T n){ 
    if(n==2||n==3)return 1;
    if(n<2||(n%2==0)||(n%3==0))return 0;
    T d=n-1;int t=0;
    while(d%2==0)d/=2,++t;
    srand(time(0));
    for(int i=0;i<8;++i)if(!M_R(rand()%(n-1)+1,d,t,n))return 0;
    return 1;
}

Pollard-Rho 分解质因数

需要 mul __gcd

返回 $n$ 的一个大于 $1$ 的真因子

不存在就死循环，寄！因此在调用前务必检查 $n$ 是大于 $1$ 的合数

时间复杂度猜想为 $O(n^{1/4})$

mt19937 mt(time(0));  //随机化
//Pollard_Rho  返回n的一个>1的约数  调用前务必判断素性
inline LL PR(LL n) {
    LL x = uniform_int_distribution<LL>(0, n - 1)(mt), s, t,
       c = uniform_int_distribution<LL>(1, n - 1)(mt);  //随机化
    for (int gol = 1; 1; gol <<= 1, s = t, x = 1) {
        for (int stp = 1; stp <= gol; ++stp) {
            LL u = mul(t, t, n) + c;
            t = (u >= n ? u - n : u);
            x = mul(x, abs(s - t), n);
            if ((stp & 127) == 0) {
                LL d = __gcd(x, n);
                if (d > 1) return d;
            }
        }
        LL d = __gcd(x, n);
        if (d > 1) return d;
    }
}

vector<LL> get_prime_divisor(LL n) {
    vector<LL> res; if (n <= 1) return res;
    queue<LL> fac; fac.push(n);
    while (fac.size()) {
        LL u = fac.front();
        fac.pop();
        if (isPrime(u)) res.push_back(u);
        else {
            LL d = PR(u);
            while (d == u) d = PR(u);
            while (u % d == 0) u /= d;
            fac.push(d);
            if (u > 1) fac.push(u);
        }
    }
    sort(res.begin(), res.end());
    unique(res.begin(), res.end());
    return res;
}

原根

只有 $2,4,p^a,2p^a$ 有原根，其中 $p$ 是奇质数

（王元）质数 $p$ 的最小正原根是 $O(p^{0.25+\epsilon})$ 的，因此可以暴力查找质数的原根

设 $g$ 为 $m$ 的原根，对任意 $x$ 满足 $(x,\phi(m))=1$ ，$g^x$ 也是原根，故原根共有 $\phi(\phi(m))$ 个

（原根判定定理）设 $(g,m)=1$ ，则 $g$ 是 $m$ 的原根 $\iff$ 对 $\phi(m)$ 的任意质约数 $p$ 有 $g^{\frac{\phi(m)}{p}}\neq 1(mod\ m)$

寻找质数原根

求任意质数 $p$ 的原根

前置算法 __gcd 快速乘 快速幂

需要对 $\phi(p)$ 分解，一般要配合 Pollard-rho

时间复杂度猜测为 $O(p^{0.25+\epsilon})$

int Primitive_Root(LL p){
    if(p<=3)return p-1;

    vector<LL>g=Pollard_Rho::get_prime_divisor(p-1);
    for(auto &it:g)it=(p-1)/it;

    for(int q=1;q<p;++q){
        bool t=1;
        for(const auto&it:g)if(fpl(q,it,p)==1){t=0;break;}
        if(t)return q;
    }
    return -1;
}

高次同余方程

以下 p 均为正整数

离散对数（BSGS/exBSGS）

给定 $a,b,p$ 求（最小非负） $x$ 满足 $$ a^x\equiv b\pmod p $$ 其中 $p\leq 10^9$ ，不必是质数

时间复杂度 $O(\sqrt p)$

前置算法 快速幂 map/unorded_map __gcd() 求逆元
- 或者自定义哈希表 typedef unordered_map<int,int,neal> HASH;

视 p 范围换用快速乘，然而由于复杂度，实则难以对一般的long long做BSGS ，此时多半是想法歪了或者性质没有观察充分

若多组询问 a,p 为常量且始终满足 (b,p)=1 则应复用哈希表 h 以降低常数
接口说明传入非负整数 a b 和正整数 p
- log() 要求 $(a,p)=1$ 无解返回-1 只求最小非负整数解
- operator() 无特殊要求无解返回-1 可以修改使得求最小正数解，见注释

typedef unordered_map<int,int,neal> HASH;

struct{
	int log(int a,int b,int p){
        if((b-1)%p==0)return 0;
        HASH h;const int t=(int)sqrt(p)+1;
		for(int j=0;j<t;++j,b=(LL)b*a%p)h[b]=j;
		int at=(a=fp(a,t,p));
		for(int i=1;i<=t;++i,a=(LL)a*at%p){
            auto it=h.find(a);
            if(it!=h.end())return i*t-(*it).second;
        }return -1;
	}
	int operator()(int a,int b,int p){//(接口)
        if((b-1)%p==0)return 0;//求正整数解则注释掉
        if((b-a)%p==0)return 1;
        a%=p;b%=p;
        const int _a=a,_b=b,_p=p;
		int d=__gcd(a,p),ax=1,t=0;
        while(d>1){
            if(b%d!=0)break;
            b/=d;p/=d;ax=(LL)ax*(a/d)%p;++t;
            d=__gcd(a,p);
        }
        for(int i=2,at=(LL)_a*_a%_p;i<=t;++i,at=(LL)at*_a%_p)if(at==_b)return i;
        if(d>1)return -1;
        int res=log(a,(LL)b*getinv(ax,p)%p,p);//res为-1说明无解 有解则非负
        return res<0?-1:res+t;//求正整数解则要修改
	}
}BSGS;

模任意数开任意次根/N次剩余

给定 $a,b,m$ ，求（最小非负） $x$ 满足 $$ x^a\equiv b\pmod m $$

~~对于一般的情形，暂时不会~~

特殊情形一

若 $(b,m)=1$ 且 $(a,\phi(m))=1$ 且易求得 $\phi(m)$

则存在非负整数 $u,v$ 满足 $au-\phi(m)v=1$，又 $b^{\phi(m)}\equiv1 \pmod m$

得 $b^{au}\equiv b^{au-\phi(m)v}\equiv b \pmod m$

即 $x=b^u$

特殊情形二

若 $(b,m)=1$ ，且易求得 $m$ 的原根 $g$ 和 $\phi(m)$ ，则步骤如下：

计算离散对数 $y=\log_gb$ （时间瓶颈所在）

将问题转化为 $x^a=(g^{x'})^a\equiv g^y(mod\ m)$

等价于求解 $x'$ 满足线性同余方程 $ax'\equiv y(mod\ \phi(m))$

模奇质数开平方根/立方根

二次剩余判定定理

设 $n$ 是奇质数 $p$ 的一个既约剩余类

$$ n^{\frac{p-1}{2}}\bmod p=\begin{cases}1 \qquad & n是二次剩余 \ -1 \qquad & n是非二次剩余 \end{cases} $$

求平方根（Cipolla算法）

求最小非负整数 $x$ 满足 $$ x^2\equiv n \pmod p $$ 其中 $p$ 为不超过 $10^9$ 的质数， $n$ 为非负整数

期望复杂度 $O(2\log\ p)$

若有两根则另一根为 $p-x$
若无解返回 -1
对 $10^{18}$ 以内的 $p$ 需用快速乘见注释
接口 int ans1 = Square_root(p)(n);

struct Square_root{
	typedef int DAT; //操作数类型
	const DAT P;
	DAT I2;
	Square_root(DAT p = 1) : P(p) {}
	DAT mul(DAT a, DAT b) const{ return (LL)a * b % P;} // DAT=int
	/*DAT mul(DAT a, DAT b) const{
		u64 r = (u64)a * b - (u64)((LD)a / P * b + 0.5L) * P;
		return r < P ? r : r + P;
	}*/ // DAT=long long

	#define X first
	#define Y second
	typedef pair<DAT, DAT> pii;
	pii mul(pii a, pii b) const { return {
		(mul(a.X, b.X) + mul(mul(I2, a.Y), b.Y)) % P, 
		(mul(a.X, b.Y) + mul(a.Y, b.X)) % P};}
	
	template <class T>
	T pow(T a, DAT b, T x){
		for (; b; b /= 2, a = mul(a, a)) if (b & 1)x = mul(x, a);
		return x;
	}

	DAT operator()(DAT n){ 
		if ((n %= P) <= 1) return n;
		if (pow(n, (P - 1) / 2, (DAT)1) == P - 1) return -1;
		DAT a;
		do a = rand(); // 随机方法可能寄
		while (pow(I2 = (mul(a, a) - n + P) % P, (P - 1) / 2, (DAT)1) == 1); 
		DAT x = pow(pii{a, 1}, (P + 1) / 2, {1, 0}).X;
		return min(x, P - x);
	}
	#undef X
	#undef Y
};

立方根（三次剩余）

求一个非负整数 $x$ 满足 $$ x^3\equiv n\pmod p $$ 其中 $p$ 是质数

复杂度期望 $O(9log\ p)$

参考https://blog.csdn.net/skywalkert/article/details/52591343

$n=0$ 或 $p<=3$ 需要特判
若 $p\equiv 2\pmod 3$，则 $x\equiv n^{(2p-1)/3} \pmod p$ 是唯一解
若 $p \equiv 1\pmod 3$，则 $\epsilon =(\sqrt{-3}-1)/2$ 是三次单位根，其满足 $\epsilon ^3 \equiv 1\pmod p$
- 如果找到一个解 $x$，则另外两个解为 $x\epsilon$，$x\epsilon^2$

接口 int ans2 = Cube_root(n, p)(); 返回其中一个非负整数解，无解返回 -1

换 long long 见注释

struct Cube_root{
    typedef int DAT;//操作数类型
    const DAT n,p;Cube_root(DAT _n=1,DAT _p=1):n(_n%_p),p(_p){};
    struct Z3{DAT x,y,z;};
    DAT mul(DAT a,DAT b)const{return (LL)a*b%p;}  //DAT=int
    //DAT mul(DAT a,DAT b)const{u64 r=(u64)a*b-(u64)((LD)a/p*b+0.5L)*p;return r<p?r:r+p;}  //DAT=long long
    Z3 mul(Z3 a,Z3 b)const{return(Z3){
        (mul(a.x,b.x)+mul((mul(a.y,b.z)+mul(a.z,b.y))%p,n))%p,
        ((mul(a.x,b.y)+mul(a.y,b.x))%p+mul(mul(a.z,b.z),n))%p,
        ((mul(a.x,b.z)+mul(a.y,b.y))%p+mul(a.z,b.x))%p
    };}
    template<class T>T pow(T a,DAT b,T x)	
    	{for(;b;b/=2,a=mul(a,a))if(b&1)x=mul(x,a);return x;}
    DAT operator()(){ 
        if(n==0||p<=3)return n;
        if(p%3==2)return pow(n,(2*p-1)/3,(DAT)1);
        if(pow(n,(p-1)/3,(DAT)1)!=1)return -1;
        Z3 r;
        do r=pow((Z3){rand(),rand(),rand()},(p-1)/3,(Z3){1,0,0});
        while(r.x!=0||r.y==0||r.z!=0); // 随机寄了就换均匀随机
        return pow(r.y,p-2,(DAT)1);
    }
};

线性筛

~~这也要板子?~~

int ntp[N],pri[N],tot;
inline void linear_sieve(int n){
    ntp[1]=1;for(int i=2;i<=n;++i)ntp[i]=0;tot=0;//初始化 一般可省
    for(int i=2;i<=n;++i){
        if(!ntp[i])pri[++tot]=i;
        for(int j=1;j<=tot;++j){
            const LL nex=(LL)i*pri[j];if(nex>n)break;//n<=1e6则不必开LL
            ntp[nex]=1;
            if(i%pri[j]==0)break;
        }
    }
}

欧拉函数 $\phi(x)$

求 $\phi(x)$ 的值可以归约为找$x$的质因子，于是可以 $O(\sqrt{x})$ 求点值

template <typename T>
inline T phi(T x) {
	T res = x;
	for (T i = 2; i * i <= x; ++i)
		if ((x % i) == 0) {
			res = res / i * (i - 1);
			while ((x % i) == 0) x /= i;
		}
	if (x > 1) res = res / x * (x - 1);
	return res;
}

或者 $O(n)/O(n log log n)/O(nlogn)$ 筛出 $\phi (1..n)$

整除分块

对于给定的正整数 $n$ 和 $i$，使得 $$ \lfloor \frac{n}{i}\rfloor =\lfloor \frac{n}{j}\rfloor $$ 成立的最大正整数（$i<=j<=n$）$j$ 的值为 $$ j=\lfloor \frac{n}{\lfloor \frac{n}{i}\rfloor}\rfloor $$

即 $\lfloor \frac{n}{i}\rfloor$ 所在块的右端点的标号为 $\lfloor \frac{n}{\lfloor \frac{n}{i}\rfloor}\rfloor$

莫比乌斯变换（数论）/狄利克雷卷积

P5495

给定 $a[1..n]$ ，求 $b[1..n]$ 满足 $$ b_k=\sum_{i|k}a_i $$

先筛出所有质数，视每个质数为一个维度，做高维前缀和，本质上与 Sm of Subset DP 相同

时间复杂度 $O(n\log\log n)$

for (int i = 1; i <= tot; ++i)
    for (int j = 1; pri[i] * j <= n; ++j) 
        a[pri[i] * j] += a[j];

莫比乌斯反演（数论）

狄利克雷卷积 $*$ 的代数性质
- 交换律、结合律
- 对加法(函数值逐项相加)分配律
- 单位元存在 $\epsilon=[n==1]$
- 对 $f(1)\neq0$ 的 $f$ 存在逆元 $g$
- 两个积性函数的狄利克雷卷积仍是积性函数
$\mathbb 1*\mu=\epsilon$
$\mathbb 1 *\phi=id$
设 $f$ 和 $g$ 是数论函数(N->R)，且当 $n>M$ 时有 $f(n)=g(n)=0$，则
- $f(n)=\sum_{n|m}g(m)\iff g(n)=\sum_{n|m}\mu(\frac m n)f(m)$

杜教筛

给定数论函数 $f$ 和一个较大的 $n$，要计算其前缀和 $$ S(n)=\sum_{i=1}^nf(i) $$ 考虑构造一个函数 $g$ ，令其与 $f$ 作狄利克雷卷积，有

$$ \sum_{i=1}^n(f*g)(i)=\sum_{i=1}^n\sum_{xy=i}f(x)g(y)=\sum_{y=1}^ng(y)S(\lfloor\frac{n}{y}\rfloor) $$

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

其中，$f*g$ 的前缀和与 $g$ 的区间和要求可以快速计算（不高于计算 $f$ 的复杂度即可，$O(1)$ 推式子或者也用杜教筛算都行）

时间复杂度为 $O(n^{3/4})$

如果能线性预处理出 $S(1..n)$ 的前 $n^{2/3}$ 项，可以将复杂度平衡至 $O(n^{2/3})$

关于记忆化：

对于单组询问，可以不用 map/unordered_map ，在线性预处理的前提下直接索引 $\lfloor \frac{N}{x} \rfloor$ 中的 $x$ 即可
对于多组询问，用unordered_map 可以利用不同询问的信息，可能更快

特别注意对 `n` 数论分块时，最后一次循环之前，执行 `i=j+1` 之后，此时 $i=n+1$，如果 `n` 刚好为数据类型的上界，会导致溢出！！！

rqy大佬的【模板】杜教筛（Sum）代码

其中 1300，maxm，maxn 分别是 $N$，$N^{2/3}$，$N^{1/3}$ 的最大值

typedef long long LL;
const int maxn = 2147483647;
const int maxm = 2000000;

LL S1[maxm], S2[1300];
bool vis[1300];
int N;

LL S(int n) {
	if (n < maxm) return S1[n];
	int x = N / n;  // 如果存在某个 x 使得 n = floor(N / x)，
					// 选 x = floor(N / n) 一定可以。
					// 且对于 x > N^{1/3} 的情形会直接查线性表
	if (vis[x]) return S2[x];
	vis[x] = true;
	LL &ans = S2[x];
	ans = (LL)n * (n + 1) / 2; // 对 (f*g)(n) = n 求和
    //在加强了数据后，对于n=(2^31)-1，下面这里i会溢出成负2^32，导致接下来发生除0 RE，应该开 long long
	for (int i = 2, j; i <= n; i = j + 1) { 
        j = n / (n / i);
        ans -= (j - i + 1) * S(n / i);
  	}
  	return ans;
}

本人完整代码

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

typedef long long LL;

const int N=2000009;

LL phi[N];
int miu[N],pri[N],tot;
int n;
bool vis[2021];
pair<LL,int> S2[2021];

pair<LL,int> S(int t){
    if(t<=2000000)return {phi[t],miu[t]};
    int x=n/t;if(vis[x])return S2[x];
    vis[x]=1;
    pair<LL,int>&res=S2[x];
    res.first=t*(t+1ll)/2;res.second=1;
    for(LL i=2,j;i<=t;i=j+1){
        j=t/(t/i);
        pair<LL,int> u=S(t/i);
        res.first-=(j-i+1)*u.first;
        res.second-=(j-i+1)*u.second;
    }
    return res;
}


inline void solve(int T){
    scanf("%d",&n);
    memset(vis,0,sizeof vis);
    pair<LL,int> ans=S(n);
    printf("%lld %d\n",ans.first,ans.second);
}

signed main(){
    miu[1]=1;phi[1]=1;
    for(int i=2;i<=2000000;++i){
        if(phi[i]==0)pri[++tot]=i,phi[i]=i-1,miu[i]=-1;
        for(int j=1;j<=tot;++j){
            const LL m=i*pri[j];if(m>2000000)break;
            phi[m]=phi[i]*(pri[j]-1);
            if(i%pri[j]){
                miu[m]=-miu[i];
            }else{
                phi[m]+=phi[i];
                break;
            }
        }
        miu[i]+=miu[i-1];
        phi[i]+=phi[i-1];
    }
    int t;scanf("%d",&t);
    for(int i=1;i<=t;++i)solve(i);
    return 0;
}

Min_25筛

给定积性函数 $f(x)$ 和较大的 $n$ 求 $$ S(n)=\sum_{i=1}^nf(i) $$

令第 $i$ 小的质数为 $p_i$ （特殊地， $p_0=1$），先线性筛出 $\sqrt n$ 以内的所有质数

Step.1

先将原来的积性函数 $f(x)$ 拆成若干个完全积性函数之和

对于特定的函数比如素数数量/素数前缀和，可以用long long存，在dp过程中不取模，算完再取模

对完全积性函数 $f(x)$ 求 $$ g(n)=\sum_{p\in prime}f(p) $$

考虑DP，设 $g(n,j)=\sum_{i=1}^n{[i是质数或其最小质因子>p_j]f(i)}$

则需要的结果是 $g(n,k)$ ，其中 $p_k$ 为 $\sqrt n$ 以内最大质数

从 $j-1$ 到 $j$ 有转移 $$ g(n,j)=g(n,j-1)-f(p_j)(\ g(\lfloor \frac{n}{p_j} \rfloor,j-1)-g(p_{j-1},j-1)\ ) $$

初始值 $g(n,0)=\sum_{i=2}^nf(i)$ ，一般而言这是个自然数已知幂的和，可以 $O(1)$ 算

时间复杂度 $O(\frac{n^{3/4}}{log n})$

只用到了形如 $\lfloor \frac{n}{x} \rfloor$ 的约 $2\sqrt n$ 处的DP值（形如 $p_{j-1}$ 的点已经包含在内），对于 $x\leq \sqrt n$，映射到下标 x ，对于 $x>\sqrt n$ ，映射到下标 n/x
可以滚动掉第二维

const int N=2e5; // sqrt(n)  

LL n;
int p[N],tot; 
... // 输入n，预处理sqrt(n)以内质数p[1..tot]，
    // 注意一定记得 p[0]=1

int id1[N],id2[N],m=0; // 存离散化的 n/x
LL v[2*N],g[2*N];

int getid(LL x){ // 找每个点值的下标
    return x<N?id1[x]:id2[n/x];
}
    
for(LL i=1,j;i<=n;i=j+1){ // 初始化
    LL u=n/i; j=n/u;
    v[++m]=u;
    if(u<N)id1[u]=m;else id2[n/u]=m;
    g[m]=...; // g[n/i][0] = sum(f(2..n))
}

for(int j=1;j<=tot;++j){
	for(int i=1;i<=m && 1ll*p[j]*p[j]<=v[i];++i){
        g[i]=g[i]-(.../*f(p_j)*/)*( g[getid(v[i]/p[j])] - g[getid(p[j-1])] );
	}
}

求 $n$ 以内素数数量

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

const int N=316333;

LL n; //1e11
int m;

int id1[N],id2[N],cnt;
LL v[2*N+1],g[2*N+1];
int getid(LL u){
    if(u<N)return id1[u];else return id2[n/u];
}

int pri[N],tot; bool ntp[N];

signed main(){
    scanf("%lld",&n);m=sqrt(n)+1;

    pri[0]=1;
    for(int i=2;i<=m;++i){
        if(!ntp[i])pri[++tot]=i;
        for(int j=1;j<=tot;++j){
            const int nex=pri[j]*i;if(nex>m)break;
            ntp[nex]=1;
            if(i%pri[j]==0)break;
        }
    }

    for(LL i=1,j;i<=n;i=j+1){
        LL t=n/i;j=n/t;
        v[++cnt]=t;
        if(t<N)id1[t]=cnt;else id2[j]=cnt;
        g[cnt]=t-1;
    }

    for(int j=1;j<=tot;++j){
        const LL pj=pri[j],lim=pj*pj;
        const LL gpre=g[getid(pri[j-1])];
        for(int i=1;i<=cnt&&lim<=v[i];++i)
            g[i]+=gpre-g[getid(v[i]/pj)];
    }

    printf("%lld\n",g[1]);
    return 0;
}

Step.2

此时已求出积性函数 $f(x)$ 在质数处的点值和 $g(n)=\sum_{p\in prime}f(p)$ ，且需要满足 $f(p^a)$ 可以快速计算

法一

令 $S(n,j)=\sum[i的最小质因子>p_j]f(i)$ ，则有 $$ S(n,j)=g(n)-g(p_j)+\sum_{j<k \and p_k\leq \sqrt n}\sum_{e=1}^{p_k^e \leq n} f(p_k^e)(S(\lfloor n/p_k^e \rfloor,k)+[e\neq 1]) $$

特殊地，当 $p_j\geq n$ 时 $S(n,j)=0$

考虑 $S(\lfloor n/p_k^e \rfloor,k)$ 在 $\lfloor n/p_k^e \rfloor \geq p_k$ 时才有值，也可变形为 $$ S(n,j)=g(n)-g(p_j)+\sum_{j<k \and p_k\leq \sqrt n}\sum_{e=1}^{p_k^e * p_k \leq n}(f(p_k^e)S(\lfloor n/p_k^e \rfloor,k)+ f(p_k^e * p_k)) $$

无需记忆化直接暴力递归，最终答案即为 $S(n,0)+f(1)$

时间复杂度为 $O(n^{1-\epsilon})$ ，当 $n\leq 10^{13}$ 时跑得比 $O(\frac{n^{3/4}}{log n})$ 快

P5325 【模板】Min_25筛

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

const int N=100009;
const int M=1000000007,i2=(M+1)/2,i6=(M+1)/6;

LL n;
int m;

int id1[N],id2[N],cnt;
LL v[2*N+1];
int g1[2*N+1],g2[2*N+1];
int getid(LL u){
    if(u<N)return id1[u];else return id2[n/u];
}

int pri[N],tot;
bool ntp[N];

LL f(LL pa){return pa%M*((pa-1+M)%M)%M;}

int S(LL t,int j){
    if(pri[j]>=t)return 0;
    int u=getid(t),v=getid(pri[j]);
    int res=(LL(g2[u])-g1[u]-g2[v]+g1[v]+M+M)%M;
    for(int k=j+1;k<=tot&&LL(pri[k])*pri[k]<=t;++k){
        for(LL pa=pri[k];pa*pri[k]<=t;pa*=pri[k]){
            res=(res+f(pa)*S(t/pa,k)+f(pa*pri[k]))%M;
        }
    }
    return res;
}

signed main(){
    scanf("%lld",&n);m=sqrt(n)+1;

    pri[0]=1;
    for(int i=2;i<=m;++i){
        if(!ntp[i])pri[++tot]=i;
        for(int j=1;j<=tot;++j){
            int nex=pri[j]*i;if(nex>m)break;
            ntp[nex]=1;
            if(i%pri[j]==0)break;
        }
    }

    for(LL i=1,j;i<=n;i=j+1){
        LL t=n/i;j=n/t;
        v[++cnt]=t;
        if(t<N)id1[t]=cnt;else id2[j]=cnt;
        g1[cnt]=(t+2)%M*((t-1)%M)%M*i2%M;
        g2[cnt]=(t%M)*((t+1)%M)%M*((2*t+1)%M)%M*i6%M-1;
        if(g2[cnt]<0)g2[cnt]+=M;
    }

    for(int j=1;j<=tot;++j){
        const LL pj=pri[j],lim=pj*pj;
        const int pre=getid(pri[j-1]);
        for(int i=1;i<=cnt&&lim<=v[i];++i){
            const int now=getid(v[i]/pj);
            g1[i]=(g1[i]-pj*(g1[now]-g1[pre]+M)%M+M)%M;
            g2[i]=(g2[i]-pj*pj%M*(g2[now]-g2[pre]+M)%M+M)%M;
        }
    }

    printf("%d\n",(S(n,0)+1)%M);
    return 0;
}

法二

$O(\frac{n^{3/4}}{log n})$ ，代码较长，由于常数大一般跑得不如法一快，暂没学

组合数的计算

计算 $C(n,m)\bmod p$

$m$ 非常小，$p$ 为质数，此时 $n^{m\over}/m!$ 项数很小
若 $p$ 是固定的线性大小的质数，$O(n)$ 预处理阶乘逆元
若 $n,m$ 很大但 $p$ 是很小的质数，根据 $Lucas$ 定理 $$ C_n^m\equiv C_{\lfloor \frac np\rfloor}^{\lfloor \frac mp\rfloor}*C_{n\mod p}^{m\mod p}\pmod p $$
若 $n,m$ 很大但 $p$ 是线性大小的任意数，用预处理分治的思想求 $C_n^m mod\ P$ ：

基于模数的预处理分治（exlucas）

luoguP4720

接口 exlucas(p).C(n,m)

预处理 $O(分解p耗时+\sum p_i^{a_i})$

单次询问 $O(\log^2 n)$ ，用光速幂可以做到 $O(\log n)$

不要求P为质数

struct exlucas{
    struct dat{
        const int p,pt,phi,t;
        vector<int>r;
        int Pow(int x,int y)const{
            int res=1;
            while(y>0){
                if(y%2!=0)res=LL(res)*x%pt;
                y/=2;x=LL(x)*x%pt;
            }
            return res;
        }
        dat(int _p,int _pt,int _t):p(_p),pt(_pt),phi(_pt/_p*(_p-1)),t(_t),r(pt){
            r[0]=1;
            for(int i=1;i<pt;++i)if(i%p==0)r[i]=r[i-1];
            else r[i]=r[i-1]*LL(i)%pt;
        }
        int f(LL n)const{
            if(n<p)return r[n];
            return Pow(r.back(),n/pt%phi)*LL(r[n%pt])%pt*f(n/p)%pt;
        }
        LL g(LL n)const{
            LL res=0;
            while(n>0)res+=n/p,n/=p;
            return res;
        }
        int C(LL n,LL m){
            LL v=g(n)-g(m)-g(n-m);
            if(v>=t)return 0;
            return LL(f(n))*Pow(f(m),phi-1)%pt*Pow(f(n-m),phi-1)%pt*Pow(p,v)%pt;
        }
    };
    vector<dat>mp;
    vector<LL>u;
    const int p;
    exlucas(int _p):p(_p){
        for(int i=2;i*i<=_p;++i)if(_p%i==0){
            int m=1,cnt=0;
            while(_p%i==0)_p/=i,m*=i,++cnt;
            mp.emplace_back(i,m,cnt);
        }
        if(_p>1)mp.emplace_back(_p,_p,1);
        u.resize(mp.size());
        for(int i=0;i<int(u.size());++i){
            u[i]=1;
            for(int j=0;j<int(mp.size());++j)if(i!=j){
                u[i]=u[i]*mp[j].pt%p*mp[i].Pow(mp[j].pt,mp[i].phi-1)%p;
            }
        }
    }
    int C(LL n,LL m){
        if(n<m||m<0)return 0;
        int res=0;
        for(int i=0;i<int(u.size());++i)res=(res+u[i]*mp[i].C(n,m))%p;
        return res;
    }
};

直线下整点数（类欧几里得算法）

给定 $n\geq 0,a\geq 0,c>0,b$，求 $$ f(a,b,c,n)=\sum_{i=0}^n\lfloor \frac{ai+b}{c} \rfloor $$

不妨设 $0\leq a,b<c$ （否则直接 $\bmod c$ 调整）， $m=\lfloor \frac{ai+b}{c} \rfloor \geq 1$ （否则直接return），

换一维枚举贡献即得 $f(a,b,c,n)=nm-f(c,c-b-1,a,m-1)$

相当于 $a$ 与 $c$ 在辗转相除，时间复杂度 $O(\log \max(a,c))$

LL f(LL a, LL b, LL c, LL n) {
    assert(c > 0 && n >= 0 && a >= 0);
    LL u = a / c; a -= u * c;
    LL v = b / c; b -= v * c; 
    if (b < 0) {b += c; --v;} assert(b >= 0);
    LL res = u * (n * (n + 1) / 2) + (n + 1) * v;
    if (a == 0) return res + b / c * (n + 1);
    LL m = (a * n + b) / c; if (m == 0) return res;
    return n * m - f(c, c - b - 1, a, m - 1) + res;
}

洛谷模板题

f=∑[(ai+b)/c] g=∑i[(ai+b)/c] h=∑[(ai+b)/c]^2 i=0..n a,b,n∈N c∈N*

struct dat{LL f, g, h;};
const LL i2 = 499122177, i3 = 332748118, M = 998244353; //预处理出模M意义下2和3的逆元
dat f(LL a, LL b, LL c, LL n){
    LL ac = a / c, bc = b / c;
    LL n2 = (n * (n + 1) % M) * i2 % M, n3 = n2 * (2ll * n + 1) % M * i3 % M;
    dat res = {
        (n2 * ac % M + (n + 1) * bc % M) % M, 
        (ac * n3 % M + bc * n2 % M) % M, 
        (ac * ac % M * n3 % M + 
            bc * bc % M * (n + 1) % M + ac * bc % M * n2 % M * 2ll) % M};
    a %= c; b %= c; if (a == 0)return res;
    LL m = (a * n + b) / c;
    dat p = f(c, c - b - 1, a, m - 1);
    LL fc = (n * m % M - p.f + M) % M, gc = (n2 * m % M - i2 * (p.f + p.h) % M + M) % M;
    return{(res.f + fc) % M,    (res.g + gc) % M, 
    (res.h + 2ll * (bc * fc % M + ac * gc % M) % M + 
        n * m % M * m % M - 2ll * p.g - p.f + 3ll * M) % M};}

~~loj模板题不会~~

高斯整数

四则运算与辗转相除

struct G {
    long long  r, i;
    G(long long  re = 0, long long  im = 0) : r(re), i(im) {}
    long long  norm() const { return r * r + i * i; }
    
    G operator+(const G& b) const { return G(r + b.r, i + b.i); }
    G operator-(const G& b) const { return G(r - b.r, i - b.i); }
    G operator*(const G& b) const {
        return G(r * b.r - i * b.i, i * b.r + r * b.i);
    }
    G operator/(const G& b) const {
        long long  l = b.norm();
        return G(round(1.0 * (r * b.r + i * b.i) / l),
                 round(1.0 * (i * b.r - r * b.i) / l));
    }
    G operator%(const G& b) const { return (*this) - (*this) / b * b; }
    
    void get() { scanf("%lld%lld", &r, &i); }
    void write() const { printf("%lld %lld ", r, i); }
};

G gcd(G a, G b) { return (b.norm() == 0) ? a : gcd(b, a % b); }

$4k+1$ 型质数的分解（二次剩余）

给定质数 $p$ 满足 $p\equiv 1\pmod 4$ ，解出一对共轭复数 $a\pm bi$ 满足 $(a+bi)(a-bi)=a^2+b^2=p$

先求正整数 $x$ 满足 $x\equiv -1 \pmod p$，由于 $-1$ 是 $4k+1$ 型质数的二次剩余，一定有解

即 $p|(x^2+1)$ ，由于 $x^2+1=(x+i)(x-i)$，计算 $d=gcd(p,x+i)$，则 $d$ 即为 $a\pm bi$ 其中之一

平分分解计数

设 $n$ 为正整数，记 $4*f(n)$ 为二维平面的圆 $x^2+y^2=n$ 上的整点数，则 $f(n)$ 是积性函数： $$ f(p^a)=\begin{cases}a+1\qquad &p=4k+1\ [a是偶数] \qquad &p=4k+3\ 1\qquad &p=2 \end{cases} $$ 或者对任意正整数 $d$ 定义完全积性函数 $\chi(d)$ ：

$$ \chi(d)=\begin{cases}1\qquad &d=4k+1\ -1 \qquad &d=4k+3\ 0\qquad &d=2k \end{cases} $$

易证 $$ \sum_{d|n}\chi(d)=f(n) $$

特别地，对于 $x^2+y^2=n^2$ 上的整点数 $4*g(n)$ 有 $$ g(p^a)=\begin{cases}2a+1\qquad &p=4k+1\ 1 \qquad &else \end{cases} $$ $\sum f(1..n)$ 与 $\sum g(1..n)$ 都可以用min25筛计算

数值积分

自适应辛普森算法
复杂度 O(玄学)
luogu p4542

struct IG{
    typedef double Func(double);
    const Func*f;IG(const Func&g):f(&g){}
    double simpson(double l,double r)const{
        double mid=(l+r)/2;
        return (r-l)*(f(l)+4*f(mid)+f(r))/6;
    }
    double find(double l,double r,const double&EPS,double res)const{
        double mid=(l+r)/2;
        double fl = simpson(l, mid), fr = simpson(mid, r);
        if (abs(fl + fr - res) <= 15 * EPS)return fl + fr + (fl + fr - res) / 15; 
        return find(l, mid, EPS / 2, fl) +find(mid, r, EPS / 2, fr);  
    }
    double operator()(double l,double r,double EPS=1e-8)const{return find(l,r,EPS,simpson(l,r));}
};

struct Integration{
    typedef double Func(double);
    const Func*f;Integration(const Func&g):f(&g){}
    typedef pair<double,double> pdd;
    pdd add(pdd a,pdd b)const{
        double mid=(a.first+b.first)/2;
        return (pdd){mid,f(mid)};
    }
    #define simpson(p1,p2,p3) (((p3).first-(p1).first)*((p1).second+4*(p2).second+(p3).second)/6)
    double find(pdd p1,pdd p3,pdd p5,double EPS,double res,int dep)const{
        pdd p2=add(p1,p3),p4=add(p3,p5);
        double fl=simpson(p1,p2,p3),fr=simpson(p3,p4,p5),d=(fl+fr-res)/15;
        if(abs(d)<=EPS&&dep<0)return fl+fr+d; 
        return find(p1,p2,p3,EPS/2,fl,dep-1)+find(p3,p4,p5,EPS/2,fr,dep-1);  
    }
    double operator()(double l,double r,double EPS=1e-6/*精度*/,int dep=12/*最小递归深度*/)const{
        pdd p1(l,f(l)),p3(r,f(r)),p2=add(p1,p3);
        return find(p1,p2,p3,EPS,simpson(p1,p2,p3),dep);
    }
    #undef simpson
};

日期操作

用于跳转的常量

const LL year_1[2]={365, 366};
const LL year_400=1460097;
const LL m_day[13]={(LL)0x3f3f3f, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

辅助函数

inline bool isLeap(LL t){return (t % 400 == 0) ||((t % 4 == 0) && (t % 100));}
inline bool pick(LL a, LL b){return ((isLeap(a) && b <= 2) ||(isLeap(a + 1) && b > 2));}
inline LL dayThisMonth(LL y, LL m){return m_day[m] + isLeap(y) * (m == 2);}

日期和整数的一一对应

LL可以改成int

struct MY_DATE{
    LL year, month, day;
    MY_DATE(LL y = 2021, LL m = 1, LL d = 1) : year(y), month(m), day(d){};
    LL p(MY_DATE op = {0, 0, 0}){//日期转换为整数
        LL y = year - op.year, m = month - op.month, d = day - op.day;
        if (m <= 2){ y--; m += 12;}
        return 365 * y + y / 4 - y / 100 + y / 400 + (153 * (m - 3) + 2) / 5 + d - 307;
    }
    MY_DATE run(LL k){//当前日期过k天
        k += p();
        LL x = k + 1789995, n = 4 * x / 146097, i, j, d;
        x -= (146097 * n + 3) / 4;
        i = (4000 * (x + 1)) / 1461001;
        x -= 1461 * i / 4 - 31;
        j = 80 * x / 2447;
        d = x - 2447 * j / 80;
        x = j / 11;
        return MY_DATE(100 * (n - 49) + i + x, j + 2 - 12 * x, d);
    }
};

拉格朗日插值

构造一个过给定 $n$ 个点 $(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$ 的多项式 $$ f(x)=\sum_{i=1}^{n}(y_i * \prod_{j\neq i}\frac{x-x_j}{x_i-x_j}) $$

曼哈顿距离与切比雪夫距离

设 $A(x_1,y_1)$，$B(x_2,y_2)$

曼哈顿距离 $d(A,B)=|x_1-x_2|+|y_1-y_2|=\max(|(x_1+y_1)-(x_2+y_2)|,|(x_1-y_1)-(x_2-y_2)|)$

这正是 $(x_1+y_1,x_1-y_1)$ 与 $(x_2+y_2,x_2-y_2)$ 的切比雪夫距离

因此将每个点 $(x,y)$ 转化为 $(x+y,x-y)$ ，原坐标系的曼哈顿距离就是新坐标系下的切比雪夫距离

同理，将每个点 $(x,y)$ 转化为 $(\frac{x+y}{2},\frac{x-y}{2})$，原坐标系的切比雪夫距离就是新坐标系下的曼哈顿距离

几个数论函数的数量级

约数个数 $d(n)$ 的数量级

A066150 Maximal number of divisors of any n-digit number

4, 12, 32, 64, 128, 240, 448, 768, 1344, 2304, 4032, 6720, 10752, 17280, 26880, 41472, 64512, 103680, 161280, 245760, 368640, 552960, 860160, 1290240, 1966080, 2764800, 4128768, 6193152, 8957952, 13271040, 19660800, 28311552, 41287680, 59719680, 88473600, 127401984, 181665792, 264241152, 382205952, 530841600

素数数量 $\pi(n)$ 的数量级

A006880 Number of primes < 10^n

0, 4, 25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534, 455052511, 4118054813, 37607912018, 346065536839, 3204941750802, 29844570422669, 279238341033925, 2623557157654233, 24739954287740860, 234057667276344607, 2220819602560918840, 21127269486018731928, 201467286689315906290

素数前缀和 $\sum_n[n\in prime]n$ 的数量级

A046731 a(n) = sum of primes < 10^n

0, 17, 1060, 76127, 5736396, 454396537, 37550402023, 3203324994356, 279209790387276, 24739512092254535, 2220822432581729238, 201467077743744681014, 18435588552550705911377, 1699246443377779418889494, 157589260710736940541561021, 14692398516908006398225702366

高斯整数

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

struct G{
    LL r,i;
    G(LL re=0,LL im=0):r(re),i(im){}
    LL norm()const{return r*r+i*i;}
    G operator+(const G&b)const{return G(r+b.r,i+b.i);}
    G operator-(const G&b)const{return G(r-b.r,i-b.i);}
    G operator*(const G&b)const{return G(r*b.r-i*b.i,i*b.r+r*b.i);}
    G operator/(const G&b)const{LL l=b.norm();return G(round(1.0*(r*b.r+i*b.i)/l),round(1.0*(i*b.r-r*b.i)/l));}
    G operator%(const G&b)const{return (*this)-(*this)/b*b;}
    //void get(){read(r,i);}void write()const{printf("%lld %lld ",r,i);}
};G gcd(G a,G b){return (b.norm()==0)?a:gcd(b,a%b);}

环染色问题

poj2154 Burnside引理

#include<cstdio>
using namespace std;
typedef long long LL;
#define pln putchar('\n')
const int N=100009;

int n,p,ans;

LL fp(LL a,LL b,LL Mod){
    LL res=(Mod!=1);
    for(;b;b>>=1,a=a*a%Mod)if(b&1)res=res*a%Mod;
    return res;
}

int ntp[N],pri[N],tot;

int f[33],t[33],d;

void dfs(int x,int phi,int u){
    if(x>d){
        ans=(ans+phi*fp(n,n/u-1,p))%p;
        return;
    }
    dfs(x+1,phi,u);
    u*=f[x];
    phi*=(f[x]-1);
    dfs(x+1,phi,u);
    for(int i=2;i<=t[x];++i){
        u*=f[x];
        phi*=f[x];
        dfs(x+1,phi,u);
    }
}

void solve(int T){
    scanf("%d%d",&n,&p);d=0;
    int x=n;
    for(int i=1;i<=tot;++i){
        if(pri[i]*pri[i]>x)break;
        if(x%pri[i])continue;
        int k=0;while(x%pri[i]==0)x/=pri[i],++k;
        f[++d]=pri[i];
        t[d]=k;
    }
    if(x>1){
        f[++d]=x;
        t[d]=1;
    }
    ans=0;
    dfs(1,1,1);
    printf("%d\n",ans);
}

signed main(){
    for(int i=2;i<=40000;++i){
        if(!ntp[i])pri[++tot]=i;
        for(int j=1;j<=tot;++j){
            int nex=pri[j]*i;
            if(nex>40000)break;
            ntp[nex]=1;
            if(i%pri[j]==0)break;
        }
    }
    int t;scanf("%d",&t);
    for(int i=1;i<=t;++i)solve(i);
    return 0;
}

置换群

Burnside引理

设有限置换群 $G$ 作用在有限集 $X$ 上，则 $X$ 上的 $G-轨道$ 数量为 $$ \frac{1}{|G|}\sum_{g\in G}|{x|g(x)=x}| $$

若将 $x\in X$ 理解为对 $1..n$ 每个元素映射到某种颜色的方案，$g$ 的不动点数量可以理解为：

对于 $g$ 生成的图（可分解为若干不相交简单环的并），满足同一个环上所有点颜色相同，的染色方案数

可旋转涂色问题

给一个 $n$ 元环涂色，$m$ 种颜色，旋转后得到的方案算同一种，求不同的方案数

$\frac{1}{n}\sum^{n-1}{i=0}m^{gcd(n,i)}=\frac{1}{n}\sum{d|n}\phi(d)m^{n/d}$

轮换指标

$x_i^j$ 表示置换中 $i$ 元环有 $j$ 个

正 $n$ 边形的旋转群 $$ \frac{1}{n}\sum_{d|n}\phi(d)x_d^{n/d} $$
正 $n$ 边形的二面体群 $$ \frac{1}{2n}\sum_{d|n}\phi(d)x_d^{n/d}+ \begin{cases}\frac{1}{2}x_1x_2^{\frac{n-1}{2}}\qquad & {n\ is\ odd} \ \frac{1}{4}(x_2^{\frac{n}{2}}+x_1^2x_2^{\frac{n-2}{2}})\qquad & {n\ is\ even} \end{cases} $$
正方体的顶点置换群 $$ \frac{1}{24}(x_1^8+8x_1^2x_3^2+9x_2^4+6x_4^2) $$
正方体的边置换群 $$ \frac{1}{24}(x_1^{12}+8x_3^4+6x_1^2x_2^5+3x_2^6+6x_4^3) $$
正方体的面置换群 $$ \frac{1}{24}(x_1^6+8x_3^2+6x_2^3+3x_1^2x_2^2+6x_1^2x_4) $$

多项式快速幂

给定 $F(x)(mod\ x^n)$ ，求 $F^k(x)(mod\ x^n)$

$k,n$ 不大，系数集不成域，$O(n \log n\log k)$ 分治NTT

$k$ 很大，且系数模奇质数 $p$ ($p>n$)：

$F[0]=1$，根据 $f(x^p)\equiv f^p(x)(mod\ p)$ ，则 $F^p(x)\equiv F[0]=1(mod\ x^p)=1(mod x^n)$，直接令 $k%=p$ 后计算 $\exp(k\ln F(x))$

$F[0]!=1$ ，提出最低次非 $0$ 系数 $ux^t$ 后化为上一种情形 $F^k(x)=u^kx^{tk}G^k(x)$，其中 $G[0]=1$，注意 $u$ 的幂模 $p-1$ ，$G$ 的幂模 $p$

广义二项式

对于任意实数 $\alpha$ 和非负整数 $k$ ： $$ C(\alpha,k)=\frac{\alpha^{k\over}}{k!} $$ $$ C(-\alpha,k)=(-1)^kC(\alpha+k-1,k) $$

有牛顿二项式： $$ (1+x)^\alpha=\sum_{k\geq 0}C(\alpha,k)x^k $$ $$ \frac{1}{(1-x)^\alpha}=\sum_{k\geq 0}C(\alpha+k-1,k)x^k $$

Fibonacci数列

$F_1=1,F_2=1,F_{n+2}=F_{n+1}+F_n$

$F_n=\frac{\sqrt 5}{5}((\frac{1+\sqrt 5}{2})^n-(\frac{1-\sqrt 5}{2})^n)$

Catalan数

$H_0=1,H_1=1$

$H_n=\sum_{i=0}^{n-1}H_iH_{n-i-1}$ $n\geq 2$

$H(x)=\sum_{n\geq 0}\frac{C(2n,n)}{n+1}x^n$

错排数

容斥原理 $n! \sum_ {k=0}^n\frac{(−1)^k}{k!}$

盒子放球方案数

luoguP5824 记得检查 n<=m 之类的合法性条件

给定正整数 $n,m$，问 $n$ 个球全部放入 $m$ 个盒子的方案数，不同的限制条件如下

划分数(与第二类Stirling数相关)

球之间互不相同，盒子之间互不相同。
- $m^n$
球之间互不相同，盒子之间互不相同，每个盒子至多装一个球。
- $m^{n\over}$
球之间互不相同，盒子之间互不相同，每个盒子至少装一个球。
- $m!S_2(n,m)$
球之间互不相同，盒子全部相同。
- $\sum_{i=1}^m S_2(n,i)$
球之间互不相同，盒子全部相同，每个盒子至多装一个球。
- $[n\leq m]$
球之间互不相同，盒子全部相同，每个盒子至少装一个球。
- $S_2(n,m)$

一类不定方程解数

球全部相同，盒子之间互不相同。
- $C(n+m-1,m-1)$
球全部相同，盒子之间互不相同，每个盒子至多装一个球。
- $C(m,n)$
球全部相同，盒子之间互不相同，每个盒子至少装一个球。
- $C(n-1,m-1)$

盒子有限的整数分拆

球全部相同，盒子全部相同。
- $x^n$
球全部相同，盒子全部相同，每个盒子至多装一个球。
- $[n\leq m]$
（整数的m分拆）球全部相同，盒子全部相同，每个盒子至少装一个球。
- OGF： $n<m\ ?\ 0\ :\ x^{n-m}$
- 记答案为 $p(n,m)$ ，递推式 $p(n,m) = p(n-1,m-1)+p(n-m,m)$

第二类Stirling数

$n$ 个带标号球全部放入 $m$ 个无标号盒子，且所有盒子非空的方案数

$S_2(n,0)=[n=0]$
$S_2(n,m)=S_2(n-1,m-1)+mS_2(n-1,m)$
$S_2(n,m)=\sum_{i=0}^{m}\frac{(-1)^i\ (m-i)^n}{i!\ (m-i)!}$ 可以卷积算一行
考虑 $n$ 个带标号球全部放入 $m$ 个带标号盒子共有 $m^n$ 种方案，有和式 $m^n=\sum_{k=0}^mS_2(n,k)m^{k\over}$
列EGF $\sum_{n\geq 0}S_2(n,m)\frac{x^n}{n!}=\frac{1}{m!}(e^x-1)^m$

分拆数

$n$ 个无标号球全部放入一些(个数无限制)无标号盒子，要求每个盒子非空

记方案数为 $P(n)$

$P(n)=\sum_{m=1}^n p(n,m)$
OGF：$\sum_{n\geq 0}P(n)x^n=\prod_{n\geq 1}\frac{1}{1-x^n}$

五边形数定理

广义五边形数 $g_n=\frac{1}{2}n(3n-1)$

对于 $n=0,1,-1,2,-2,3,-3,...$

$g_n$ 的前几项为 $0,1,2,5,7,12,15,22,26,...$

欧拉函数的OGF $$ \phi(x)=\prod_{i=1}^{\inf}(1-x^i) =1-x-x^2+x^5+x^7-x^{12}-x^{15} \ =1+\sum_{i=1}^{\inf}(-1)^i(x^{i(3i-1)/2}+x^{-i(-3i-1)/2}) $$

其中 $i(3i-1)/2$ 与 $-i(-3i-1)/2$ 恰好是相邻的广义五边形数

由 $P(x)\phi(x)=1$ 得

$P(1)=1$

$P(n)=P(n-1)+P(n-2)-P(n-5)-P(n-7)+...$

有时可以用来 $O(n\sqrt n)$ 预处理分拆数，比约 $O(50nlogn)$ 的龟速MTT省事(HDU6042)

第一类Stirling数

求有 $k$ 个轮换的 $n$ 元置换的方案数，记为 $S_1(n,k)$

$S_1(n,0)=[n=0]$
$S_1(n,k)=S_1(n-1,k-1)+(n-1)S_1(n-1,k)$
行OGF（无符号） $x^{\over n}=\sum_{k=0}^nS_1(n,k)x^k$
行OGF（符号） $x^{n \over }=\sum_{k=0}^nS_1(n,k)x^k$
列EGF（有符号） $\sum_{n\geq 0}S_1(n,k)\frac{x^n}{n!}=\frac{1}{k!}ln^k(1+x)$
列EGF（无符号）$\sum_{n\geq 0}S_1(n,k)\frac{x^n}{n!}=\frac{1}{k!}ln^k(1/(1-x))$

Stirling反演(未验证)

$\sum_{k\geq 0}S_2(n,k)S_1(k,m)=[n=m]$

$f_n=\sum_{i=0}^nS_2(n,i)g_i\iff g_n=\sum_{i=0}^nS_1(n,i)f^i$

Files

1-数学.md

Latest commit

History

1-数学.md

File metadata and controls

数学

模乘

快速乘

龟速乘

快速幂

光速幂

最大公约数（欧几里得算法）

binary-gcd

乘法逆元

exgcd求逆元

线性递推求逆元

线性同余方程组的合并（ex中国剩余定理）

大指数运算（欧拉降幂）

Miller-Rabin 素性检测

Pollard-Rho 分解质因数

原根

寻找质数原根

高次同余方程

离散对数（BSGS/exBSGS）

模任意数开任意次根/N次剩余

特殊情形一

特殊情形二

模奇质数开平方根/立方根

二次剩余判定定理

求平方根（Cipolla算法）

立方根（三次剩余）

线性筛

欧拉函数 $\phi(x)$

整除分块

莫比乌斯变换（数论）/狄利克雷卷积

莫比乌斯反演（数论）

杜教筛

特别注意 对 n 数论分块时，最后一次循环之前，执行 i=j+1 之后，此时 $i=n+1$，如果 n 刚好为数据类型的上界，会导致溢出！！！

Min_25筛

Step.1

Step.2

法一

法二

组合数的计算

基于模数的预处理分治（exlucas）

直线下整点数（类欧几里得算法）

洛谷模板题

高斯整数

四则运算与辗转相除

$4k+1$ 型质数的分解（二次剩余）

平分分解计数

数值积分

日期操作

用于跳转的常量

辅助函数

日期和整数的一一对应

拉格朗日插值

曼哈顿距离与切比雪夫距离

几个数论函数的数量级

约数个数 $d(n)$ 的数量级

素数数量 $\pi(n)$ 的数量级

素数前缀和 $\sum_n[n\in prime]n$ 的数量级

高斯整数

环染色问题

置换群

Burnside引理

可旋转涂色问题

轮换指标

多项式快速幂

广义二项式

Fibonacci数列

Catalan数

错排数

盒子放球方案数

划分数(与第二类Stirling数相关)

一类不定方程解数

盒子有限的整数分拆

第二类Stirling数

分拆数

特别注意对 `n` 数论分块时，最后一次循环之前，执行 `i=j+1` 之后，此时 $i=n+1$，如果 `n` 刚好为数据类型的上界，会导致溢出！！！