扩展卢卡斯
本身当模数为质数时,可以直接使用卢卡斯定理:
\[\begin{pmatrix}a\\ b\end{pmatrix}
\equiv
\begin{pmatrix}\left\lfloor \frac{a}{p} \right\rfloor \\ \left\lfloor\frac{b}{p}\right\rfloor\end{pmatrix}
\times
\begin{pmatrix}a\mod p\\ b\mod p\end{pmatrix}
\pmod p
\]
然后今天做到一道模数非质数的题,需要扩卢。
套娃
学了下发现并不是很难。
首先可以将模数进行质因数拆分,然后列出同余方程组,将右边项求得便可以用 CRT/EXCRT 求解。
\[\begin{cases}
\begin{pmatrix}a\\ b\end{pmatrix}\mod {p_1}^{x_1} =s_1
\\
\begin{pmatrix}a\\ b\end{pmatrix}\mod {p_2}^{x_2} =s_2
\\
\cdots \end{cases}
\]
然后我们就需要去求每一个 \(s_i\)。
式子过几天自己推一下,防止遗忘,今天就咕咕咕。
(今天因为调模板题调了一晚上,已经把式子刻在 DNA 里了)
调了一个晚上,结果发现快速幂写成了这个鬼东西:
if(x>=pk){sum=sl[pk];for(LL i=x/pk,s=sl[pk];i;i>>=1){if(i&1)sum=sum*s%pk;s=s*s%pk;}}
wssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssb
#include <stdio.h>
#define LL long long
using namespace std;
#define sort random_shuffle
const int P=1e6+3;
const int D=8e4+3;
const int G=23;
int p,k,pk;
LL n,m;
inline LL rin()
{
LL s=0;
bool bj=false;
char c=getchar();
for(;(c>'9'||c<'0')&&c!='-';c=getchar());
if(c=='-')bj=true,c=getchar();
for(;c>='0'&&c<='9';c=getchar())s=(s<<1)+(s<<3)+(c^'0');
if(bj)s=-s;
return s;
}
struct gyq
{
int x,k;
int pk;
LL s;
}a[G];
int nam;
bool vit[P];
int d[D];
int tail;
inline void prime_()
{
int i,j;
for(i=2;i<=p;i++)
{
if(!vit[i])d[++tail]=i;
for(j=1;j<=tail;j++)
{
int now=i*d[j];
if(now>p)break;
vit[now]=true;
if(i%d[j]==0)break;
}
}
for(i=1;i<=tail;i++)
{
if(p%d[i]==0)
{
a[++nam].x=d[i];
for(;p%d[i]==0;p/=d[i])a[nam].k++;
if(p==1)break;
}
}
return;
}
LL sl[P];
LL cuts=0;
inline LL f(LL x)
{
if(x<p)return sl[x];
cuts+=x/p;
LL sum=1;
if(x>=pk){sum=1;for(LL i=x/pk,s=sl[pk];i;i>>=1){if(i&1)sum=sum*s%pk;s=s*s%pk;}}
for(LL i=(x/pk)*pk+1;i<=x;i++)if(i%p)sum=sum*(i%pk)%pk;
sum=sum*f(x/p)%pk;
return sum;
}
inline void Ex_Gcd(LL a,LL b,LL &x,LL &y){(!b)?(x=y=1):(Ex_Gcd(b,a%b,y,x),y-=x*(a/b));return;}
inline void Ex_Lucus()
{
sl[0]=sl[1]=1;
for(int i=1;i<=nam;i++)
{
p=a[i].x;k=a[i].k;pk=1;
for(int j=k;j;j--)pk*=p;a[i].pk=pk;
for(int j=2;j<=pk;j++){sl[j]=sl[j-1];if(j%p)sl[j]=sl[j]*j%pk;}
LL cutt=0;
LL sr_1,sr_2;
cuts=0;a[i].s=f(n);cutt+=cuts;
cuts=0;sr_1=f(m);cutt-=cuts;
cuts=0;sr_2=f(n-m);cutt-=cuts;
Ex_Gcd(sr_1,pk,sr_1,cuts);(sr_1=sr_1%pk+pk)%pk;
Ex_Gcd(sr_2,pk,sr_2,cuts);(sr_2=sr_2%pk+pk)%pk;
a[i].s=(a[i].s*sr_1%pk)*sr_2%pk;
a[i].s=(a[i].s%pk+pk)%pk;
if(cutt>=k)a[i].s=0;
else for(LL s=p;cutt;cutt>>=1){if(cutt&1)a[i].s=a[i].s*s%pk;s=s*s%pk;}
}
return;
}
inline void Ex_Crt()
{
LL ans=a[1].s;
LL sum=a[1].pk;
for(int i=2;i<=nam;i++)
{
Ex_Gcd(sum,a[i].pk,m,n);m=m*sum;
sum*=a[i].pk;
ans+=m*(a[i].s-ans)%sum;
ans=(ans%sum+sum)%sum;
}
printf("%lld\n",ans);
}
int main()
{
int i,j;
n=rin();m=rin();p=rin();
prime_();
Ex_Lucus();
Ex_Crt();
return 0;
}
$$\texttt{Dirty Deeds Done Dirt Cheap}$$