NumberTheory

快速幂

• 模数在 $int$ 范围

  ll qPow(ll a, ll b, ll c) {  //求(a^b) % c
      ll ret = 1;
      while (b) {
          if (b & 0x1) ret = ret * a % c;
          a = a * a % c;
          b >>= 1;
      }
      return ret;
  }

• 模数超出 $int$ 范围

  ll qMulti(ll x,ll y,ll n){ //求(x*y)%n
  	ll ret=0,tmp=x%n; 
  	while(y) { 
          if(y&0x1)if((ret+=tmp)>n)ret-=n;
          if((tmp<<=1)>n)tmp-=n; 
          y>>=1;
      }
      return ret; 
  }
  //O(1)快速乘
  ll qMulti(ll x,ll y,ll mod){
      return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;     
  }
  ll qPow(ll a,ll b,ll c){ //求(a^b) % c ，c超出int范围
      ll ret=1;
  	while(b) {
  		if(b&0x1)ret=qMulti(ret,a,c);
  		a=qMulti(a,a,c);
  		b>>=1;
  	}
  	return ret;
  }

• $\mathcal{O}(1)$ 快速幂

  int step;
  ll pre1[N];//pre1[i]表示a^{i*step}
  ll pre2[N];//pre2[i]表示a^{i}
  ll qPow(ll a,ll b,ll c){ //求(a^b) % c
  	return pre1[c/step]*pre2[c%step]%c;
  }

欧几里得算法

gcd

ll gcd(ll x,ll y){
    return y?gcd(y, x%y):abs(x);
}

exgcd

ll exgcd(ll a, ll b, ll &x, ll &y){ //求解ax+by=gcd(a,b),返回gcd(a,b)
    if(b==0){
    	x = 1;
		y = 0;
	    return a;
    }
	ll r = exgcd(b, a%b, y, x);
	y -= a/b*x; //这里已经是递归，回溯的过程了，x,y已经颠倒了
	return r;
}
/*
x = c/gcd(a,b)*x0 + b/gcd(a,b)*t;
y = c/gcd(a,b)*y0 - a/gcd(a,b)*t; 
(x0, y0 为方程的一组特解, t为整数)
*/

裴蜀定理

$ax+by=c$ 有解的充要条件为 $gcd(a,b) | c$

推广: $a_1x_1+a_2x_2+\cdots+a_nx_n=c$ 有解的充要条件为 $gcd(a_1,a_2,\cdots,a_n) | c$

类欧问题

素数和因子

打表

bool notP[N];
int prime[N],num_prime=0;
void get_prime(){
  	notP[1]=1;
    for(int i=2;i<N;i++){
        if(!notP[i]) prime[num_prime++]=i;
        for(int j=0;j<num_prime&&(ll)i*prime[j]<N;j++){
            int k = i*prime[j];
            notP[k] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}

区间素数

memset(visab,0,sizeof(visab));
for (int i = 0; i < cnt && prime[i] <= b; i++) {
	LL k = a / prime[i];
  	if (k * prime[i] < a) k++;
  	for (LL j = k * prime[i]; j <= b; j += prime[i]) {
  		visab[j - a] = 1;
 	}
}

唯一分解定理

pair<ll,int> factor[50];
int num_factor = 0;
void decomposition(ll n) {
    num_factor = 0;
    for (int i = 0; i < num_prime && (ll)prime[i] * prime[i] <= n && n > 1;++i) {
        if (n % prime[i] == 0) {
            int cnt=0;
            while (n % prime[i] == 0 ) {
                n /= prime[i];
                cnt++;
            }
            factor[num_factor++] = pair<ll,int>(prime[i],cnt);
        }
    }
    if (n > 1) factor[num_factor++] = pair<ll,int>(n,1);
}

pair<ll,int> factor[50];
int num_factor = 0;
void decomposition(ll n) {
    num_factor = 0;
    int m=(int)sqrt(n+0.5);
    for(int i=2;i<=m;i++){
        if(n%i==0){
            factor[num_factor]=pair<ll,int>(i,0);
            while(n%i==0){
                factor[num_factor].se++;
                n/=i;
            }
            num_factor++;
        }
    }
    if (n > 1) factor[num_factor++] = pair<ll,int>(n,1);
}

n 的因子数

分解 n，$\large n=p_1^{a_1}*p_2^{a_2}\cdots p_m^{a_m}$，n 的因子和为

$$\Large\tau(n)=\prod\limits_{i=1}^ma_i+1$$

n 的因子和

分解 n，$\large n=p_1^{a_1}*p_2^{a_2}\cdots p_m^{a_m}$，n 的因子和为

$$\Large\sigma(n)=\prod\limits_{i=1}^m(p_i^0+p_i^1+\cdots+p_i^{a_i})=\prod\limits_{i=1}^m\frac {p_i^{a_i+1}-1} {p_i-1}$$

1-n的因子数

$\Large\sum\limits_{i=1}^n\lfloor\frac n i\rfloor$

n=12

$\large\lfloor\frac n i\rfloor$	[l,r]
12	[1,1]
6	[2,2]
4	[3,3]
3	[4,4]
2	[5,6]
1	[7,12]

每个 $l$ 等于上一个 $r+1$，$r=n/(n/l)$

//solution1
ll res = 0;
for (int l = 1; l <= n; l = r + 1) {
    int r = n / (n / l);
    res += (n / l) * (r - l + 1);
}
//solution2
int k=sqrt(n+0.1);
for (int i = 1; i <= k; i++) {
	ans+=n/i;
	if (n/(i+1)<i)break;
	else ans+=(n/i-n/(i+1))*i;
}
//solution3
for (int i = 1; i <= k; i++)  ans+=n/i;
ans*=2;
ans-=k*k;

1-n的因子之和

统计 1-n 每个数取的次数$\Large\frac n i$，即$\Large\sum\limits_{i=1}^n\lfloor\frac n i\rfloor*i$

ll res = 0;
for (int l = 1; l <= n; l = r + 1) {
    int r = n / (n / l);
    res += (n / l) * (r + l) * (r - l + 1) / 2;
}

费马小定理

$$x^{p-1}\equiv1(mod\ p)$$

威尔逊定理

$$(p-1)!\equiv-1(mod\ p)$$

RSA

两个大素数密钥 $p,q$ ，令 $n=p\times q$ ，随机取加密密钥 $e$ ，使得 $gcd(e,(p-1)(q-1))=1$ ，使用 $exgcd$ 计算解密密钥 $d=e^{-1}\mod (p-1)(q-1)$ ， $d$ 与 $n$ 也是互质的。 $e$ 和 $n$ 是公开密钥， $d$ 是私人密钥。

设 $A$ 为明文， $B$ 为密文

$$B=A^e \mod n\\ A=B^d \mod n$$

积性函数

1. $f[1]=1$
2. $设n=\prod p_i^{k_i},则f(n)=f(\prod p_i^{k_i})=\prod f(p_i^{l_i})$
3. $设n=\prod p_i^{k_i},则f(n)=f(\prod p_i^{k_i})=\prod f(p_i)^{k_i}$ 若满足这一项称为完全积性函数。

欧拉函数

$${\phi(n)=\left\{\begin{array}{c}1&n=1\\\prod(p_i-1)p_i^{k_i-1}=n\prod\left(1-\frac1 p_i\right)& n=\prod p_i^{k_i}\\\end{array}\right.}$$

• $1\sim n$中与 n 互质的个数(包括 1 )

• 积性

  $$\phi(mn)=\phi(m)\phi(n)\ , \ gcd(m,n)=1$$

• $prime[i]\le x<prime[i+1]$，必定有$\phi(x)\le prime[i]$

• 欧拉定理

  $$a^{\phi(n)}\equiv1(mod\ n)\ ,\ gcd(n,a)=1$$

• 扩展欧拉定理

  $$\large A^K \equiv \left\{\begin{array}{c} A^{K\mod\phi(n)}&gcd(n,A)=1 \\ A^K&K< \phi(n)\\ A^{K\mod\phi(n)+\phi(n) }&K\ge \phi(n) \end{array}\right.\qquad(mod\ n)$$

• $1\sim n$中与 n 互质的数的和为$\large \frac {n*\phi(n)} 2+\epsilon(n)$

  设 gcd(n,p)=1，由更相减损术得 gcd(n-p,p)=1，且一定有p≠n-p，故 1~n 中与 n 互质的数可以由此两两配对，而每一对的和都为 n，所以总和为 n*φ(n)/2

• 求欧拉函数

  int euler(int n){
      int ans=1;
      for (int i=0; i<prime.size()&&prime[i]*prime[i]<=n&&n!=1; i++) {
          int e=0;
          while (n%prime[i]==0&&n!=1) {
              e++;
              n/=prime[i];
          }
          if (e) ans*=pow(prime[i], e-1)*(prime[i]-1);
      }
      if (n>1) ans*=(n-1);
      return ans;
  }
  O(sqrt(n))
  int euler(int n){
      int m=(int)sqrt(n+0.5);
      int ans=n;
      for(int i=2;i<=m;i++){
          if(n%i==0){
              ans=ans/i*(i-1);
              while(n%i==0) n/=i;
          }
      }
      if(n>1) ans=ans/n*(n-1);
      return ans;
  }

• 欧拉函数打表

  • 直接打表

    for(i=1; i<=maxn; i++) p[i]=i;
    for(i=2; i<=maxn; i+=1)
        if(p[i]==i){
            for(j=i; j<=maxn; j+=i)
                p[j]=p[j]/i*(i-1);
        }

  • 利用线筛打表O(n)

    bool notPrime[N+1];
    int prime[N+1],num_prime=0,phi[N];
    void get_prime(){
      	notPrime[1]=1;
        for(int i=2;i<=N;i++){
            if(!notPrime[i]) {
    			prime[num_prime++]=i;
              	phi[i]=i-1;
            }
            for(int j=0;j<num_prime&&(ll)i*prime[j]<=N;j++){
                notPrime[i*prime[j]] = 1;
                if(i%prime[j] == 0) {
                  	phi[i*prime[j]]=phi[i]*prime[j];
                  	break;
                }else phi[i*prime[j]]=phi[i]*(prime[j]-1);
            }
        }
    }

莫比乌斯函数

$$\mu(n)= \left\{ \begin{array}{c} 1&n=1\\ (-1)^r& n=\prod\limits_{i=1}^r p_i\\ 0&other \end{array} \right.$$

• 打表

  bool notPrime[N+1];
  int prime[N+1],num_prime=0,miu[N];
  void get_prime(){
    	notPrime[1]=1;
      for(int i=2;i<=N;i++){
          if(!notPrime[i]) {
  			prime[num_prime++]=i;
            	miu[i]=-1;
          }
          for(int j=0;j<num_prime&&(ll)i*prime[j]<=N;j++){
              notPrime[i*prime[j]] = 1;
              if(i%prime[j] == 0) {
                	miu[i*prime[j]]=0;
                	break;
              }else miu[i*prime[j]]=miu[i]*(-1);
          }
      }
  }

除数函数

$$\large\sigma_x(n)=\sum\limits_{d|n}d^x$$

一些完全积性函数

• 常函数

  $$\forall n\in N,1(n)=1$$

• 单位函数

  $$\forall n\in N,id(n)=n$$

• 幂函数

  $$\forall n\in N,id_k(n)=n^k$$

• 狄利克雷卷积单位函数

  $$\epsilon(n)= \left\{ \begin{array}{c} 1&n=1\\ 0&n>1 \end{array} \right.$$

卷积

乘积函数

$$(f*g)(n)=f(n)g(n)$$

Dirichlet卷积函数

$$(f\times g)(n)=\sum\limits_{d|n}f(d)g(\frac n d)$$

一些卷积

$$\epsilon \times 1=1\\ \mu\times1=\epsilon\\ \ \\ \phi\times1=id\\ \mu\times id=\phi$$

莫比乌斯反演

$$f\times 1=g\Longleftrightarrow \mu\times g=f$$

取模

除法取模

• b和m不互质

  $\frac a b \%m =\frac{a \% (bm)}b$

  证:

  设$\frac a b \%m=r$

  $\frac a b =km+r$

  $a=kbm+br$

  $br=a-kbm$

  $r=\frac{a \% (b*m)}b$

• 逆元( inv )

  $\frac a b \%m$时，因 $b$ 可能会过大，会出现爆精度的情况，所以需变除法为乘法，设 $c$ 是 $b$ 的逆元，则有 $b*c≡1(mod m)$，即 $\frac a b$ 的模等于 $a\times inv(b)$ 的模；

  • $m$ 为素数，根据费马小定理 $c=b^{m-2}(\mod m)$ ;

  • $b*c≡1(\mod m)$ 使用扩展欧几里得算法

  • 逆元打表

    //求1,2,...,N关于mod的逆元（mod为质数）
    ll inv[N];
    inv[1] = 1;
    for (int i = 2; i < N; i++) inv[i] = inv[mod % i] * (mod - mod / i) % mod;
    /*证明 : i* inv[i] = inv[mod % i] * ((i - 1) * mod + mod - mod / i * i) % mod;
    i* inv[i] = inv[mod % i] * ((i - 1) * mod + mod % i) % mod;
    i* inv[i] = (inv[mod % i] * (i - 1) * mod + 1) % mod;
    i* inv[i] = 1 % mod;*/
    
    //求N!关于mod的逆元（mod为质数）
    ll inv[N];
    inv[N-1] = qPow(fac[N-1], mod - 2, mod);
    for (int i = N - 2; i >= 0; i--) inv[i] = (inv[i + 1] * (i + 1)) % mod;
    /*证明 : inv[i + 1] * (i + 1) != 1 % mod;
    inv[i + 1] * (i + 1) * i != 1 % mod;*/

组合数取模

• n,m 1e5 预处理出逆元和阶乘

  ll fac[N],inv[N];
  ll C(ll a,ll b){
      if(b>a)return 0;
      return fac[a]*inv[b]%mod*inv[a-b]%mod;
  }
  void init(){//快速计算阶乘的逆元
      fac[0]=fac[1]=1;
      for(int i=2;i<N;i++)fac[i]=fac[i-1]*i%mod;
      inv[N-1] = qPow(fac[N-1], mod - 2, mod);
    	for (int i = N - 2; i >= 0; i--) inv[i] = (inv[i + 1] * (i + 1)) % mod;
  }

• n 1e9 , m 1e5

  for(int i=1;i<=m;i++)
      ans=ans*(n-i+1)%mod*qPow(i,mod-2,mod)%mod;

• n,m 1e18 , mod 1e5,素数 卢卡斯定理(Lucas)

  $$\left(\begin{array}{c}{m} \\ {n}\end{array}\right) \equiv \prod_{i=0}^{k}\left(\begin{array}{c}{m_{i}} \\ {n_{i}}\end{array}\right) \quad(\bmod p)\\ m=m_{k} p^{k}+m_{k-1} p^{k-1}+\cdots+m_{1} p+m_{0}\\ n=n_{k} p^{k}+n_{k-1} p^{k-1}+\cdots+n_{1} p+n_{0}$$

  ll Lucas(ll n,  ll m){
      if(m==0) return C(n, m);
      return C(n % mod, m % mod) * Lucas(n / mod, m / mod)%mod;
  }
  /*表达式：C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p。（可以递归）
  递归方程：(C(n%p, m%p)*Lucas(n/p, m/p))%p。（递归出口为m==0，return 1）*/

• exLucas

  求解 $C_n^m \mod p$ ， $p$ 不是素数

  • 分解 $\large p=\prod p_i^{k_i}$ ，分别求 $C_n^m\mod p_i^{k_i}$ ，使用中国剩余定理合并

  • 求解 $C_n^m\mod p_i^{k_i}$

    $$C_n^m=\frac {n!}{m!\times (n-m)!}\\ \text{提出 }n!,m!,(n-m)!\text{ 中的素因子 }p\\ \frac {\frac {n!}{p^{k1}}}{\frac {m!}{p^{k2}}\times\frac {(n-m)!}{p^{k3}}}\times p^{k1-k2-k3}$$

    易证 $k1-k2-k3\ge0$

  • 求解 $f(n)=\large\frac {n!}{p^a}\mod p^k$ ，分三部分

    • $1\sim n$ 中拥有 $p$ 因子的数的贡献：$f(\lfloor\frac n p\rfloor)$
    • $\Large\prod\limits_{i=1,(i,p)=1}^{p^k}i^{\lfloor\frac n {p^k}\rfloor}$
    • $\Large\prod\limits_{i=1,(i,p)=1}^{n\mod p^k}i$

    ll exgcd(ll a, ll b, ll &x, ll &y){ //求解ax+by=gcd(a,b),返回gcd(a,b)
        if(b==0){
        	x = 1;
    		y = 0;
    	    return a;
        }
    	ll r = exgcd(b, a%b, y, x);
    	y -= a/b*x; //这里已经是递归，回溯的过程了，x,y已经颠倒了
    	return r;
    }
    ll Chinese_Remainder(ll a[],ll m[],int len){//中国剩余定理  a存放余数  m存放两两互质的数
        ll M=1,ret=0;
        for(int i = 0; i < len; i++)  M *= m[i];
        for(int i = 0; i < len; i++){
            ll Mi = M/m[i],x,y;
            exgcd(m[i],Mi,x,y);
            ret = (ret+y*Mi*a[i])%M;
        }
        return (M+ret%M)%M;
    }
    ll qPow(ll a, ll b, ll c) {  //求(a^b) % c
        ll ret = 1;
        while (b) {
            if (b & 0x1) ret = ret * a % c;
            a = a * a % c;
            b >>= 1;
        }
        return ret;
    }
    ll pre[N],cnt;
    ll fac(ll n,ll pk,ll p){
        if(n<p)return pre[n];
        cnt+=n/pk;
        return pre[n%pk]*fac(n/p,pk,p)%pk;
    }
    ll vp(ll n,ll p){
        ll ans=0;
        while(n){
            ans+=n/p;
            n/=p;
        }
        return ans;
    }
    ll C(ll n,ll m,ll pk,ll p){
        pre[0]=1;
        for(int i=1;i<=pk;i++){
            pre[i]=pre[i-1];
            if(i%p)pre[i]=pre[i-1]*i%pk;
        }
        cnt=0;
        ll ans=fac(n,pk,p)*qPow(pre[pk],cnt,pk)%pk*qPow(p,vp(n,p)-vp(m,p)-vp(n-m,p),pk)%pk;
        ll x,y;
        cnt=0;
        exgcd(fac(m,pk,p)*qPow(pre[pk],cnt,pk)%pk,pk,x,y);
        ans=ans*x%pk;
        cnt=0;
        exgcd(fac(n-m,pk,p)*qPow(pre[pk],cnt,pk)%pk,pk,x,y);
        ans=ans*x%pk;
        return ans;
    }
    ll exLucas(ll n,ll m,ll p){
        ll a[100],b[100];
        int sq=sqrt(p+0.1);
        int tot=0;
        for(int i=2;i<=sq&&p>1;i++){
            int pk=1;
            while(p%i==0){
                pk*=i;
                p/=i;
            }
            if(pk>1){
                a[tot]=C(n,m,pk,i);
                b[tot++]=pk;
            }
        }
        if(p>1){
            a[tot]=C(n,m,p,p);
            b[tot++]=p;
        }
        return Chinese_Remainder(a,b,tot);
    }

• $n!$ 中 $p$ (素数) 的幂次

  $$v_p(n!)=\sum_{i=1}^{\infty}\lfloor\frac n {p^i}\rfloor$$

• 库默尔定理

  设 $m,n(\ge0)$，$p$ 为素数，则 $v_p(C_{m+n}^m)$ 等于在 $p$ 进制下 $m+n$ 时进位次数。

连分数

$$\Huge a_0+\frac 1 {a_1+\frac 1{a_2+\frac 1 {\cdots+\frac 1 {a_n}}}}$$

记为

$$[a_0,a_1,\cdots,a_n]=\frac {p_n}{q_n}\\ p_0=a_0,p_1=a_1a_0+1,p_k=a_kp_{k-1}+p_{k-2}\ (2\le k\le n)\\ q_0=1,q_1=a_1,q_k=a_kq_{k-1}+q_{k-2}\ (2\le k\le n)\\ p_nq_{n-1}-p_{n-1}q_n=(-1)^{n-1}\ (n\ge1)$$

佩尔方程

$x^2-dy^2=1$

• 定理1: 对于正整数 p,q ,如果有 $$|p^2-a^2q^2|<a$$ 则 $\frac p q$ 必为 a 的一个渐进值 $$x_{n+1}=x_0x_n+d*y_0y_n\\ y_{n+1}=y_0x_n+x_0y_n\\ .\\ \begin{pmatrix} x_n\\y_n \end{pmatrix}=\begin{pmatrix} x_0&d*y_0\\y_0&x_0 \end{pmatrix}^n\begin{pmatrix} x_0\\y_0 \end{pmatrix}$$

//1
//求出p2 - D * q2 = 1的基解(最小正整数解),这个可能溢出,有必要的话,用java, 写成类比较好
bool PQA(LLI D, LLI &p, LLI &q) {//来自于PQA算法的一个特例
    LLI d = sqrt(D);
    if ((d + 1) * (d + 1) == D) return false;
    if (d * d == D)             return false;
    if ((d - 1) * (d - 1) == D) return false;//这里是判断佩尔方程有没有解
    LLI u = 0, v = 1, a = int(sqrt(D)), a0 = a, lastp = 1, lastq = 0;
    p = a, q = 1;
    do {
        u = a * v - u;
        v = (D - u * u) / v;
        a = (a0 + u) / v;
        LLI thisp = p, thisq = q;
        p = a * p + lastp;
        q = a * q + lastq;
        lastp = thisp;
        lastq = thisq;
    } while ((v != 1 && a <= a0));//这里一定要用do~while循环
    p = lastp;
    q = lastq;
    //这样求出后的(p,q)是p2 – D * q2 = (-1)k的解，也就是说p2 – D * q2可能等于1也可能等于-1,如果等于1,(p,q)就是解,如果等于-1还要通过(p2 + D * q2,2 * p * q)来求解，如下
    if (p * p - D * q * q == -1) {
        p = lastp * lastp + D * lastq * lastq;
        q = 2 * lastp * lastq;
    }
    return true;
}

//2
typedef long long ll;  
ll a[20000];  
bool pell_minimum_solution(ll n,ll &x0,ll &y0){  
    ll m=(ll)sqrt((double)n);  
    double sq=sqrt(n);  
    int i=0;  
    if(m*m==n)return false;//当n是完全平方数则佩尔方程无解  
    a[i++]=m;  
    ll b=m,c=1;  
    double tmp;  
    do{  
        c=(n-b*b)/c;  
        tmp=(sq+b)/c;  
        a[i++]=(ll)(floor(tmp));  
        b=a[i-1]*c-b;  
        //printf("%lld %lld %lld\n",a[i-1],b,c);  
    }while(a[i-1]!=2*a[0]);  
    ll p=1,q=0;  
    for(int j=i-2;j>=0;j--){  
        ll t=p;  
        p=q+p*a[j];  
        q=t;  
        //printf("a[%d]=%lld %lld %lld\n",j,a[j],p,q);  
    } 
    if((i-1)%2==0){x0=p;y0=q;} 
    else{x0=2*p*p+1;y0=2*p*q;} 
    return true;  
}  

int main(){  
    ll n,x,y;  
    while(~scanf("%lld",&n)){  
        if(pell_minimum_solution(n,x,y)){  
            printf("%lld^2-%lld*%lld^2=1\t",x,n,y);  
            printf("%lld-%lld=1\n",x*x,n*y*y);  
        }  
    }

中国剩余定理

• 互质

  $$(S):\left\{\begin{array}{} x\equiv a_1 \ (mod\ m_1)\\ x\equiv a_2 \ (mod\ m_2)\\\quad\vdots\\x\equiv a_n \ (mod\ m_n)\\ \end{array}\right.$$

  假设整数$m_1,m_2,\cdots,m_n$两两互质，则对任意的整数$a_1,a_2,\cdots,a_n$ 方程组$(S)$有解

  设

  $$M=m_1\times m_2\times \cdots \times m_n=\prod\limits_{i=1}^nm_i\\ M_i=\frac M {m_i}\\ t_i=M_i^{-1}\\t_iM_i\equiv1\ (mod\ m_i)$$

  通解形式为

  $$x=a_1t_1M_1+a_2t_2M_2+\cdots+a_nt_nM_n+kM=kM+\sum\limits_{i=1}^na_it_iM_i\ ,k\in\mathbb{Z}$$

  ll Chinese_Remainder(ll a[],ll m[],int len){//中国剩余定理  a存放余数  m存放两两互质的数
      ll M=1,ret=0;
      for(int i = 0; i < len; i++)  M *= m[i];
      for(int i = 0; i < len; i++){
          ll Mi = M/m[i],x,y;
          exgcd(m[i],Mi,x,y);
          ret = (ret+y*Mi*a[i])%M;
      }
      return (M+ret%M)%M;
  }

• 非互质

  /*采用两两合并的思想
  x=a1+m1*x1 
  x=a2+m2*x2
  得到 
  a1+m1*x1 = a2+m2*x2 → m1*x1+m2*x2 = a2-a1
  再通过扩展欧几里得算法解出x1的最小正整数解，代入 
  x=a1+m1*x1
  得到x后合并为一个方程的结果为 
  y ≡ x(mod lcm(m1,m2))
  */
  pair<ll, ll> linear(ll B[], ll M[], ll n){//求解 x = B[i] (mod M[i]),总共n个线性方程组
      ll x = 0, m = 1;
      for(int i = 0; i < n; i ++) {
          ll a =  m, b = B[i] - x, d = gcd(M[i], a);
          if(b % d != 0)  return pair<ll, ll>(-1, -1);//答案不存在，返回-1
          ll q ,w;
          exgcd(a/d, M[i]/d, q, w);
          ll t = b/d * q%(M[i]/d);
          x = x + m*t;
          m *= M[i]/d;
      }
      x = (x % m + m ) % m;
      return pair<ll, ll>(x, m);//返回的x就是答案，m是最后的lcm值
  }

整数拆分

1. 将n划分成若干正整数(可相同)之和的划分数。

   • dp

     $dp[n][m]= dp[n][m-1]+ dp[n-m][m]$, $dp[n][m]$表示整数 n 的划分中，每个数不大于 m 的划分数。 划分分为两种情况:

     a.划分中每个数都小于 m，相当于每个数不大于 m - 1, 故划分数为$dp[n][m-1]$

     b.划分中有一个数为 m. 那就在 n中减去 m ,剩下的就相当于把 n-m 进行划分， 故划分数为 $dp[n-m][m]$;

     $$f(n,m)=\left\{\begin{array}{} 1&(n=1\parallel m=1)\\ f(n,n)&(n<m)\\ 1+f(n,m-1)&(n=m)\\f(n-m,m)+f(n,m-1)&(n>m)\end{array}\right.$$

   • dp

     $dp[n][k]= dp[n-k][k]+ dp[n-1][k-1]$ ,$dp[n][k]$表示将n划分成k个正整数之和的划分数，划分分为两种情况:

     a.不含 1 的分法，为保证每份都 >= 2，可以先拿出 k 个 1 分到每一份，然后再把剩下的 n-k 分成 k 份即可，分法有: $dp[n-k][k]$

     b.至少有一份为 1 的分法，可以先那出一个 1 作为单独的1份，剩下的 n- 1 再分成 k- 1 份即可，分法有：$dp[n-1][k-1]$

   • 五边形数

     $$P(n)=\sum\limits_{k=1}^n(-1)^{k+1}\left [P\left( n-\frac 1 2k(3k-1)\right)+P\left( n-\frac 1 2k(3k+1)\right )\right ]$$

     其中 $n-\frac 1 2k(3k-1)>=0 , n-\frac 1 2k(3k+1)>=0$

2. 将n划分成若干正整数(不同)之和的划分数

   • O(n^2)

     $dp[n][m]= dp[n][m-1]+ dp[n-m][m-1]$,$dp[n][m]$表示整数 n 的划分中，每个数不大于m且两两不同的划分数。 划分分为两种情况：

     a.划分中每个数都小于m,相当于每个数不大于 m-1,划分数为 $dp[n][m-1]$.
     b.划分中有一个数为 m.在n中减去m,剩下相当对n-m进行划分,并且每一个数不大于m-1，故划分数为 $dp[n-m][m-1]$

• O(n^1.5)

  $dp[i][j]$代表i个数相加等于j ,由于有不同的数最多不超过O(sqrt(n))算法个，则可以优化 $dp[i][j] = dp[i-1][j-i] + dp[i][j-i]$ ,这个递推方程的转移含义是i-1个数每个数都加1,最后再添上一个1，就从$dp[i-1][j-i]$转到$dp[i][j]$,还有就是i个数每个数都加1，就从$dp[i][j-i]$转到$dp[i][j]$

1. 将n划分成若干奇正整数之和的划分数.

   • 法一

     $g[i][j]$: 将i划分为j个偶数

     $f[i][j]$: 将i划分为j个奇数

     $g[i][j] = f[i - j][j]$;

     $f[i][j] = f[i - 1][j - 1] + g[i - j][j]$;

   • 法二

     奇数和$dp[i][j]$为划分最大的奇数不超过j的数，$dp[i][j]=dp[i-n(j)][j]+dp[i][j-2]$,n(j)为不超过j的最大奇数,初始条件，$dp[i][1]=1$,j为偶数时候$dp[i][j]=dp[i][j-1]$;当i==n(j) ,出现$dp[0][j]$，也就是当i为奇数时候，$dp[0][j]=1$;

2. 其他

   n的分拆数中最大部分为m的个数=把n分拆成m部分的个数

   

   n的分拆数中每一部分都小于等于m的个数=把n分成m份或更小
   n的分拆数中每部分的数两两不同的个数等于全部都是奇数的个数
   n的分拆数中每部分的数都相等的个数=n 的因子个数
   n的分拆数中每部分都是1或2(或者把n分拆成1或2部分)的个数=floor(n/2+1)
   n的分拆数中每部分都是1或2或3(或者把n分拆成1或2或3部分)的个数=(n+3)^2/12
   n的分拆数中每部分都>=k的个数为：f[n]-f[n-1]-f[n-2] · · · +f[n-k] (f[n]表示将n划分成若干正整数(可相同)之和的划分数)

约瑟夫环

令$f[1]$表示$n-m+1$个人玩游戏报 k 退出，第 1 个人退出的编号(下标从 0 开始)；

$f[m]$表示 n 个人玩游戏报 k 退出，第 m 个人退出的编号，

注:当 m<<n 的时候，可以用等差数列的性质，加速运算.

$$f(x)=\left\{\begin{array}{c} (k-1)\%(n-m+1)&x=1\\ (f(x-1)+k)\%(n-m+x)&1<x\le m \end{array} \right.$$

高斯消元

整数

const int MAXN = 300; //有 equ 个方程，var 个变元。增广矩阵行数为 equ, 列数为 var+1, 分别为 0 到 var 
int equ,var; 
int a[MAXN][MAXN]; //增广矩阵 
int x[MAXN]; //解集 
int free_x[MAXN];//用来存储自由变元（多解枚举自由变元可以使用） 
int free_num;//自由变元的个数

//返回值为 -1 表示无解，为 0 是唯一解，否则返回自由变元个数 
int Gauss(){
    int max_r,col,k; 
    free_num = 0; 
    for(k = 0, col = 0 ; k < equ && col < var ; k++, col++){ 
        max_r = k; 
        for(int i = k+1;i < equ;i++){ 
            if(abs(a[i][col]) > abs(a[max_r][col])) max_r = i; 
        } 
        if(a[max_r][col] == 0){ 
            k--; 
            free_x[free_num++] = col;//这个是自由变元 
            continue; 
        } 
        if(max_r != k){ 
            for(int j = col; j < var+1; j++) swap(a[k][j],a[max_r][j]); 
        } 
        for(int i = k+1;i < equ;i++){ 
            if(a[i][col] != 0){ 
                for(int j = col;j < var+1;j++)a[i][j] ^= a[k][j];
            }
        }
    }
    for(int i = k;i < equ;i++) if(a[i][col] != 0)return -1;//无解 
    if(k < var) return var-k;//自由变元个数 
    //唯一解，回代 
    for(int i = var-1; i >= 0;i--){
        x[i] = a[i][var];
        for(int j = i+1;j < var;j++)
        x[i] ^= (a[i][j] && x[j]); 
    }
    return 0;
}

实数

//1
const int num = 100; 
double matrix[num][num + 1]; //系数矩阵， 从 0 开始
double ans[num]; //结果数组
void exchange_col(int p1,int p2,int n) {//交换 p1 行和 p2 行的所有数据 
    double t; 
    int i;
    for(i = 0 ; i <= n ; i++) t = matrix[p1][i],matrix[p1][i] = matrix[p2][i],matrix[p2][i] = t;
}
bool gauss(int n){ //求解系数矩阵为 n 的线性方程组
    int i,j,k,p; 
    double r;
    for(i = 0 ; i < n - 1 ; i++) {
        p = i; 
        for(j = i + 1 ; j < n ; j++) { //寻找 i 列绝对值最大值位置 
            if(abs(matrix[j][i]) > abs(matrix[p][i])) p = j; 
        }
        if(p != i) exchange_col(i,p,n); 
        if(matrix[i][i] == 0) return false; 
        for(j = i + 1 ; j < n ; j++) { //剩余列进行消元 
            r = matrix[j][i] / matrix[i][i]; 
            for(k = i ; k <= n ; k++) matrix[j][k] -= r * matrix[i][k]; 
        }
    }
    for(i = n - 1 ; i >= 0 ; i--) { //获得结果
        ans[i] = matrix[i][n];
        for(j = n - 1 ; j > i ; j--) ans[i] -= matrix[i][j] * ans[j]; 
        if(matrix[i][i] == 0) return false;
        ans[i] /= matrix[i][i]; 
    }
    return true;
}

//kuangbin
#define eps 1e−9 
const int MAXN=220; 
double a[MAXN][MAXN],x[MAXN];//方程的左边的矩阵和等式右边的值, 求解之后 x 存的就是结果 
int equ,var;//方程数和未知数个数 
int Gauss(){ /* * 返回 0 表示无解,1 表示有解 */ 
    int i,j,k,col,max_r;
    for(k=0,col=0;k<equ&&col<var;k++,col++){
        max_r=k;
        for(i=k+1;i<equ;i++) if(fabs(a[i][col])>fabs(a[max_r][col])) max_r=i;
        if(fabs(a[max_r][col])<eps)return 0;
        if(k!=max_r){ 
            for(j=col;j<var;j++) swap(a[k][j],a[max_r][j]); 
            swap(x[k],x[max_r]);
        }
        x[k]/=a[k][col];
        for(j=col+1;j<var;j++)a[k][j]/=a[k][col];
        a[k][col]=1;
        for(i=0;i<equ;i++)
        if(i!=k){
            x[i]-=x[k]*a[i][col]; 
            for(j=col+1;j<var;j++)a[i][j]-=a[k][j]*a[i][col]; 
            a[i][col]=0; 
        }
    }
    return 1;
}

母函数

$$(1+a_1x)(1+a_2x) · · · (1+a_nx) =\\ 1+(a_1+a_2+ · · · +a_n)x+(a_1a_2+a_1a_3+ · · · +a_{n-1}a_n)x^2+ · · · +a_1a_2 · · · a_nx^n$$

可以看出：$x^2$项的系数$a_1a_2+a_1a_3+ · · · +a_{n-1}a_n$中所有的项包括n个元素$a_1,a_2,· · · ,a_n$中取两个组合的全体，同理：$x^3$项系数包含了从n个元素$a_1,a_2,· · · ,a_n$中取3个元素组合的全体，以此类推。

栗子：
整数拆分
由a,b,c,d….组成n

$G(x)=(1+x^a+x^{2a}+x^{3a}+ · · · )(1+x^b+x^{2b}+x^{3b}+ · · · ) · · · =q_0+q_1x+q_2x^2+q_3x^3+ · · · $

$q_i$表示由a,b,c,d….组成 $i$ 的种类

for (int i=0;i<=n;i++){
	ans[i]=1;
	tmp[i]=0;
}//element[0]=1时这样写，视情况而定
for (int i=1; i<elememt_count; i++){
	for (int j=0;j<=n;j++)
	for (int k=0;k+j<=n; k+=element[i]) tmp[k+j]=(tmp[k+j]+ans[j])%mod;
	for (int j=0;j<=n;j++){
		ans[j]=tmp[j];
		tmp[j]=0;
	}
}

• 分治FFT

  计算 $\large\prod(1+a_ix)$ ，复杂度 $\mathcal{O}(n*log^2n)$

  ll qPow(ll a, ll b, ll c) {  //求(a^b) % c
      ll ret = 1;
      while (b) {
          if (b & 0x1) ret = ret * a % c;
          a = a * a % c;
          b >>= 1;
      }
      return ret;
  }
  typedef complex <db> cp;
  int nn;
  cp A[N],omg[M], inv[M];
  void init(int a,int b){
      nn=1;
      while(nn < a + b) nn <<= 1;
      for(int i = 0; i < nn; i++){
          omg[i] = cp(cos(2 * pi * i / nn), sin(2 * pi * i / nn));
          inv[i] = conj(omg[i]);
      }
  }
  void fft(cp *a,cp  *type){//type=omg 系数->点值 ; type=inv 点值->系数
      int lim = 0;
      while((1 << lim) < nn) lim++;
      for(int i = 0; i < nn; i++){
          int t = 0;
          for(int j = 0; j < lim; j++)
              if((i >> j) & 1) t |= (1 << (lim - j - 1));
          if(i < t) swap(a[i], a[t]); // i < t 的限制使得每对点只被交换一次（否则交换两次相当于没交换）
      }
      for(int l = 2; l <= nn; l *= 2){
          int m = l / 2;
          for(cp *p = a; p != a + nn; p += l)
          for(int i = 0; i < m; i++){
              cp t = type[nn / l * i] * p[i + m];
              p[i + m] = p[i] - t;
              p[i] += t;
          }
      }
  }
  int a[M]; 
  int main(){
      int n,trn=1;
      in(n);
      while(trn<n)trn<<=1;
      for0(i,n){
          in(a[i]);
      }
      for0(i,trn){
          A[4*i].real(1);
          A[4*i+1].real(a[i]);
      }
      int len=2;
      for(int i=4;i<4*trn;i<<=1){
          init(len,len);
          int nxt_len=2*len-1;
          for(int j=0;j<4*trn;j+=2*i){
              for(int k=len;k<nn;k++)A[j+k]=A[j+k+i]=cp(0,0);
              fft(A+j,omg);fft(A+j+i,omg);
              for(int k=0;k<nn;k++)A[j+k]*=A[j+k+i];
              fft(A+j,inv);
              for(int k=0;k<nxt_len;k++)A[j+k]=cp(ll(A[j+k].real()/nn+0.1)%mod,0);
          }
          len=nxt_len;
      }
      for(int i=0;i<=n;i++){
          a[i]=A[i].real()+0.1;
          out(a[i],0);
      }
  }

傅立叶变换

DFT

离散傅立叶变换 $O(n^2)$

• 复数相乘的规则：模长相乘，幅角相加。
• 取n个复数，单元根$\large \omega_n=\cos\frac{2\pi}n+i\sin\frac{2\pi}n=e^{i\frac {2\pi}n}$，取$\omega_n^0,\omega_n^1,\omega_n^2,\cdots,\omega_n^{n-1}$
  • $\omega_{2n}^{2k}=\omega_{n}^{k}$
  • $\omega_{n}^{k+\frac n2}=-\omega_{n}^{k}$
  • $\omega_n^0=\omega_n^n=1$
• 设$(y_0,y_1,\cdots,y_{n-1})$为多项式$A(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}$的离散傅立叶变换，设多形式$B(x)=y_0+y_1x+\cdots+y_{n-1}x^{n-1}$，取$\omega_n^0,\omega_n^{-1},\omega_n^{-2},\cdots,\omega_n^{-(n-1)}$代入$B(x)$得$(z_0,z_1,\cdots,z_{n-1})$，易证$a_i=\frac {z_i} n$

FFT

快速傅立叶变换，在$O(n log n)$的时间内将一个多项式转换成它的点值表示(或相反)的算法

设$A(x) = a_0 + a_1x + a_2x^2 +…+a_{n-1}x^{n-1}$，将$A(x)$拆分成奇偶两部分：

$$A(x) = (a_0 + a_2x^2 + ... + a_{n - 2}x^{n - 2}) + (a_1x + a_3x^3 + ... + a_{n-1}x^{n-1})$$

设：

$$A_1(x) = a_0 + a_2x + ... + a_{n - 2}x^{\frac{n}{2} - 1}\\ A_2(x) = a_1 + a_3x + ... + a_{n - 1}x^{\frac{n}{2} - 1}\\ A(x)=A_1(x^2)+xA_2(x^2)$$

当$k<\frac n 2$时，

$$A(\omega_n^k)= A_1(\omega_n^{2k}) + \omega_n^kA_2(\omega_n^{2k})= A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$$

当$>\frac n 2$时，

$$A(\omega_n^{k + \frac{n}{2}}) = A_1(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}}A_2(\omega_n^{2k + n}) \\ = A_1(\omega_{\frac{n}{2}}^{k} \times \omega_n^n) + \omega_n^{k + \frac{n}{2}} A_2(\omega_{\frac{n}{2}}^{k} \times \omega_n^n) \\ = A_1(\omega_{\frac{n}{2}}^{k}) - \omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$$

蝴蝶操作

typedef complex <double> cp;
int n;
cp a[M], b[M], omg[M], inv[M];
void init(int a,int b){
    n=1;
    while(n < a + b) n <<= 1;
    for(int i = 0; i < n; i++){
        omg[i] = cp(cos(2 * pi * i / n), sin(2 * pi * i / n));
        inv[i] = conj(omg[i]);
    }
}
void fft(cp *a, cp *type){//type=omg 系数->点值 ; type=inv 点值->系数
    int lim = 0;
    while((1 << lim) < n) lim++;
    for(int i = 0; i < n; i++){
        int t = 0;
        for(int j = 0; j < lim; j++)
            if((i >> j) & 1) t |= (1 << (lim - j - 1));
        if(i < t) swap(a[i], a[t]); // i < t 的限制使得每对点只被交换一次（否则交换两次相当于没交换）
    }
    for(int l = 2; l <= n; l *= 2){
        int m = l / 2;
        for(cp *p = a; p != a + n; p += l)
        for(int i = 0; i < m; i++){
            cp t = type[n / l * i] * p[i + m];
            p[i + m] = p[i] - t;
            p[i] += t;
        }
    }
}
//2
typedef complex <double> cp;
struct FFT{
    int n, R[N];
    cp a[N], b[N];
    const double pi=acos(-1);
    void init(int bound){//bound是积多项式的最高次幂
        int L(0);
        for(n=1;n<=bound;n<<=1,L++);
        for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1)),a[i]=b[i]=0;
    }
    void fft(cp* a, int opt){
        int i, j, k;
        cp wn, w, x, y;
        for(i=0;i<n;i++)if(i>R[i])swap(a[i],a[R[i]]);
        for(i=1;i<n;i<<=1){
            wn=cp(cos(pi/i),opt*sin(pi/i));
            for(j=0;j<n;j+=i<<1)
                for(w=cp(1,0),k=0;k<i;k++,w=w*wn){
                    x=a[k+j], y=a[k+j+i]*w;
                    a[k+j]=x+y, a[k+j+i]=x-y;
                }
        }
        if(opt==-1)for(i=0;i<n;i++)a[i]/=n;
    }
    void run(){
        fft(a,1), fft(b,1);
        for(int i=0;i<n;i++)a[i]*=b[i];
        fft(a,-1);
    }
}fft;

NTT

• 特定模数

  将单元根$\omega _n=e^{i\frac {2\pi}n} $变为$g^{\frac {P-1} n}$，需要满足模数 P 是质数且 P−1 是 n (2的次幂) 的倍数。即：

  $$P=r\cdot2^k+1\quad(n\le2^k)$$

  显然对$\omega _n=g^{\frac {P-1} n}$有$\omega_{2n}^{2k}=\omega_{n}^{k}$，而且由于$g^{\frac {P-1} {2n}}\equiv -1 (mod P)$所以$\omega_{n}^{k+\frac n2}=-\omega_{n}^{k}$

int n,r[N];
ll g=3,inv_g,inv_n,mod=998244353;
ll qPow(ll a, ll b, ll c) {  //求(a^b) % c
    ll ret = 1;
    while (b) {
        if (b & 0x1) ret = ret * a % c;
        a = a * a % c;
        b >>= 1;
    }
    return ret;
}
void init(int a,int b){
    int L=0;n=1;
    while(n < a + b) n <<= 1,L++;
    inv_n = qPow(n,mod-2,mod);
    inv_g = qPow(g,mod-2,mod);
    for(int i = 0; i < n; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));  
}
inline void NTT(int *A, int type) {//type=1 系数->点值 ; type=-1 点值->系数
    for(int i = 0; i < n; i++) if(i < r[i]) swap(A[i], A[r[i]]);
    for(int mid = 1; mid < n; mid <<= 1) {
        ll Wn = qPow( type == 1 ? g : inv_g , (mod - 1) / (mid << 1),mod);
        for(int j = 0; j < n; j += (mid << 1)) {
            ll w = 1;
            for(int k = 0; k < mid; k++, w = (w * Wn) % mod) {
                int x = A[j + k], y = w * A[j + k + mid] % mod;
                A[j + k] = (x + y) % mod,
                A[j + k + mid] = (x - y + mod) % mod;
            }
        }
    }
    if(type==-1) for0(i,n)A[i]=A[i]*inv_n%mod;
}
int main(){
    init(len_a,len_b);
    NTT(a, 1); NTT(b, 1);
    for(int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%mod;
    NTT(a, -1); //结果放入a中,有效值范围为[0,len_a+len_b-2]
}

• 三模数

  两个多项式相乘后，每位的答案一定小于1e9*1e9*n，即1e24，所以我们选3个乘积大于1e24的NTT模数，分别做一次，得到每个位上模意义下的答案，再做CRT，先将两个模数意义下的答案合并，得到两个模数，设模数为P1(longlong) ,P2(int), 余数为a1,a2。设答案$ANS=P1\times K+a1=P2\times Q+a2$，那么$K\times P1=P2\times Q+(a2-a1)$，$K\times P1 \% P2=a2-a1$，$a1-a2$为常数，用同余方程的解法即可解出K模P2(int)意义下的值，又有$ANS<P1\times P2$，所以 $K\times P1+a1<P1\times P2$，显然$K<P2$，所以原本答案K的值只能为模P2意义下的值。

#include<bits/stdc++.h>
using namespace std;
#define RI register int
int read() {
    int q=0;char ch=' ';
    while(ch<'0'||ch>'9') ch=getchar();
    while(ch>='0'&&ch<='9') q=q*10+ch-'0',ch=getchar();
    return q;
}
typedef long long LL;
const int N=262150,mm[3]={7*(1<<26)+1,998244353,479*(1<<21)+1},G=3;
int n,m,kn=1,len,mod;
int a[N],b[N],k1[N],k2[N],ans[3][N],rev[N];
int ksm(int x,int y,int p) {
    int re=1;
    for(;y;y>>=1,x=1LL*x*x%p) if(y&1) re=1LL*re*x%p;
    return re;
}
LL ksc(LL x,LL y,LL p) {return (x*y-(LL)((long double)x/p*y+1e-8)*p+p)%p;}
void NTT(int *a,int n,int p,int x) {
    for(RI i=0;i<n;++i) if(rev[i]>i) swap(a[i],a[rev[i]]);
    for(RI i=1;i<n;i<<=1) {
        int gn=ksm(G,(p-1)/(i<<1),p);
        for(RI j=0;j<n;j+=i<<1) {
            int g=1,t1,t2;
            for(RI k=0;k<i;++k,g=1LL*g*gn%p) {
                t1=a[j+k]%p,t2=1LL*g*a[j+i+k]%p;
                a[j+k]=(t1+t2)%p,a[j+i+k]=(t1-t2+p)%p;
            }
        }
    }
    if(x==1) return;
    int inv=ksm(n,p-2,p);reverse(a+1,a+n);
    for(RI i=0;i<n;++i) a[i]=1LL*a[i]*inv%p;
}
void work(int o) {
    for(RI i=0;i<kn;++i) k1[i]=a[i],k2[i]=b[i];
    NTT(k1,kn,mm[o],1),NTT(k2,kn,mm[o],1);
    for(RI i=0;i<kn;++i) ans[o][i]=1LL*k1[i]*k2[i]%mm[o];
    NTT(ans[o],kn,mm[o],-1);
}
int main(){
    n=read(),m=read(),mod=read();
    for(RI i=0;i<=n;++i) a[i]=read();
    for(RI i=0;i<=m;++i) b[i]=read();
    while(kn<=n+m) kn<<=1,++len;
    for(RI i=1;i<kn;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
    for(RI i=0;i<3;++i) work(i);
    LL M=1LL*mm[0]*mm[1];
    LL kl1=ksc(mm[1],ksm(mm[1]%mm[0],mm[0]-2,mm[0]),M);
    LL kl2=ksc(mm[0],ksm(mm[0]%mm[1],mm[1]-2,mm[1]),M);
    for(RI i=0;i<=n+m;++i) {
        int t0=ksm(ans[0][i],mm[1]-2,mm[1]),t1=ksm(ans[1][i],mm[0]-2,mm[0]);
        LL A=(ksc(kl1,ans[0][i],M)+ksc(kl2,ans[1][i],M))%M;
        LL k=((ans[2][i]-A)%mm[2]+mm[2])%mm[2]*ksm(M%mm[2],mm[2]-2,mm[2])%mm[2];
        printf("%lld ",((M%mod)*(k%mod)%mod+A%mod)%mod);
    }
    return 0;
}

• MTT

  将多项式的每个系数拆成$k\times M+b$的形式，其中 M 为常数,考虑卷积过程：

  $$(k_1M+b_1)*(k_2M+b_2)=k_1k_2M^2+b_1k_2M+b_2k_1M+b_1b_2$$

  当$M=\sqrt P$时,这些东西都是 P 的级别，所以FFT 出来的结果都是1e14级别的。所以需要4次 FFT 和3次 IFFT。进一步优化，定义

  $$P(x)=A(x)+B(x)i,Q(x)=A(x)−B(x)i$$

  则会有

  $$Q(\omega _n^k)=conj(P(\omega_n^{-k}))=conj(P(\omega_n^{n-k}))$$

  于是我们只用DFT这个P(x)，就能得到Q(x)被DFT之后的结果,还可以求得 A，B 的 FFT。这样我们就可以处理出$a=k_1k_2,b=k_1b_2,c=k_2b_1,d=b_1b_2$的值，然后构造M1=a+bi和M2=c+di两个多项式，IDFT回去。则分别取这两个IDFT后的多项式的实部和虚部（第一个的实部为a，虚部为b以此类推），就是IDFT后的a,b,c,d节约2的常数，这样只要做四次DFT即可。

#include<bits/stdc++.h>
using namespace std;
#define RI register int
int read() {
    int q=0;char ch=' ';
    while(ch<'0'||ch>'9') ch=getchar();
    while(ch>='0'&&ch<='9') q=q*10+ch-'0',ch=getchar();
    return q;
}
typedef long long LL;
const int N=262150,mm[3]={7*(1<<26)+1,998244353,479*(1<<21)+1},G=3;
int n,m,kn=1,len,mod;
int a[N],b[N],k1[N],k2[N],ans[3][N],rev[N];
int ksm(int x,int y,int p) {
    int re=1;
    for(;y;y>>=1,x=1LL*x*x%p) if(y&1) re=1LL*re*x%p;
    return re;
}
LL ksc(LL x,LL y,LL p) {return (x*y-(LL)((long double)x/p*y+1e-8)*p+p)%p;}
void NTT(int *a,int n,int p,int x) {
    for(RI i=0;i<n;++i) if(rev[i]>i) swap(a[i],a[rev[i]]);
    for(RI i=1;i<n;i<<=1) {
        int gn=ksm(G,(p-1)/(i<<1),p);
        for(RI j=0;j<n;j+=i<<1) {
            int g=1,t1,t2;
            for(RI k=0;k<i;++k,g=1LL*g*gn%p) {
                t1=a[j+k]%p,t2=1LL*g*a[j+i+k]%p;
                a[j+k]=(t1+t2)%p,a[j+i+k]=(t1-t2+p)%p;
            }
        }
    }
    if(x==1) return;
    int inv=ksm(n,p-2,p);reverse(a+1,a+n);
    for(RI i=0;i<n;++i) a[i]=1LL*a[i]*inv%p;
}
void work(int o) {
    for(RI i=0;i<kn;++i) k1[i]=a[i],k2[i]=b[i];
    NTT(k1,kn,mm[o],1),NTT(k2,kn,mm[o],1);
    for(RI i=0;i<kn;++i) ans[o][i]=1LL*k1[i]*k2[i]%mm[o];
    NTT(ans[o],kn,mm[o],-1);
}
int main()
{
    n=read(),m=read(),mod=read();
    for(RI i=0;i<=n;++i) a[i]=read();
    for(RI i=0;i<=m;++i) b[i]=read();
    while(kn<=n+m) kn<<=1,++len;
    for(RI i=1;i<kn;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
    for(RI i=0;i<3;++i) work(i);
    LL M=1LL*mm[0]*mm[1];
    LL kl1=ksc(mm[1],ksm(mm[1]%mm[0],mm[0]-2,mm[0]),M);
    LL kl2=ksc(mm[0],ksm(mm[0]%mm[1],mm[1]-2,mm[1]),M);
    for(RI i=0;i<=n+m;++i) {
        int t0=ksm(ans[0][i],mm[1]-2,mm[1]),t1=ksm(ans[1][i],mm[0]-2,mm[0]);
        LL A=(ksc(kl1,ans[0][i],M)+ksc(kl2,ans[1][i],M))%M;
        LL k=((ans[2][i]-A)%mm[2]+mm[2])%mm[2]*ksm(M%mm[2],mm[2]-2,mm[2])%mm[2];
        printf("%lld ",((M%mod)*(k%mod)%mod+A%mod)%mod);
    }
    return 0;
}

FWT

快速沃尔什变换，利用类似 $FFT$ 的东西解决多项式的位运算卷积

定义：

$$A+B=(A_0+B_0,A_1+B_1,\cdots\cdots)\\ A-B=(A_0-B_0,A_1-B_1,\cdots\cdots)\\ A\oplus B=\left(\sum_{i \oplus j=0} A_{i} * B_{j}, \sum_{i \oplus j=1} A_{i} * B_{j}, \cdots\cdots\right)\\ A\times B=(A_0*B_0,A_1*B_1,\cdots\cdots)\\$$

寻找 $FWT(A\oplus B)=FWT(A)\times FWT(B)$

• or

  $$C_k=\sum_{i|j=k}A_i*B_j$$

  对于长度为 $2^n$ 的多项式，$A_0$ 是前 $2^{n-1}$ 前项，$A_1$ 是后 $2^{n-1}$ 前项

  $$F W T(A)=\left\{\begin{array}{ll}{\left(\ F W T\left(A_{0}\right), F W T\left(A_{0}+A_{1}\right)\ \right)} & {n>0} \\ {A} & {n=0}\end{array}\right.\\ IFWT(A)=(\ IFWT(A_0),IFWT(A_1)-IFWT(A_0)\ )$$

• and

  $$C_k=\sum_{i\&j=k}A_i*B_j$$

  对于长度为 $2^n$ 的多项式，$A_0$ 是前 $2^{n-1}$ 前项，$A_1$ 是后 $2^{n-1}$ 前项

  $$F W T(A)=\left\{\begin{array}{ll}{\left(\ F W T\left(A_{0}+A_1\right), F W T\left(A_{1}\right)\ \right)} & {n>0} \\ {A} & {n=0}\end{array}\right.\\ IFWT(A)=(\ IFWT(A_0)-IFWT(A_1),IFWT(A_1)\ )$$

• xor

  $$C_k=\sum_{i\land j=k}A_i*B_j$$

  对于长度为 $2^n$ 的多项式，$A_0$ 是前 $2^{n-1}$ 前项，$A_1$ 是后 $2^{n-1}$ 前项

  $$F W T(A)=\left\{\begin{array}{ll}{\left(\ F W T(A_0)+FWT(A_1), F W T(A_0)-FWT(A_1)\ \right)} & {n>0} \\ {A} & {n=0}\end{array}\right.\\ IFWT(A)=(\frac {IFWT(A_0)+IFWT(A_1)}2,\frac {IFWT(A_0)-IFWT(A_1)}2)$$

  ll qPow(ll a, ll b, ll c) {  //求(a^b) % c
      ll ret = 1;
      while (b) {
          if (b & 0x1) ret = ret * a % c;
          a = a * a % c;
          b >>= 1;
      }
      return ret;
  }
  int n,inv_2;
  void init(int maxx){//maxx=max(len_a,len_b)
      n=1;
      while(n<maxx)n<<=1;
      inv_2=qPow(2,mod-2,mod);
  } 
  void FWT_or(int *a,int opt){//type=1 系数->点值 ; type=-1 点值->系数
      for(int i=1;i<n;i<<=1)
          for(int p=i<<1,j=0;j<n;j+=p)
              for(int k=0;k<i;++k)
                  if(opt==1)a[i+j+k]=(a[j+k]+a[i+j+k])%mod;
                  else a[i+j+k]=(a[i+j+k]+mod-a[j+k])%mod;
  }
  void FWT_and(int *a,int opt){
      for(int i=1;i<n;i<<=1)
          for(int p=i<<1,j=0;j<n;j+=p)
              for(int k=0;k<i;++k)
                  if(opt==1)a[j+k]=(a[j+k]+a[i+j+k])%mod;
                  else a[j+k]=(a[j+k]+mod-a[i+j+k])%mod;
  }
  void FWT_xor(int *a,int opt){
      for(int i=1;i<n;i<<=1)
          for(int p=i<<1,j=0;j<n;j+=p)
              for(int k=0;k<i;++k){
                  int X=a[j+k],Y=a[i+j+k];
                  a[j+k]=(X+Y)%mod;
                  a[i+j+k]=(X+mod-Y)%mod;
                  if(opt==-1) a[j+k]=1ll*a[j+k]*inv_2%mod , a[i+j+k]=1ll*a[i+j+k]*inv_2%mod;
              }
  }

FMT

快速莫比乌斯变化

阶乘最末尾非零

设 $G(n)$ 为除去 $5$ 的倍数的前 $n$ 项乘积的末尾位

$\frac 2 2=6, \frac 4 2=2, \frac 6 2=8, \frac 8 2=4$

$$\large F(n)=F(\lfloor\frac n 5\rfloor)\frac {G(n)} {2^{\lfloor\frac n 5\rfloor\%4}}$$

$\huge\frac {G(n)} {2^{\lfloor\frac n 5\rfloor\%4}}$有一个长度为20的循环节6,2,6,4,2,2,4,2,8,4,4,8,4,6,8,8,6,8,2,6

//普通版
int G[10]={4,4,8,4,6,6,6,2,6,4},G5[5]={0,1,2,6,4},ce[4]={6,8,4,2};
int f(int n){
    if(n<5)return G5[n];
    return G[(n-5)%10]*ce[n/5%4]*f(n/5);
}
int main() {
    int n;
    ind(n);
    printf("%d\n",f(n)%10);
    return 0;
}
//高精度
#include <stdio.h>
#include <string.h>
#define MAXN 10000
int lastdigit(char* buf) {
    const int mod[20] = {1, 1, 2, 6, 4, 2, 2, 4, 2, 8, 4, 4, 8, 4, 6, 8, 8, 6, 8, 2};
    int len = strlen(buf), a[MAXN], i, c, ret = 1;
    if (len == 1) return mod[buf[0] - '0'];
    for (i = 0; i < len; i++) a[i] = buf[len - 1 - i] - '0';
    for (; len; len -= !a[len - 1]) {
        ret = ret * mod[a[1] % 2 * 10 + a[0]] % 5;
        for (c = 0, i = len - 1; i >= 0; i--)c = c * 10 + a[i], a[i] = c / 5, c %= 5;
    }
    return ret + ret % 2 * 5;
}
int main() {
    char n[10000];
    while (scanf("%s", n) != EOF) printf("%d/n", lastdigit(n));
}

二次剩余

当存在某个 $X$，式子 $X^2 = n (\% p)$ 成立时，称 $n$ 是模 $p$ 的二次剩余
当对任意 $X$ 不成立时，称 $n$ 是模 $p$ 的二次非剩余

• 定理1: 对于 $x^2 = n (\% p)$ 总共有 $\frac {p-1} 2$ 种 n 能使方程有解(排除n=0)

• 引理: p为奇素数

  $$n^{\frac {p-1} 2}\equiv\pm1(\%\ p)$$

• 定理2: 勒让德符号(Legendre symbol)

  $$\left (\frac q p\right)=q^{\frac {p-1}2}(\%\ p)= \left\{ \begin{array}{c} 1&q在模p 意义下是二次剩余\\ -1&q在模p 意义下是非二次剩余\\ 0&q\equiv0(\%\ p) \end{array} \right.$$

• 二次互反律

  设 p,q 是不同的奇素数, 则

  $$\left (\frac {-1} p\right)= \left\{ \begin{array}{c} 1&p\equiv1(\%\ 4)\\ -1&p\equiv3(\%\ 4) \end{array} \right.\\ \left (\frac {2} p\right)= \left\{ \begin{array}{c} 1&p\equiv1\ or\ 7(\%\ 8)\\ -1&p\equiv3\ or \ 5(\%\ 8) \end{array} \right.\\ \left (\frac q p\right)= \left\{ \begin{array}{c} \left (\frac p q\right)&p\equiv1(\%\ 4)\ or\ q\equiv1(\%\ 4)\\ -\left (\frac p q\right)&p\equiv3(\%\ 4)\ or\ q\equiv3(\%\ 4) \end{array} \right.\\$$

• 广义二次互反律

  设 a,b 为正奇数, 则

  $$\left (\frac {-1} b\right)= \left\{ \begin{array}{c} 1&b\equiv1(\%\ 4)\\ -1&b\equiv3(\%\ 4) \end{array} \right.\\ \left (\frac {2} b\right)= \left\{ \begin{array}{c} 1&b\equiv1\ or\ 7(\%\ 8)\\ -1&b\equiv3\ or \ 5(\%\ 8) \end{array} \right.\\ \left (\frac a b\right)= \left\{ \begin{array}{c} \left (\frac b a\right)&a\equiv1(\%\ 4)\ or\ b\equiv1(\%\ 4)\\ -\left (\frac b a\right)&a\equiv b\equiv3(\%\ 4) \end{array} \right.\\$$

• Cipolla算法

  求解 $x^2 = n (\% p)$
  当 $\omega=a^2-n$ 不是模 p 的二次剩余 $x=(a+\sqrt {\omega})^{\frac {p+1}2}$

  证明:

$$(a+\sqrt{\omega})^p=a^p+\omega^{\frac {p-1}2}\sqrt{\omega}=a-\sqrt{\omega}\\ (a+\sqrt{\omega})^{p+1}\equiv(a+\sqrt{\omega})^p(a+\sqrt{\omega}) \equiv(a-\sqrt{\omega})(a+\sqrt{\omega})\equiv a^2-\omega\equiv n\ (\%\ p)$$

//方案1:
ll sqrtMod(ll a, ll p){ // x^2 == a (mod P)
    ll b,k,i;
    if(p==2)return a%p;
    if(qPow(a,(p-1)/2,p)==1){
        if(p%4==3)return qPow(a,(p+1)/4,p);
        for(b=1;qPow(b,(p-1)/2,p)==1;b++);
        i=(p-1)/2;k=0;
        do{
            i/=2;
            k/=2;
            if((qPow(a,i,p)*qPow(b,k,p)+1)%p==0)k+=(p-1)/2;
        }while(i%2==0);
        return (qPow(a,(i+1)/2,p)*qPow(b,k/2,p))%p;
    }
    return -1;
}

//方案2:
ll w, P;
struct Complex {
    ll r, i;
    Complex(ll _r = 0, ll _i = 0):r(_r), i(_i) {}
    Complex operator +(ll x) {
        return Complex((x + r) % P, i);
    }
    Complex operator *(Complex a) {
        ll rr = (a.r * r % P + w * a.i % P * i % P + P * 3) % P,
        ii = (a.i * r % P + a.r * i % P) % P;
        return Complex(rr, ii);
    }
};
ll complexPow(Complex a, ll b) {
    Complex c(1);
    while (b) {
        if (b & 1) c = c * a;
        a = a * a; b >>= 1;
    }
    return c.r;
}
ll Solve(ll a) { // x^2 == a (mod P)
    if (qPow(a, (P - 1) / 2 , P) == P - 1) return -1;
    if (qPow(a, (P - 1) / 2 ,P) == 0) return 0;
    ll tmp;
    Complex k(0, 1);
    while (true) {
        tmp = rand();
        if (qPow((tmp * tmp - a + P) % P, (P - 1) / 2,P) == P - 1) break;
    }
    w = (tmp * tmp - a + P) % P;
    return complexPow(k + tmp, (P + 1) / 2);
}

• 求解 $x^2 = n (\% p^a)$ , 其中 $gcd(n , p) = 1$, 并且 p 为奇质数

  若 r 是 $x^2 = n (\% p)$ 的一个解

  $$r^2-n=kp\Rightarrow (r^2-n)^a\equiv k^ap^a\equiv 0(\%\ p^a)\\ 设\ (r+\sqrt{n})^a=t+u\sqrt{n}\ ,\ (r-\sqrt{n})^a=t-u\sqrt{n}\\ (r^2-n)^a\equiv t^2-u^2n\equiv 0(\%\ p^a)\\ 可以证明\ gcd(t,p)=1\ ,\ gcd(u,p)=1\\ t^2*Inv(u^2)\equiv n(\%\ p^a)\\ x=t*Inv(u)$$

• 二次同余方程模合数

  $$x^2\equiv n(\%\ y)\Longrightarrow \left\{ \begin{array}{c} x^2\equiv n(\%\ p_1^{a_1})\\ x^2\equiv n(\%\ p_2^{a_2})\\ \cdots\\ x^2\equiv n(\%\ p_m^{a_m})\\ \end{array} \right.\\ y=\prod_{i=1}^mp_i^{a_i}$$

  求出 $x_1,x_2,\cdots,x_m$ 再用中国剩余定理求解

阶

定义

设 $m>1 , (a,m)=1$，满足 $a^x\equiv 1(\mod m)$ 的最小的 $x$，称为 $a$ 对 $m$ 的阶，记为 $\delta_m(a)$

定理与推论

1. 若 $m>1,(a,m)=1$ 且 $a^n\equiv 1(\mod m),n>0$ ，则 $\delta_m(a) | n$
2. $\delta_{m}(a) | \phi(m)$
3. 若 $p,q$ 为奇素数，且 $q | (a^p-1)$，则或有 $q | (a-1)$，或有 $q=2kp+1$，其中 $k$ 为某整数
4. $2^p−1$ 的任何因子必取 $2kp+1$ 的形式。
5. $\large a^{0}, a^{1}, \dots, a^{\delta_m(a)-1} $ 构成模 m 的既约剩余系
6. $a^{x} \equiv a^{y}(\bmod m) \Leftrightarrow x \equiv y\left(\bmod \delta_{m}(a)\right)$
7. $\Large \delta_{m}\left(a^{d}\right)=\frac{\delta_{m}(a)}{\left(\delta_{m}(a), d\right)}$

求法

分解 $\Large \phi(m)=\prod p_i^{k_i}$ ， 逐渐减小 $k_i$ ，记为 $P$，直到再减一 $[a^P\ne 1(\mod m)]$ 为止

ll qPow(ll a, ll b, ll c) {  //求(a^b) % c
    ll ret = 1;
    while (b) {
        if (b & 0x1) ret = ret * a % c;
        a = a * a % c;
        b >>= 1;
    }
    return ret;
}
ll gcd(ll a,ll b){
    return b?gcd(b,a%b):a;
}
int euler(int n){
    int m=(int)sqrt(n+0.5);
    int ans=n;
    for(int i=2;i<=m;i++){
        if(n%i==0){
            ans=ans/i*(i-1);
            while(n%i==0) n/=i;
        }
    }
    if(n>1) ans=ans/n*(n-1);
    return ans;
}
pair<ll,int> factor[50];
int num_factor = 0;
void decomposition(ll n) {
    num_factor = 0;
    int m=(int)sqrt(n+0.5);
    for(int i=2;i<=m;i++){
        if(n%i==0){
            factor[num_factor]=pair<ll,int>(i,0);
            while(n%i==0){
                factor[num_factor].se++;
                n/=i;
            }
            num_factor++;
        }
    }
    if (n > 1) factor[num_factor++] = pair<ll,int>(n,1);
}
int order(int a,int m){//求a在模m意义下的阶
    if(gcd(a,m)!=1)return 0;
    int phi=euler(m);
    decomposition(phi);
    for0(i,num_factor){
        while(qPow(a,phi/factor[i].fi,m)==1&&factor[i].se){
            factor[i].se--;
            phi/=factor[i].fi;
        }
    }
    return phi;
}

原根

定义

当 $\delta_m(a)=\phi (m)$ 时称为 $a$ 为 $m$ 的原根

定理与推论

1. 若 $g$ 是 $m$ 的一个原根，则 $g, g^{2}, \cdots, g^{\phi(m)}$ 各数对模 $m$ 的最小剩余，恰是小于 $m$ 且与 $m$ 互素的 $\phi(m)$ 个正整数的一个排列。
2. $m$ 的原根(若存在)的数目为 $\large\phi(\phi(m))$
3. $m$ 有原根 $\Longleftrightarrow$ $\large m=1,2,4,p^k,2p^k$ （ $p$ 为奇素数）
4. 若 $d|(p−1)$ ,则 $x^d\equiv 1(\mod p)$ 恰有 $d$ 个解
5. 若 $p$ 为素数,$d|(p−1)$,则阶为 $d$ 的最小剩余 $(\mod p)$ 的个数为 $\phi(d)$
6. 一个数的最小原根为 $O(n^{0.25})$
7. 若 $g$ 是 $m$ 的原根，则 $g^d$ 为 $m$ 的原根的充要条件为 $(d,\phi(n))=1$

求法

分解 $\Large \phi(m)=\prod p_i^{k_i}$ ，枚举 $g$ ，直到 $\Large g^{\frac {\phi(m)} {p_i}}\ne 1(\mod m)$

ll qPow(ll a, ll b, ll c) {  //求(a^b) % c
    ll ret = 1;
    while (b) {
        if (b & 0x1) ret = ret * a % c;
        a = a * a % c;
        b >>= 1;
    }
    return ret;
}
int euler(int n){
    int m=(int)sqrt(n+0.5);
    int ans=n;
    for(int i=2;i<=m;i++){
        if(n%i==0){
            ans=ans/i*(i-1);
            while(n%i==0) n/=i;
        }
    }
    if(n>1) ans=ans/n*(n-1);
    return ans;
}
int factor[50];
int num_factor = 0;
void decomposition(ll n) {
    num_factor = 0;
    int m=(int)sqrt(n+0.5);
    for(int i=2;i<=m;i++){
        if(n%i==0){
            factor[num_factor]=i;
            while(n%i==0)n/=i;
            num_factor++;
        }
    }
    if (n > 1) factor[num_factor++] =n;
}
int root(int m){
    int phi=euler(m);
    decomposition(phi);
    int g=2;
    while(1){
        bool yes=1;
        for0(i,num_factor){
            if(qPow(g,phi/factor[i],m)==1){
                yes=0;
                break;
            }
        }
        if(yes)break;
        g++;
    }
    return g;
}

离散对数

定义

设 $g$ 是 $m$ 的一个原根，当 $g^t\equiv k(\mod m)$ 时，$in d_gk\equiv t(\mod \phi(m))$ ，此处的 $ind_gk$ 表示 $k$ 以 $g$ 为底，模 $\phi(m)$ 的离散对数值

定理与推论

1. $ind_gab\equiv ind_ga+ind_gb (\mod \phi(m))$
2. $ind_ga^n\equiv n\times ind_ga (\mod \phi(m))$

BSGS 算法

$a^x\equiv C(\% m) , gcd(a,m)=1$

$0\le x< \delta_m(a)$, 设 $x=uT-v$ , 则原方程转化为 $\large a^{uT}\equiv Ca^v(a^{-1}存在)$, 我们可以预处理出每一个 $Ca^v$ 对每个 $a^{uT}$ 找一下有没有等于 $Ca^v$

//map
ll BSGS(ll a, ll b, ll p) {  // a^x == b (mod p)
    if(a%p==0)return -1;
    if(b==1)return 0;
    map<ll, int> h;
    ll m = ceil(sqrt(p)), t = b%p, tt = qPow(a,m,p), v = 1;
    for (ll i = 0; i < m; ++i) {
        h[t] = i+1;
        t = t * a % p;
    }
    for (ll i = 1; i <= m; ++i) {
        v = v * tt % p;
        if (h[v]) return i * m - h[v]+1;
    }
    return -1;
}

//hash_table
struct Hash {
    static const int MOD = 2e6+7;
    static const int N = 2e6;
    ll head[MOD + 10], nx[N], top;
    ll hs[N], id[N];
    void init() {
        memset(head, -1, sizeof head);
        top = 0;
    }
    ll find(ll x) {
        ll k = x % MOD;
        for (int i = head[k]; i != -1; i = nx[i]) {
            if (hs[i] == x) {
                return id[i];
            }
        }
        return -1;
    }
    bool insert(ll x, ll y) {//x-key y-val
        if(find(x)!=-1)return 0;
        ll k = x % MOD;
        hs[top] = x; id[top] = y; nx[top] = head[k]; head[k] = top++;
        return 1;
    }
}hs;
ll BSGS(ll a, ll b, ll p) {  // a^x == b (mod p)
    if(a%p==0)return -1;
    if(b==1)return 0;
    hs.init();
    ll m = ceil(sqrt(p)), t = b%p , tt = qPow(a,m,p),v=1;
    for (ll i = 0; i < m; ++i) {
        hs.insert(t,i);
        t = t * a % p;
    }
    for (ll i = 1; i <= m; ++i) {
        v = v * tt % p;
        if (hs.find(v) >= 0) return i * m - hs.find(v);
    }
    return -1;
}

扩展 BSGS

初始方程为

$$a^x\equiv b\ (mod\ p)$$

设 $b’=b,p’=p$，初始化 $v= 1$，每次约分拿出一个 $a$ ，令

$$d = gcd(a, p’ )$$

如果
$d=1$，则约分过程结束，约分完成。否则的话，如果 $d$ 不能整除
$b’$，说明无法完成约分，$x$ 无解，整个算法结束。
再否则的话，更新
$b’，p’，v$ 为

$$b'=\frac {b'}d,p'=\frac {p'}d,v=v*\frac a d$$

然后继续循环进行下一轮约分。
其中 $v$ 即为 $a$ 每次约分剩下的没有被约掉的因子的乘积。约分完成后，假设约分了 $cnt$ 次，则程变为

$$v\times a^{x-cnt}\equiv b'\ (mod\ p')$$

此时
$(a，p’) = 1$，显然我们可用基本的 $BGGS$ 算法求得 $x -cnt$ 的值，算法至此结束。
具体可参考下方的实现理解。

ll exBSGS(ll a, ll b, ll p) {
  a %= p;b %= p;
  if (a == b) return 1;
  if(b==1)return 0;
  if(a==0)return -1;
  ll ret = 1;
  for (ll i = 0; i <= 50; ++i) {
    if (ret == b) return i;
    ret = (ret * a) % p;
  }  //枚举比较小的i
  ll x, y, d, v = 1, cnt = 0;
  while ((d = gcd(a, p)) != 1) {
    if (b % d) return -1;
    b /= d, p /= d;
    v = (v * (a / d)) % p;
    ++cnt;
  }  //约分直到(a, p) == 1
  if(b==1)return cnt;
  init();
  ll m = ceil(sqrt(p)), t = b , tt = qPow(a,m,p);
  for (ll i = 0; i < m; ++i) {
    Insert(i, t);
    t = t * a % p;
  }
  for (ll i = 1; i <= m; ++i) {
    v = v * tt % p;
    if (Find(v) >= 0) return i * m - Find(v)+cnt;
  }
  return -1;
}

全错位排列

$$D_n=(n-1)(D_{n-1}+D_{n-2})\\ D_n=n!M_n=n!(\frac 1 {2!}-\frac 1 {3!}+\cdots+(-1)^n\frac 1 {n!})\\$$

其他

康托展开

计算当前排列在所有由小到大全排列中的顺序
$X = A[0] \times (n-1)! + A[1] \times (n-2)! + … + A[n-1] \times 0! $

• $A[i]$ 指的是位于位置 $i$ 后面的数小于 $A[i]$ 值的个数,后面乘的就是后面还有多少个数的阶乘

• 说明 ：这个算出来的数康拖展开值，是在所有排列次序 $- 1$的值，因此 $X+1$ 即为在全排列中的次序

  //返回数组a中当下顺序的康拖映射
  int cantor(int *a,int n)
  {
  	int ans=0;
  	for(int i=0;i<n;i++)
  	{
  		int x=0;int c=1,m=1;//c记录后面的阶乘
  		for(int j=i+1;j<n;j++)
  		{
  			if(a[j]<a[i])x++;
  			m*=c;c++;
  		}
  		ans+=x*m;
  	}
  	return ans;
  }

逆康拖展开
例子 ：
在（1，2，3，4，5) 给出61可以算出起排列组合为34152
具体过程如下：
用 61 / 4! = 2余13，说明 ,说明比首位小的数有2个，所以首位为3。
用 13 / 3! = 2余1，说明 ，说明在第二位之后小于第二位的数有2个，所以第二位为4。
用 1 / 2! = 0余1，说明 ，说明在第三位之后没有小于第三位的数，所以第三位为1。
用 1 / 1! = 1余0，说明 ，说明在第二位之后小于第四位的数有1个，所以第四位为5。

static const int FAC[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};   // 阶乘
  
//康托展开逆运算
void decantor(int x, int n){
    vector<int> v;  // 存放当前可选数
    vector<int> a;  // 所求排列组合
    for(int i=1;i<=n;i++)v.push_back(i);
    for(int i=m;i>=1;i--){
        int r = x % FAC[i-1];
        int t = x / FAC[i-1];
        x = r;
        sort(v.begin(),v.end());// 从小到大排序
        a.push_back(v[t]);      // 剩余数里第t+1个数为当前位
        v.erase(v.begin()+t);   // 移除选做当前位的数
    }
}

高斯整数

高斯整数是实数和虚数部分都是整数的复数。所有高斯整数组成了一个整域，写作 $Z[i]$, 是个不可以转成有序环的欧几里德域。$Z[i] = \{a+bi | a,b \in \mathbb Z \}$

$Z[i]$ 的单位元是 $1,-1,i,-i$

对 $Z[sqrt(-k)], k \in \{1,2,3,7,11,19,43,67,163\}$ 才满足唯一分解定理

高斯素数

• $a ,b$ 中有一个是 $0$,另一个是形如 $4n+3$ 或 $-(4n+3)$ 的素数
• $a,b$ 均不为 $0$，$a^2+b^2$ 为素数

四平方和定理

每个正整数均可表示为4个整数的平方和

费马平方和定理

奇素数能表示为两个平方数之和的充分必要条件是该素数被4除余1。

如果两个整数都能表示为两个平方数之和，则它们的积也能表示为两个平方数之和。

如果一个能表示为两个平方数之和的整数被另一个能表示为两个平方数之和的素数整除，则它们的商也能表示为两个平方数之和。

如果一个能表示为两个平方数之和的整数被另一个不能表示为两个平方数之和的整数整除，则它们的商也必有一个不能表示为两个平方数之和的因子。

如果 $a$ 和 $b$ 互素，则 $a^2+b^2$ 的所有因子都能表示为两个平方数之和。

阶乘近似

$$n!\approx\sqrt{2\pi n}(\frac n e)^n$$

Ｎ以内相邻素数最大间距

$ln N \times (ln N - 2\times ln ln N)+2$

调和级数

$$f(n)\approx ln(n)+C+\frac 1 2n\\ Euler-constant:C\approx0.57721566490153286060651209$$

快速平方根倒数

$\frac 1 {\sqrt x}$

float fast_inverse_sqrt(float x) {
    float half_x = 0.5 * x;
    int i = *((int *)&x);              // 以整数方式读取X
    i = 0x5f3759df - (i >> 1);         // 神奇的步骤
    x = *((float *)&i);                // 再以浮点方式读取i
    x = x * (1.5 - (half_x * x * x));  // 牛顿迭代一遍提高精度
    return x;
}

N 次剩余

$x^a\equiv b(\% p)$

• $p$ 为素数

  设 $g$ 为 $P$ 的一原根，记 $k$ 关于 $g$ 的离散对数 $(mod p)$ 为 $ind_dk$，则有

  $$a\times ind_gx\equiv ind_gb(\mod \phi(p))$$

  设
  $t_1=ind_gb$，则有
  $g^{t_1}\equiv b (mod p)$ 。那么 $t_1$ 即与情形2—致，可由 $BSGS$ 算法求出
  $t_1$。

  设
  $t_2=ind_gx$，则原式变为
  $a\times t_2\equiv t_1 (mod \phi (p))$
  其中只有
  $t_2$ 是未知的，可由扩展欧几里德求出所有 $t_2$
  的值。
  又由
  $t_2=ind_gx$
  知 $g^{t_2}\equiv x (mod p)$ 与情形1 一致，可由快速幂求出所有 $x$ 值。
  至此解决 $P$ 为素数时的问题。

• p为合数

  首先将 $p$ 分解素因子为
  $p=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}$

  那么我们对每一个 $i$ 属于 $[1,\cdots,k]$，都有
  $x^a \equiv b (mod p_i^{e_i})$ 我们可利用上面的知识去解出这里的 $x$，不过由于模数不一定为素数的原因，不能单纯地使用上面 $p$ 为素数的方法去解，关于这个式子的解法是解决整个问题的关键，下面会给出细节，此处暂且认为已经解出 $x$ 的一
  组解 $x_{i1},x_{i2},\cdots,x_{it_i}$。

  那么上面 $k$ 个式子就会有 $k$ 组这样的解，我们从每组里取一个 $x$ 值出来，记为 $\{X_1,X_1,\cdots,X_k\}$，那么如果我们能找到一个 $z$ 满足

  $$\left\{ \begin{array}{c} z\equiv X_1\ (mod\ p_1^{e_1})\\ z\equiv X_2\ (mod\ p_2^{e_2})\\ \vdots\\ z\equiv X_k\ (mod\ p_k^{e_k})\\ \end{array} \right.$$

  那么 $z$ 必为原式的一个解。而要解这个，只需使用 $CRT$ 即可，同时我们需要一个 $dfs$ 来枚举所有 $X$ 可能的组合。
  至此解决问题。

大素数判断

/*
1. a^(p-1) ≡ 1(mod p)
2. 二次探次定理  若a^2≡1(mod p)，则a≡±1(mod p)
3. Miller-Rabin质数测试 p-1=d*2^r   a^d mod n=1或者存在a^(d*2^i) mod n=n-1 ( 0<=i<r ) 
如果n<2^64，只用选取a=2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37做测试即可
以a为基,n-1=x*2^t      a^(n-1)=1(mod n)  验证n是不是合数
*/
bool check(ll a,ll n,ll x,ll t){//一定是合数返回true,不一定返回false
	ll res = qPow(a,x,n);
	ll last = res;
	for(int i = 1; i <= t; ++i){
		res=qMulti(res,res,n);
		if(res==1 && last !=1 && last != n-1)return true;//合数
		last = res;
	}
	if(res != 1)return true;
	return false;
}
bool Miller_Rabin(ll n){
	if(n < 2) return false;
	if(n == 2) return true;
	if((n&1) == 0) return false;//偶数
	ll x = n-1,t = 0;
	while((x&1)==0){
		x >>= 1;
		++t;
	}
	for(int i = 0; i < 20; ++i){
		ll a = rand()%(n-1)+1;
		if(check(a,n,x,t))return false;//合数
	}
	return true;
}

//2
bool isPrime(ll n){
    ll d,u,t,p[]={2,3,5,7,11,13,17,19,23,29,31,37};
    for(d=n-1 ; !(d&1) ; )d>>=1;
    for(int i=0;i<12;i++){
        if(!(n%p[i]))return n==p[i];
        for(t=qPow(p[i],u=d,n);t!=n-1&&t!=1&&u!=n-1;)t=qMulti(t,t,n),u<<=1;
        if(t!=n-1&&!(u&1))retrun 0;
    }
    retrun 1;
}

Pollard-Rho 质因数分解

ll factor[100];//质因数分解结果（刚返回时是无序的）
int tol;//质因数的个数。数组小标从0开始
ll gcd(ll a,ll b){
	if(a==0)return 1;//???????
	if(a<0) return gcd(-a,b);
	while(b){
		ll t=a%b;
		a=b;
		b=t;
	}
	return a;
}
ll Pollard_rho(ll x,ll c){
	ll i=1,k=2,x0=rand()%x, y=x0;
	while(1){
		i++;
		x0=(qMulti(x0,x0,x)+c)%x;
		ll d=gcd(y-x0,x);
		if(d!=1&&d!=x) return d;
		if(y==x0) return x;
		if(i==k){
			y=x0;
			k+=k;
		}
	}
}
void findfac(ll n){ //对n进行素因子分解
	if(Miller_Rabin(n)){ 
		factor[tol++]=n;
		return;
	}
	ll p=n;
	while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
	findfac(p);
	findfac(n/p);
}