exBSGS算法

解决模数、底数不互质的离散对数问题。

（1）为何(BSGS)算法不再适用：(A)不一定存在逆元，而且无法保证解的循环性。
（2）无解的结论：
设方程为(A^x=B pmod{P})
当 ((A,P) mid B)且(B e 1) 时，原方程无自然数解。
还有就是(A=0,B≠0)这种。
（3）算法流程：
先判无解。
然后若(B=1)，显然(x=0)，特判之。
否则 (G=(A,P))
若(G>1)，
两边同除以(G),得(A^{x-1}=frac{B}{G}*(frac{A}{G})^{-1} pmod{frac{P}{G}})
迭代至((A,P)==1)时即可套用(BSGS)算法了。

int BSGS(int a,int b,int P)
{
    int s,o;
    S.clear();
    int m=(int)ceil(sqrt(P));
    o=a;
    s=b;
    for(int i=1; i<=m; i++)
    {
        s=(ll)s*o%P;
        S[s]=i;
    }
    o=fpow(a,m,P);
    s=1;
    for(int i=1; i<=m; i++)
    {
        s=(ll)s*o%P;
        if((it=S.find(s))!=S.end())
            return i*m-it->second;
    }
    return -1;
}

int exbsgs(int a,int b,int c)
{
    a%=c; b%=c;
    if(b==1||c==1)return 0;
    if(!a)return b?-1:1;
    int d=gcd(a,c);
    if(d==1)return BSGS(a,b,c);
    if(b%d!=0)return -1;
    int o=exbsgs(a,(ll)b/d*Inv(a/d,c/d)%(c/d),c/d);
    return o==-1?-1:o+1;
}

快速计算阶乘

即快速计算：
对于质数(p)，把(N!)去掉所有(p)因子的部分对(p^e)取模，(p^ele 10^5)
（1）为了便于统计出现了多少个p的次幂，先将阶乘中所有p的倍数提出来。可以简单算出：共有
$displaystyle lfloor frac{n}{p} floor $ 个。
这中间每一项都除去p，
可以得到 (displaystyle lfloor frac{n}{p} floor !) 。该部分可以选择递归求解。
（2）那么接下来只剩下非p的倍数的几项了。通过简单观察可以知道，剩余几项具有循环节，循环节长度小于(p^e) 。
原因是剩余项关于p具有循环节，而
(x+p^e equiv x pmod{p^e})，
所以可以一起计算。
结果会剩下几项凑不齐一个循环节，但是这几项长度已经小于(p^e)了，可以选择暴力乘法求解。
有必要举个例子模拟一下：

(N=22,p=3,e=2,)
(22!=1*2*3*......*22)
(=(1*2*4*5*7*8*10*11*13*14*16*17*19*20*22)*3^7*7!)
(7!)递归处理；
而
((1*2*4*5*7*8*10*11*13*14*16*17*19*20*22))
(equiv (1*2*4*5*7*8)^2*(19*20*22) pmod {3^2})
两个括号里的暴力计算；
考虑对 (p^e) 以内做预处理，查询复杂度可以做到 (O(log_p n)) 左右。
代码在跟下面一起贴

组合数取模（扩展卢卡斯）

即快速计算(dbinom{N}{M} pmod{P})，(N,Mle 10^9,P不一定是质数,Ple 10^9),但(P)中任何一个质因子的总积不超过(10^5)

考虑中国剩余定理：
令(Ans=dbinom{N}{M})

对质因子分开计算有：
(Ans equiv a_1 pmod{p_1^{q_1}})
(Ans equiv a_2 pmod{p_2^{q_2}})
(......)
对于一个式子：
不含(p_i)的部分使用上述快速阶乘方法求解。
含有(p_i)的重数使用(displaystyle sum_{i=1} lfloor frac{n}{p^i} floor)计算。
然后两部分分开算，乘起来就是对应的(a_i)。

最后CRT合并起来就可以了。

[Ansequiv sum a_i imes (prod_{j eq i}b_j^{-1} mod {b_i}) imes prod_{j eq i}b_j pmod {prod b_i} ]

struct bnm {
	int p, Mod, fac[1 << 20];
	void init(int x, int y)
	{
		p = x;
		Mod = y;
		fac[0] = 1;
		for(int s = 1; s != Mod; ++s)
			if(s % p)
				fac[s] = 1LL * fac[s - 1] * s % Mod;
			else
				fac[s] = fac[s - 1];
	}
	int Fac(ll x)
	{
		if(!x)
			return 1;
		return 1LL * Fac(x / p) * fexp(fac[Mod - 1], x / Mod, Mod) % Mod * fac[x % Mod] % Mod;
	}
	int operator()(ll n, ll m) {
		if(n < m || m < 0) return 0;
		int res = Fac(n) * Inv(Fac(m), Mod) % Mod * Inv(Fac(n - m), Mod) % Mod;
		int z = 0;
		for(ll d = n; d; z += d /= p);
		for(ll d = m; d; z -= d /= p);
		for(ll d = n - m; d; z -= d /= p);
		res = 1LL * res * fexp(p, z, Mod) % Mod;
		return res;
	}
}Solve;

普通二次剩余

即解方程(x^2=a pmod{P})，P是质数。

（1）勒让德符号的定义
((frac{a}{p})=1)：a在模p意义下有二次剩余。
((frac{a}{p})=-1)：a在模p意义下无二次剩余。
计算公式：((frac{a}{p}) = a^{frac{p-1}{2}} pmod p)
（2）求解二次剩余
随机一个(a)使得(b=sqrt{a^2-n})（此时 (a^2-n) 无二次剩余，就直接写成根号形式），然后把(b)作为二次域的单位元。
那么(x=(a+b)^{frac{p+1}{2}})（证明需要二项式定理展开，具体参考这里 !ACdreamer）
写一个扩域乘法类就好了。
注意这里(p)得是奇质数，但(p=2)也是很容易特判的。
二次剩余有一正一负两解，下面代码随便取了一个。

namespace Cipolla
{
	inline int Le(int x){return fpow(x,(P-1)/2);}
	
	int w;
	struct Cp
	{
		int x,y;
		Cp(int a=0,int b=0){x=a;y=b;}
	};
	Cp operator+(Cp a,Cp b){return Cp((a.x+b.x)%P,(a.y+b.y)%P);}
	Cp operator-(Cp a,Cp b){return Cp((a.x-b.x+P)%P,(a.y-b.y+P)%P);}
	Cp operator*(Cp a,Cp b){return Cp(((ll)a.x*b.x+(ll)w*a.y%P*b.y)%P,((ll)a.x*b.y+(ll)a.y*b.x)%P);}
	Cp operator^(Cp a,int b){
		Cp o(1,0);
		for(;b;b>>=1) {
			if(b&1)o=o*a;
			a=a*a;
		}
		return o;
	}
	
	int calc(int x)
	{
		x%=P;
		int a;
		while(1) {
			a=rand();
			w=((ll)a*a-x+P)%P;
			if(Le(w)==P-1) break;
		}
		Cp s=Cp(a,1)^((P+1)/2);
		return min(s.x,P-s.x);
	}
}

原根和阶

奇素数与奇素数的方幂有原根，原根 (g) 是阶恰好等于 (p-1) 的数，也就是 (g^x(xin [0,p-2])) 与整数 (yin [1,p-1]) 一一对应。
我们一般用暴力枚举的方法求原根，原理类似试除法，判断是否存在 (d|p-1) 使得 (g^dequiv 1 pmod p) 即可。

bool judge(int x)
{
	for(int j = 0; j < tot; j ++) //只需对质因子做
		if(fexp(x, (P - 1) / c[j], P) == 1)
			return false;
	return true;
}

对一个任意模数求阶也不是很难，只要对 (varphi (m)) 不断试除保证 (x^kequiv 1 pmod m)即可。

原根可以转换为单位根，即若 (n|(varphi(p)-1)) 则 (omega_n=g^{frac{varphi(p)-1}{n}})。
常用于单位根反演 ([n|i]=dfrac{1}{n}sum_{k=0}^{n-1} omega_n^{ki})
单位根反演可以应用于 ( ext{NTT}) 以及结合“二项式定理”优化复杂度。

快速乘法取模

某些情况下机器不支持int128，或者毒瘤出题人卡int128常数，就需要写这种 (O(1)) 快速乘。

inline ll mul(ll a,ll b,ll p){
	ll res=a*b-((ll)((long double)a*b/p+0.5))*p;
	return res<0?res+p:res;
}

一定要保证 (a,blt p) 的时候计算，因为它是由这个保证精度的。

合并同余方程组

(x equiv now pmod m)
(x equiv a pmod b)

由第一式得
(x=now+k*m)
代入二式
(now+km equiv a pmod b)
(mk equiv a-now pmod b)
拿个扩欧就能解出来最小的(k)
然后(now=now+k imes m)

有一个坑点就是你的(k)不要对(m)取模，
归纳易证你这个(now)一定是最小的正整数解，只要题目保证了啥啥，(now)就肯定不会溢出。

但是好多题目中(m)都会溢出，变成负数，再取模就会出bug。

当然还是不能排除求出假答案的情况，最后可以 check 一遍。

void excrt(ll a,ll b)
{
    ll c=((a-now)%b+b)%b;
    ll x,y;
    ll d=exgcd(M,b,x,y);
    if(c%d!=0) err();
    c/=d;
    b/=d;
    x=(x%b+b)%b;
    x=(__int128)x*c%b;
    now+=x*M;
    M*=b;
}

Miller-Rabin素数测试

二次探测原理：
当p是质数时，则一定有

若(x^2 equiv 1 pmod p)，则(x=1 pmod p)或(x=-1 pmod p)。

这样我们就考虑反过来想办法验证(n)是质数。

令(n=2^k imes m + 1)（m是奇数），

拿个小质数来进行二次探测，进行(k)次迭代即可。
取10-12个质数就足够稳了。

代码21行是二次探测，25行是费马小定理。

int ispr(ll n)
{
    if(n<2)
        return 0;
    for(int i=0; i<12; i++)
        if(n%Pr[i]==0)
            return n==Pr[i];
    ll m=n-1,x,y;
    int k=0;
    while(~m&1)
    {
        m>>=1;
        k++;
    }
    for(int i=0; i<12; i++)
    {
        x=fpow(Pr[i],m,n);
        for(int j=0; j<k; j++)
        {
            y=Mul(x,x,n);
            if(y==1&&x!=1&&x!=n-1)
                return 0;
            x=y;
        }
        if(x!=1)
            return 0;
    }
    return 1;
}

质因数分解 Pollard-rho

先说一下为大整数 (N(Nle 10^{18})) 分解质因子的流程：

（1）弄出来一个 (N) 的非平凡因子 (x)。
（2）递归分解 (x) 和 (N/x)。

而Pollard-rho就是高效地寻找非平凡因子的算法。

Pollard-rho基于随机化实现，然后算法的效率保证是“生日悖论”，这里的含义如下：

如果随机一个数 (x)，那么 (x) 是 (N) 约数的概率极其微小。
但是我们如果随机出 (x,y)，那么 (|x-y|) 是 (N) 约数的概率将显著提高。

（1）不断随机数字(x)，假如(gcd(x,n)>1)直接返回这个(gcd)

额样例都跑不出来，没有任何前途。

（2）按照上面讲的这么写，并按照倍增的形式不断调整：

x=y=0;
for(int k=2; ; k<<=1)
{
    for(int i=1; i<=k; i++)
    {
        x=(x*x+sd)%n;
        v=gcd(abs(x-y),n);
        if(v!=1 && v!=n)
            return v;
        if(x==y)
            return n;
    }
    y=x;
}

效率玄学地大幅提升，但是依然TLE飞起。。。
（3）然而你发现这个gcd常数巨大，我们把改成用一个模 (N) 意义下的变量记一下，这样每隔一会查一次就行！

static int Size(1<<7);
x=y=0;
for(int k=2; ;k<<=1)
{
    v=1;
    y=x;
    for(int i=1; i<=k; i++)
    {
        x=((u128)x*x+sd)%n;
        v=((u128)v*_abs(x-y))%n;
        if(i%Size==0)
        {
            z=gcd(v,n);
            if(z>1) return z;
        }
    }
    z=gcd(v,n);
    if(z>1) return z;
}

就能过了，常数极小。
Pollard-rho时间复杂度是(O(n^{0.25}))
实测每次找质因子内层乘法能稳定在(10^5)之内。
其实感觉特别暴力。