标签:
Tom放学回家的路上,看到天空中出现一个矩阵。Tom发现,如果矩阵的行、列从0开始标号,第i行第j列的数记为ai,j,那么ai,j=Cji
如果i < j,那么ai,j=0
Tom突发奇想,想求一个矩形范围内所有数的和。Tom急着回家,当然不会自己算,所以就把任务交给你了。
因为数可能很大,答案对一个质数p取模。
输入包含多组数据(大约8组)。每组数据只有一行五个非负整数,x1、y1、x2、y2、p,你要求的是∑x2i=x1∑y2j=y1ai,j模p后的值。
x1≤x2≤105,y1≤y2≤105,2≤p≤109
对于每组数据输出一行,答案模p。
0 0 1 1 7 1 1 2 2 13 1 0 2 1 2
3 4 1
对于一个 矩阵的排列组合,C(X1,Y1) +C(X1,Y1+1)+....+(C(X1,Y2)=C(X1+1,Y2)-C(X1,Y1+1);
每一列都可以这样化,
所以后面就是:C(X,Y)%P的问题,这里证明LUCAS定律
C(X,Y)=【C(X/P,Y/P)+(X%P,Y%P)】%P;
来自百度百科:
复习lucas
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <algorithm> 5 #include <queue> 6 #include <vector> 7 #include <cmath> 8 #define LL long long 9 using namespace std; 10 const int N = 100010; 11 const int mod = 1000000007; 12 LL f[N],num[N]; 13 LL p; 14 void init() 15 { 16 f[0]=1; 17 for(int i=1;i<=N-5;i++) 18 { 19 int tmp=i; 20 num[i]=num[i-1]; 21 while(tmp%p==0) 22 { 23 num[i]++; 24 tmp/=p; 25 } 26 f[i]=f[i-1]*tmp%p; 27 } 28 } 29 30 LL inv(LL a,LL n) 31 { 32 LL res=1; 33 a%=p; 34 while(n) 35 { 36 if(n&1)res=res*a%p; 37 a=a*a%p; 38 n>>=1; 39 } 40 return res; 41 } 42 43 LL calc(int a,int b) 44 { 45 if(a<b)return 0; 46 if(num[a]-num[b]-num[a-b])return 0; 47 else return f[a]*inv(f[b],p-2)%p*inv(f[a-b],p-2)%p; 48 } 49 //先提取出阶乘里面的P 后面才能求逆元 50 int main() 51 { 52 int x1,y1,x2,y2; 53 54 while(scanf("%d%d%d%d%I64d",&x1,&y1,&x2,&y2,&p)>0) 55 { 56 init(); 57 LL ans=0; 58 for(int i=y1;i<=y2;i++) 59 { 60 ans=((ans+calc(x2+1,i+1)-calc(x1,i+1))%p+p)%p; 61 } 62 printf("%I64d\n",ans); 63 } 64 }
因为p是素数,所以a!=0 关于p的逆为 a^-1=a^(p-2)%p;小费马定理
‘因为p 很小 所以要先预处理 阶乘可不可过能能mod p==0;
标签:
原文地址:http://www.cnblogs.com/forgot93/p/4643912.html