exCRT&exLucas

前置：

数学公式整理

exCRT（扩展中国剩余定理）

P4777 【模板】扩展中国剩余定理（EXCRT）

中国在剩余定理不能处理模数不互质的情况，这个看似小的问题，却无法进行优化中国剩余定理来克服这个问题。

因为中国剩余定理的核心就是通过每组同余式的通解联立。那么只能通过别的方法，来实现这个目的。

我们用的方法是每回将两个同余式合并，找到通解，最后将同余式组合并成一个同余式后，就可以直接找到最小解了。

我们现在处理两个同余式合并：

$$\begin{cases} x\equiv b_1(\mod a_1)\\ x\equiv b_2(\mod a_2) \end{cases}$$

如果这两个同余式不冲突，那么一定可以找到 $k_1,k_2$ 使下面的等式成立：

$$x=k_1a_1+b_1=k_2a_2+b_2$$

通过移项得到：

$$k_1a_1-k_2a_2=b_2-b_1$$

假设 $\gcd(a_1,a_2)\not =1$，那么就有：

$$k_1\frac{a_1}{\gcd(a_1,a_2)}-k_2\frac{a_2}{\gcd(a_1,a_2)}=\frac{b_2-b_1}{\gcd(a_1,a_2)}$$

其中，如果 $\gcd(a_1,a_2)\nmid (b_2-b_1)$ 那么就可以直接返回无解，因为 $k_1,k_2$ 可以凑成同余式的最小的解就是 $\gcd(a_1,a_2)$，如果不能整除，那不可能有解。

假设 $q_1=\frac{a_1}{\gcd(a_1,a2)}$，$q_2=\frac{a_2}{\gcd(a_1,a_2)}$，$r=\frac{b_2-b_1}{\gcd(a_1,a_2)}$

$$k_1q_1-k_2q_2=r$$

我们可以用扩展欧几里得（exgcd）求出 $k_i,k_j$。

$$k_iq_1+k_jq_2=1$$ $$\begin{cases} k_1=rk_i\\ k_2=-rk_j \end{cases}$$ $$k_1\equiv k_i^r\equiv \frac{b_2}{\gcd(a_1,a_2)}(\mod \frac{a_2}{\gcd(a_1,a_2)})$$

带入：

$$x=k_1a_1+b_1=rk_ia_1+b_1=\frac{b2-b1}{\gcd(a_1,a_2)}k_ia_1+b_1$$

那么我们就有了一个 $x$ 同时满足两个同余式。

定理：

如果有一个特解 $x_k$，那么

$$\begin{cases} x\equiv b_1(\mod a_1)\\ x\equiv b_2(\mod a_2) \end{cases}$$

的通解：

$$x_k+k\times lcm(m_1,m_2)$$

也就是：

$$x\equiv x_k(\mod lcm(a_1,a_2))$$

所有同余式合并就剩一个后，此时的 $x_k$ 就是最小解，输出即可。

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
using namespace std;
long long r_r(){//快读 
  long long x=0,f=1;
  char c=getchar();
  while(!isdigit(c)){
    if(c=='-')f=-1;
    c=getchar();
  }
  while(isdigit(c)){
    x=(x<<1)+(x<<3)+(c^48);
    c=getchar();
  }
  return x*f; 
}
const int o_o=1e5+10;
long long a_i[o_o],b_i[o_o];
int n;
long long k_m(long long a,long long b,long long p){//快速幂 
  long long a_s=0;
  while(b){
    if(b&1)a_s=(a_s+a)%p;
    a=(a+a)%p;
    b>>=1;
  }
  return a_s;
}
long long g_d(long long a,long long b,long long &x,long long &y){//扩展欧几里得 
  if(!b){
    x=1;y=0;
    return a;
  }
  long long g_g=g_d(b,a%b,y,x);
  y-=(a/b)*x;
  return g_g;
}
long long g_a(){
  long long x,y;
  long long a_1=a_i[0],b_1=b_i[0];//“初始化” 
  for(int i=1;i<n;i++){//挨个合并 
    long long a_a=a_i[i];//新的要合并的模数 
    long long k_k=(b_i[i]-b_1%a_a+a_a)%a_a;//移项后的值 
    long long g_g=g_d(a_1,a_a,x,y);//最大公因数 
    //x,y 即是 k_i,k_j 
    long long a_2=a_a/g_g;//约分 
    if(k_k%g_g!=0)return -1;//右边不能整除，无解 
    x=k_m(x,k_k/g_g,a_2);//求出 k_1 
    b_1+=x*a_1;//特解 
    a_1*=a_2;//最小公倍数 
    b_1=(b_1%a_1+a_1)%a_1;//通解 
  }
  return (b_1%a_1+a_1)%a_1;//最小解 
}
int main(){
  n=r_r();
  for(int i=0;i<n;i++)a_i[i]=r_r(),b_i[i]=r_r();//读入每组方程 
  printf("%lld",g_a());//输出结果 
  return 0;
} 

exLucas（扩展卢卡斯定理）

P4720 【模板】扩展卢卡斯定理/exLucas

题目要求：

$$C_n^m\mod p$$

此时的 $p$ 不保证是质数，所以假设可以拆分为：$p=p_1^{a_1}p_2^{a_2}p_3^{a_3}…$

$$\begin{cases} C_n^m\mod p_1^{a_1}\\ C_n^m\mod p_1^{a_2}\\ … \end{cases}$$

最后可以靠中国剩余定理统计答案。

现在先处理其中一个：

$$C_n^m\mod p^k$$ $$=\frac{n!}{m!(n-m)!}\mod p^k$$

此时的 $m!,(n-m)!$ 是不能求逆元的，因为不能保证互质。那么我们设

$$\gcd(n!,p^k)=p^a$$ $$\gcd(m!,p^k)=p^b$$ $$\gcd((n-m)!,p^k)=p^c$$

那么就化成：

$$\frac{\frac{n!}{p^a}}{\frac{m!}{p^b}\frac{(n-m)!}{p^c}}p^{a-b-c}\mod p^k$$

这样我们就可以将分母化成逆元的形式取模了。

但是，此时仍然没有解决问题，数据范围仍然可以卡爆暴力，所以要化简。

$$\frac{n!}{p^a}\mod p^k$$

我们先看其中的一部分。

其中：$n!=1\times 2\times 3\times …\times n$

既然要约分化简，那么就要将 $n!$ 中关于 $p$ 的倍数的数找出来：

$$=(p\times 2p\times …)(1\times 2\times …)$$

再 $1$ 到 $n$ 中有 $\lfloor \frac{n}{p}\rfloor$ 个 $p$ 的倍数的数。

那么就可以写成：

$$=p^{\lfloor \frac{n}{p}\rfloor}(1\times 2\times …\times \lfloor \frac{n}{p}\rfloor)(1\times 2\times …)$$

注意此时最右边的括号中是上面剩下来的，也就是和 $p$ 互质的数。

右边要处理，它们在 $\mod p$ 下是有循环节的。例如：

$1,2,3,4,5,6,7,8,9,10,11$

在 $\mod 5$ 后：

$1,2,3,4,0,1,2,3,4,0,1$

将例子中的逗号换成乘号再将所有的 $0$ 去掉，就是上面说的循环节的情况。（去掉 $0$ 是因为已经全部提出来了。）

注意，最后可能还剩几个不能凑成一个完成的循环节的，要单独处理，也叫做余项。

而可以凑成循环节的部分就可以计算一个循环节，用快速幂的方式计算总共的值，同样有 $\lfloor \frac{n}{p}\rfloor$ 个循环节。

那么式子就变成了：

$$=p^{\lfloor{\frac{n}{p}}\rfloor}{\lfloor \frac{n}{p}\rfloor}!{(\prod_{i=1,i\not\equiv 0(\mod p)}^{p^k}i)^{\lfloor {\frac{n}{p}}\rfloor}}{ (\prod_{i=p^k\lfloor{\frac{n}{p^k}}\rfloor,i\not\equiv 0(\mod p)}^ni)}$$

最后那部分就是余项，倒数第二部分就是循环节部分。

此时我们已经做出了本质上的优化。

但是我们只是讨论了其中的一部分，现在找一些规律，让其他部分也可以直接用到。

现在定义：$f(n)=\frac{n!}{p^a}$

那么就有：

$$f(n)=\frac{1}{p^a}p^{\lfloor{\frac{n}{p}}\rfloor}{\lfloor \frac{n}{p}\rfloor}!{(\prod_{i=1,i\not\equiv 0(\mod p)}^{p^k}i)^{\lfloor {\frac{n}{p}}\rfloor}}{ (\prod_{i=p^k\lfloor{\frac{n}{p^k}}\rfloor,i\not\equiv 0(\mod p)}^ni)}$$

最前面的 $p^{\lfloor \frac{n}{p}\rfloor}$ 要除掉，而且 $(\lfloor \frac{n}{p}\rfloor!)$ 还包含 $p$，而且此时只有这部分存在 $p$ 的倍数，那么我们的 $p^a$ 只有在这一部分能消掉，而组合后又是一个新的部分推到，也就是：

$$\frac{1}{p^a}p^{\lfloor \frac{n}{p}\rfloor}\lfloor \frac{n}{p}\rfloor!=f(\lfloor \frac{n}{p}\rfloor)$$

那么，我们就有了推导式：

$$f(n)=f(\lfloor\frac{n}{p}\rfloor){(\prod_{i=1,i\not\equiv 0(\mod p)}^{p^k}i)^{\lfloor {\frac{n}{p}}\rfloor}}{ (\prod_{i=p^k\lfloor{\frac{n}{p^k}}\rfloor,i\not\equiv 0(\mod p)}^ni)}$$

此时我们就可以通过递推来求值了，将这个函数带入到原分式中：

$$=\frac{f(n)}{f(m)f(n-m)}p^{a-b-c}\mod p^k$$

好看了不少。

但是还有一个问题，$p^{a-b-c}$ 到底是多少？

我们定义的 $p^a$ 代表的是 $\gcd(n!,p^k)$ 其他几个意义相同。

而在我们处理部分分式的时候，我们在前面提出 $p$ 的一些倍数，也就是：

$$p^{\lfloor \frac{n}{p}\rfloor}\lfloor \frac{n}{p}\rfloor!$$

这一部分。

我们可以直接获得的有 $p^{\lfloor \frac{n}{p}\rfloor}$，而剩下的要继续递推 $f(\lfloor \frac{n}{p}\rfloor)$ 才能得到。

那么我们定义一个新函数 $g(n)$。

$$g(n)=\lfloor \frac{n}{p}\rfloor+g(\lfloor \frac{n}{p}\rfloor)$$

也是通过递推的方式求解。

那么最后的答案就是：

$$=\frac{f(n)}{f(m)f(n-m)}p^{g(n)-g(m)-g(n-m)}\mod p^k$$

到这里仍然没有结束！

因为这里的 $p^k$ 只是原 $p$ 中多因子的其中一个，我们还要将所有的因子处理一遍，最后还要用中国剩余定理求出最后的答案即可。

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<set>
using namespace std;
long long r_r(){//快读 
  long long x=0,f=1;
  char c=getchar();
  while(!isdigit(c)){
    if(c=='-')f=-1;
    c=getchar();
  }
  while(isdigit(c)){
    x=(x<<1)+(x<<3)+(c^48);
    c=getchar();
  }
  return x*f; 
}
const int o_o=110;
void e_x(long long a,long long b,long long&x,long long&y){//扩展欧几里得 
    if(!b){
    x=1,y=0;
    return;
  }
    e_x(b,a%b,y,x);
    y-=a/b*x;
}
long long t_j(long long a,long long p){//找特解 
    long long x,y;
    e_x(a,p,x,y);
    return(x+p)%p;
}
long long k_m(long long a,long long b,long long p){//快速幂 
  long long r_s=1;
  a%=p;
  while(b){
    if(b&1)r_s=(r_s*a)%p;
    a=(a*a)%p;
    b>>=1;
  }
  return r_s;
}
long long f_i(long long n,long long k,long long p_k){//定义的 f 函数 
    if(n==0)return 1;
    long long x_j=1;//循环节 
    long long y_v=1;//余项 
    for(long long i=1;i<=p_k;i++)
        if(i%k)x_j=x_j*i%p_k;//1 到 p 中所有无 p 因子数的乘积 
    x_j=k_m(x_j,n/p_k,p_k);//快速幂 
    for(long long i=p_k*(n/p_k);i<=n;i++)
        if(i%k)y_v=y_v*(i%p_k)%p_k;//余项的乘积 
    return f_i(n/k,k,p_k)*x_j%p_k*y_v%p_k;//累计 
}
long long g_i(long long n,long long k){//定义的 g 函数 
    if(n<k)return 0;
    return g_i(n/k,k)+(n/k);//累计 x 
}
long long c_k(long long n,long long m,long long k,long long p_k){
    long long f_z=f_i(n,k,p_k);//分子 
  long long f_1=t_j(f_i(m,k,p_k),p_k);//分母左部分 
  long long f_2=t_j(f_i(n-m,k,p_k),p_k);//分母右部分 
    long long m_i=k_m(k,g_i(n,k)-g_i(m,k)-g_i(n-m,k),p_k);//右边部分 
    return f_z*f_1%p_k*f_2%p_k*m_i%p_k;//累计结果 
}
long long a_i[o_o],b_i[o_o];
long long e_l(long long n,long long m,long long m_d){
    long long k_k=m_d,x_p=0;
    long long t_p=2;
  while(t_p*t_p<=m_d){//分解质因数 
      if(!(k_k%t_p)){//能除尽 
            long long p_k=1;//初始化 
            while(!(k_k%t_p)){
                p_k*=t_p;//累计 
        k_k/=t_p;//化简 
            }
            a_i[++x_p]=p_k;//记录 
      b_i[x_p]=c_k(n,m,t_p,p_k);//套公式 
        }
    t_p++;
  }
    if(k_k!=1){//还有一个大质数 
        a_i[++x_p]=k_k;//记录 
    b_i[x_p]=c_k(n,m,k_k,k_k);//套公式 
    }
    
    //中国剩余定理，统计答案 
    long long a_s=0;
    for(long long i=1;i<=x_p;i++){//枚举每个因数 
        long long m_i=m_d/a_i[i];//除掉这个因子 
    long long t_i=t_j(m_i,a_i[i]);//找特解 
        a_s=(a_s+b_i[i]*m_i%m_d*t_i%m_d)%m_d;//套公式累计 
    }
    return a_s;
}
int main(){
    long long n=r_r(),m=r_r(),m_d=r_r();
    printf("%lld\n",e_l(n,m,m_d));//处理答案 
    return 0;
}