（主要是各基本知识点的例题，查缺补漏，希望以后比赛不要有本来应该掌握的知识点，却不知道的情况）
1、一元线性同余方程组的一般性解法（拓展gcd、剩余定理）
1 /*poj2891*解一次线性同余方程组/ 2 /x%a[i]==r[i]*/ 3 /*通解*/ 4 #include<iostream> 5 #include<string.h> 6 #include<stdio.h> 7 #define maxn 10000+5 8 #define LL long long 9 using namespace std; 10 LL r[maxn],m[maxn]; 11 12 void Ex_gcd(LL a,LL b,LL &d,LL &x,LL &y){ 13 if (!b){ 14 x=1;y=0;d=a; 15 }else{ 16 Ex_gcd(b,a%b,d,y,x); 17 y=y-x*(a/b); 18 } 19 } 20 LL Ex_crt(LL N){ 21 LL M=m[1],R=r[1],x,y,d; 22 for(int i=2;i<=N;i++){ 23 Ex_gcd(M,m[i],d,x,y); 24 if ((r[i]-R)%d!=0) return -1;//无解 25 x=(r[i]-R)/d*x%(m[i]/d); 26 R=R+x*M; 27 M=M/d*m[i]; 28 R%=M; 29 } 30 if (R<0) R+=M; 31 return R; 32 } 33 int n; 34 int main(){ 35 while(~scanf("%I64d",&n)){ 36 for(int i=1;i<=n;i++){ 37 scanf("%I64d%I64d",&m[i],&r[i]); 38 } 39 LL ans=Ex_crt(n); 40 printf("%I64d ",ans); 41 } 42 return 0; 43 }
1 /*hdu3579解一次线性同余方程组/ 2 /x%a[i]==r[i]*/ 3 /*通解*/ 4 /*注意，题目要求的是最小“正”整数解，如果同余方程组的解为0，则解是A[i]的最小公倍数*/ 5 #include<iostream> 6 #include<string.h> 7 #include<algorithm> 8 #include<stdlib.h> 9 #include<stdio.h> 10 #define maxn 100+5 11 #define typed long long 12 using namespace std; 13 typed r[maxn],m[maxn]; 14 15 typed gcd(typed a,typed b){ 16 return b==0?a:gcd(b,a%b); 17 } 18 void Ex_gcd(typed a,typed b,typed &d,typed &x,typed &y){ 19 if (!b){ 20 x=1;y=0;d=a; 21 }else{ 22 Ex_gcd(b,a%b,d,y,x); 23 y=y-x*(a/b); 24 } 25 } 26 typed Ex_crt(typed N){ 27 typed M=m[1],R=r[1],x,y,d; 28 for(int i=2;i<=N;i++){ 29 Ex_gcd(M,m[i],d,x,y); 30 if ((r[i]-R)%d!=0) return -1;//无解 31 x=(r[i]-R)/d*x%(m[i]/d); 32 R=R+x*M; 33 M=M/d*m[i]; 34 R%=M; 35 } 36 if (R<0) R+=M; 37 return R; 38 } 39 typed n,T; 40 typed lcm(){ 41 if (n==1) return m[1]; 42 sort(m+1,m+n+1); 43 typed GCD=gcd(m[2],m[1]); 44 for(int i=3;i<=n;i++) 45 GCD=gcd(m[i],GCD); 46 typed ans=1; 47 for(int i=1;i<=n;i++) ans=ans*m[i]/GCD; 48 return ans; 49 } 50 int main(){ 51 cin>>T; 52 for(int cas=1;cas<=T;cas++){ 53 scanf("%I64d",&n); 54 for(int i=1;i<=n;i++){ 55 scanf("%I64d",&m[i]); 56 } 57 for(int i=1;i<=n;i++){ 58 scanf("%I64d",&r[i]); 59 } 60 printf("Case %d: ",cas); 61 typed ans=Ex_crt(n); 62 if (ans==0) ans=lcm(); 63 printf("%lld ",ans); 64 } 65 return 0; 66 }
2、高次同余方程babystep算法
1 /*HDU - 2815 2 题意比较有趣， 3 但化成数学模型， 4 就是求K^D=N mod P; 5 K,P,N已知，求D的最小整数解。 6 有个小坑，N<P 才可 7 所以就是个解高次同余方程的题目 8 */ 9 #include<cmath> 10 #include<iostream> 11 #include<string.h> 12 #include<stdio.h> 13 #define maxn 65535 14 #define LL long long 15 using namespace std; 16 struct hash{ 17 int a,b,next; 18 }Hash[maxn*2]; 19 int flg[maxn+70]; 20 int top,idx; 21 void ins(int a,int b){ 22 int k=b&maxn; 23 if (flg[k]!=idx){ 24 flg[k]=idx; 25 Hash[k].next=-1; 26 Hash[k].a=a; 27 Hash[k].b=b; 28 return ; 29 } 30 while(Hash[k].next!=-1){ 31 if (Hash[k].b==b) return ; 32 k=Hash[k].next; 33 } 34 Hash[k].next=++top; 35 Hash[top].next=-1; 36 Hash[top].a=a; 37 Hash[top].b=b; 38 } 39 40 int find(int b){ 41 int k=b&maxn; 42 if (flg[k]!=idx) return -1; 43 while(k!=-1){ 44 if (Hash[k].b==b) 45 return Hash[k].a; 46 k=Hash[k].next; 47 } 48 return -1; 49 } 50 int gcd(int a,int b){ 51 return b==0?a:gcd(b,a%b); 52 } 53 int ext_gcd(int a,int b,int &x,int &y){ 54 int t,ret; 55 if (!b) { 56 x=1,y=0;return a; 57 } 58 else{ 59 ret=ext_gcd(b,a%b,x,y); 60 t=x;x=y;y=t-a/b*y; 61 return ret; 62 } 63 64 } 65 int Inval(int a,int b,int n){ 66 int x,y,e; 67 ext_gcd(a,n,x,y); 68 e=(long long) x*b%n; 69 return e<0?e+n:e; 70 } 71 int pow_mod(long long a,int b,int c){ 72 long long ret=1%c;a=a%c; 73 while(b){ 74 if (b&1) 75 ret=ret*a%c; 76 a=a*a%c; 77 b=b>>1; 78 } 79 return ret; 80 } 81 int BabyStep(int A,int B,int C){ 82 top=maxn;++idx; 83 long long buf=1%C,D=buf,K; 84 int i,d=0,tmp; 85 for(i=0;i<=100;buf=buf*A%C,i++) 86 if (buf==B) return i; 87 88 while((tmp=gcd(A,C))!=1){ 89 if (B%tmp) return -1; 90 ++d; 91 C=C/tmp; 92 B=B/tmp; 93 D=D*A/tmp%C; 94 } 95 int M=(int)ceil(sqrt((double)C)); 96 for( buf=1%C,i=0;i<=M;buf=buf*A%C,i++) ins(i,buf); 97 for( i=0,K=pow_mod((long long)A,M,C);i<=M;D=D*K%C,i++){ 98 tmp=Inval((int)D,B,C);int w; 99 if (tmp>=0 && (w=find(tmp))!=-1) 100 return i*M+w+d; 101 } 102 return -1; 103 } 104 LL K,P,N; 105 int main(){ 106 while(cin>>K>>P>>N){ 107 if (N>P){ 108 cout<<"Orz,I can’t find D!"<<endl; 109 continue; 110 } 111 N=N%P; 112 int tmp=BabyStep(K,N,P);//参数是底数，被除数，除数 113 if (tmp<0) 114 cout<<"Orz,I can’t find D!"<<endl; 115 else 116 printf("%d ",tmp); 117 } 118 return 0; 119 }
1 /*Poj 3243 2 就是求K^X=N mod P; 3 K,P,N已知，求X的最小整数解 4 解高次同余方程的题目 5 */ 6 #include<cmath> 7 #include<iostream> 8 #include<string.h> 9 #include<stdio.h> 10 #define maxn 65535 11 #define LL long long 12 using namespace std; 13 struct hash{ 14 int a,b,next; 15 }Hash[maxn*2]; 16 int flg[maxn+70]; 17 int top,idx; 18 void ins(int a,int b){ 19 int k=b&maxn; 20 if (flg[k]!=idx){ 21 flg[k]=idx; 22 Hash[k].next=-1; 23 Hash[k].a=a; 24 Hash[k].b=b; 25 return ; 26 } 27 while(Hash[k].next!=-1){ 28 if (Hash[k].b==b) return ; 29 k=Hash[k].next; 30 } 31 Hash[k].next=++top; 32 Hash[top].next=-1; 33 Hash[top].a=a; 34 Hash[top].b=b; 35 } 36 37 int find(int b){ 38 int k=b&maxn; 39 if (flg[k]!=idx) return -1; 40 while(k!=-1){ 41 if (Hash[k].b==b) 42 return Hash[k].a; 43 k=Hash[k].next; 44 } 45 return -1; 46 } 47 int gcd(int a,int b){ 48 return b==0?a:gcd(b,a%b); 49 } 50 int ext_gcd(int a,int b,int &x,int &y){ 51 int t,ret; 52 if (!b) { 53 x=1,y=0;return a; 54 } 55 else{ 56 ret=ext_gcd(b,a%b,x,y); 57 t=x;x=y;y=t-a/b*y; 58 return ret; 59 } 60 61 } 62 int Inval(int a,int b,int n){ 63 int x,y,e; 64 ext_gcd(a,n,x,y); 65 e=(long long) x*b%n; 66 return e<0?e+n:e; 67 } 68 int pow_mod(long long a,int b,int c){ 69 long long ret=1%c;a=a%c; 70 while(b){ 71 if (b&1) 72 ret=ret*a%c; 73 a=a*a%c; 74 b=b>>1; 75 } 76 return ret; 77 } 78 int BabyStep(int A,int B,int C){//A：底数，B:被除数，C除数 79 top=maxn;++idx; 80 long long buf=1%C,D=buf,K; 81 int i,d=0,tmp; 82 for(i=0;i<=100;buf=buf*A%C,i++) 83 if (buf==B) return i; 84 85 while((tmp=gcd(A,C))!=1){ 86 if (B%tmp) return -1; 87 ++d; 88 C=C/tmp; 89 B=B/tmp; 90 D=D*A/tmp%C; 91 } 92 int M=(int)ceil(sqrt((double)C)); 93 for( buf=1%C,i=0;i<=M;buf=buf*A%C,i++) ins(i,buf); 94 for( i=0,K=pow_mod((long long)A,M,C);i<=M;D=D*K%C,i++){ 95 tmp=Inval((int)D,B,C);int w; 96 if (tmp>=0 && (w=find(tmp))!=-1) 97 return i*M+w+d; 98 } 99 return -1; 100 } 101 LL K,P,N; 102 //K^X=N mod P; 103 int main(){ 104 while(cin>>K>>P>>N){ 105 if (K==0 && P==0 && N==0) break; 106 int ans=BabyStep(K,N,P); 107 if (ans==-1) puts("No Solution");else cout<<ans<<endl; 108 } 109 return 0; 110 }
3、中国剩余定理
1 /*Poj1006中国剩余定理： 2 给定n组数对，使得: 3 A[i]==X(mod M[i]) 4 5 剩余定理描述： 6 设m1,m2,…,mk是两两互素的正整数，则同余 7 方程组 8 x=b1 mod m1 9 x=b2 mod m2 10 … 11 x=bk mod mk 12 有唯一解：x=（b1T1M1+…+bkTkMk ）mod M 13 其中，M=m1×m2×…mk（即lcm） 14 Mi=M/mi 15 Ti=Mi(-1) mod mi 乘法逆元 16 17 18 pps: 19 应用中国剩余定理，可以简化一些复杂的运算。 20 如计算 47317 mod 713 21 解：因为713=2331，令x=47317， 22 则计算x mod 713等价于求解同余方程： 23 x=47317 mod 23=1317 mod 23 24 x=47317 mod 31=817 mod 31 25 即： x=6 mod 23 26 x=2 mod 31 27 根据中国剩余定理可得解x=374(mod 713) 28 因此：47317 mod 713=374(mod 713) 29 30 */ 31 /*对于这道题啊，因为数据量太小，手动计算即可， 32 M=23*28*33=21252 33 ans=(b1*924*6+b2*759*19+b3*644*2-d + M)% M;//最小正整数解 34 */ 35 #include<iostream> 36 #include<string> 37 #include<algorithm> 38 #include<cstdio> 39 #include<cstring> 40 #include<cmath> 41 #include<queue> 42 #include<map> 43 #include<stack> 44 #include<set> 45 #include<cstdlib> 46 47 using namespace std; 48 49 int d,b1,b2,b3; 50 int cas=0; 51 int main(){ 52 while(cin>>b1>>b2>>b3>>d){ 53 if (d==-1 && b1==-1 && b2==-1 && b3==-1) break; 54 cas++; 55 int M=23*28*33; 56 int ans=(b1*924*6+b2*759*19+b3*644*2-d + M)% M; 57 if (ans==0) ans=M;//正整数解！ 58 printf("Case %d: the next triple peak occurs in %d days. ",cas,ans); 59 } 60 return 0; 61 }
4、二元一次不定方程
1 /*poj2142 2 解二元一次不定方程 3 题目描述：给定质量为a和b的砝码，不限天平左右可以放置砝码，称出一个c的质量的物品就可以了 4 ax+by=d列出数学方程，x和y的正负表示放置的左右，现在要求 5 |X|+|Y|尽量小，然后是|aX|+|BY|尽量小 6 7 http://116.255.173.144/index.php?c=article&a=read&id=56659 8 因为欧几里德扩展算法只能求ax + by = gcd(a, b)；所以，先求a, b的最大公约数gcd_ab，对于ax + by = d，令a, b, d同除以gcd_ab，然后可求得(a/gcd_ab)x + (b/gcd_ab)y = 1的最小整数解，然后将求得的这组解乘以d/gcd_ab得x0, y0，就是ax + by = d一组解，由数论知识得ax + by = d的所有解为： x = x0 - bt, y = y0 + at（t为参数），那么，要使|x|+|y|最小，可令x = 0求得一组解，再令y = 0求得一组解，两组比较取较小者就好。 9 10 一般方法，拓展欧几里得解出一组特解（x0，y0） 11 x0=x0*(c/d); 12 y0=y0*(c/d); 13 令d=gcd(a,b),a0=a/d,b0=b/d; 14 则通解是x=x0+a0t，y=y0-b0t（t不所谓正负） 15 16 这道题，一般都是分析函数的单调性了（需要中学数学知识）： 17 假设a0>b0，当|x|+|y|最小时，有-b0<y<b0（可以画出函数图像） 18 尝试几组解即可 19 //http://hi.baidu.com/lydrainbowcat/item/07913682130d07dcd1f8cd83 20 */ 21 #include <iostream> 22 #include <string.h> 23 #include <stdio.h> 24 #include <stdlib.h> 25 #include <cmath> 26 #define LL long long 27 using namespace std; 28 29 LL abs(LL x){ 30 if (x<0) return -x;else return x; 31 } 32 void Ex_gcd(LL a,LL b,LL &x,LL &y,LL &d){//d==gcd(a,b) 33 if (b==0) { 34 d=a;x=1;y=0; 35 }else{ 36 Ex_gcd(b,a%b,y,x,d);y-=x*(a/b); 37 } 38 } 39 int main(){ 40 LL a,b,c,d,x0,y0; 41 while(cin>>a>>b>>c){ 42 if (a==0 && b==0 && c==0) break; 43 LL chan=0; 44 if (a<b){ 45 chan=1; 46 swap(a,b); 47 } 48 Ex_gcd(a,b,x0,y0,d); 49 x0=x0*c/d; 50 y0=y0*c/d; 51 //令y0-b0t==0解出一个特解t0，然后在t0周围测试 52 LL t0=(y0)*d/a; 53 LL a0=a/d,b0=b/d; 54 LL max_xy=999999999; 55 LL ansx,ansy; 56 LL max_xayb; 57 for(LL t=t0-10;t<=t0+10;t++){ 58 LL nx=x0+b0*t; 59 LL ny=y0-a0*t; 60 if(max_xy>abs(nx)+abs(ny)){ 61 max_xy=abs(nx)+abs(ny); 62 max_xayb=abs(nx*a)+abs(ny*b);//开始时忘记更新 63 ansx=nx;ansy=ny; 64 } 65 else if (max_xy==abs(nx)+abs(ny)){ 66 if (max_xayb>abs(nx*a)+abs(ny*b)){ 67 max_xayb=abs(nx*a)+abs(ny*b); 68 ansx=nx;ansy=ny; 69 } 70 } 71 } 72 if (chan) swap(ansx,ansy); 73 cout<<abs(ansx)<<" "<<abs(ansy)<<endl; 74 } 75 return 0; 76 }
5、n元一次不定方程
1 /*poj1091跳蚤 2 很经典的一道题，隐藏在活泼的外表下的是，n元一次不定方程。= = 3 还有，代数系统，容斥原理的应用 = =！ 4 5 题目描述: 6 1、有一只跳蚤在一个横坐标上跳着，它有1+N张卡片，从大到小，最后一张上是M。每次可选择一张，按照上面的数字，卡片可以重复使用，随便向左还是向右跳。 7 2、它的目标是跳到左边一个单位的位置 8 9 分析： 10 1、数学方程：设每张卡片的值是a1，a2.....an,M (1+n张) 11 假设相应选择了xi次（正负代表左右） 12 则 a1x1+a2x2+a3x3.....anxn+Mx(n+1)=1 13 这个方程有解的充要条件是gcd(a1,a2,a3,....,M)==1 14 2、现在问题转化成了[1，M]的范围内找出几个数（包括M），使它们的gcd等于1，问这种分组有多少？ 15 16 容斥原理 17 1、因为必须包含M这个数，所以一组中的数一定没有公因数（来自M因式分解的因子） 18 2、容斥： 19 总排列数M^N-sigm（排除含p1的分组，含p2的....）+sigm（含p1*p2的...含pi*pi+1的）-sigm（p1p2p3.....）.... 20 3、举例：含p1公因子的分组数k=M/p1 是M范围内有因子p1的数的个数.选取其中的N张（有顺序，可重复），总共K^N 21 4、表示：状态压缩，很久没用了啊 22 23 这个博客不错：http://www.cnblogs.com/zhaoguanqin/archive/2012/02/25/2368058.html 24 */ 25 #include <iostream> 26 #include <string.h> 27 #include <stdio.h> 28 #include <stdlib.h> 29 #include <cmath> 30 #define maxn 3200 31 #define LL long long 32 using namespace std; 33 bool flag[maxn]; 34 LL cnt; 35 LL prim[maxn]; 36 LL dive[25]; 37 void prim_calc(){//素数打表 38 cnt=0; 39 memset(flag,0,sizeof(flag)); 40 for(LL i=2;i<maxn;i++){ 41 if (!flag[i]) { 42 prim[cnt++]=i; 43 flag[i]=true; 44 for(LL j=2*i;j<maxn;j+=i) flag[j]=false; 45 } 46 } 47 } 48 LL count; 49 LL divide(LL M){ 50 count = 0; 51 for(LL i=0;i<cnt;i++){ 52 if (M%prim[i]==0){ 53 dive[count++]=prim[i]; 54 while(M%prim[i]==0) M=M/prim[i]; 55 if (M==1) break; 56 } 57 } 58 if (M>1) dive[count++]=M; 59 // for(LL i=0;i<count;i++) cout<<dive[i]<<" ";cout<<endl; 60 return count; 61 } 62 LL N,M; 63 int main(){ 64 prim_calc(); 65 while(cin>>N>>M){ 66 LL dig=divide(M);// 67 LL ans=pow(M,N); 68 LL tmp,flag,i,j; 69 for(i=1;i<(LL)(1<<dig);i++)//状态压缩 70 { 71 tmp=1,flag=0; 72 for(j=0;j<dig;j++) 73 if(i&((LL)(1<<j)))//判断第几个因子目前被用到 74 flag++,tmp*=dive[j]; 75 if(flag&1)//容斥原理，奇加偶减 76 ans-=pow(M/tmp,N); 77 else 78 ans+=pow(M/tmp,N); 79 } 80 cout<<ans<<endl; 81 } 82 return 0; 83 }
6、二元二次不定方程（佩尔方程）
1 /*HDU3292 2 佩尔方程（二元二次方程） 3 形如X^2-d*Y^2=1（d>1） 4 如果d为不完全数，则有解：否则，无解。 5 求解佩尔方程的两种方法 6 1、暴力法（本题） 7 2、连分数法（zoj3759：http://www.cnblogs.com/little-w/p/3588285.html） 8 求佩尔方程的解（矩阵快速幂） 9 |x[k]|=|x1 d*y1|^(k-1) |x1| 10 |y[k]|=|y1 x1 | |y1| 11 思路：暴力或者连分数法求x1+矩阵快速幂 12 */ 13 #include <iostream> 14 #include <string.h> 15 #include <stdio.h> 16 #include <cmath> 17 #define LL long long 18 #define mod 8191 19 using namespace std; 20 21 int x,y,d,k; 22 bool isquare(){ 23 int sq=sqrt(d+0.0); 24 if (((int)sq*(int)sq)==d) return true; 25 else return false; 26 } 27 void search(){ 28 y=1; 29 while(1){ 30 x=(long long)sqrt(d*y*y+1); 31 if (x*x-d*y*y==1) 32 break; 33 y++; 34 } 35 } 36 struct Mat 37 { 38 LL mat[3][3]; 39 }; 40 Mat operator*(Mat a,Mat b) 41 { 42 Mat c; 43 memset(c.mat,0,sizeof(c.mat)); 44 for (int i=0;i<2;i++) 45 for (int j=0;j<2;j++) 46 for (int k=0;k<2;k++) 47 if (a.mat[i][k] && b.mat[k][j]) 48 c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod; 49 return c; 50 } 51 Mat operator+(Mat a,Mat b) 52 { 53 Mat c; 54 memset(c.mat,0,sizeof(c.mat)); 55 for (int i=0;i<2;i++) 56 for (int j=0;j<2;j++) 57 c.mat[i][j]=(a.mat[i][j]+b.mat[i][j])%mod; 58 return c; 59 } 60 Mat E; 61 void builtE() 62 { 63 memset(E.mat,0,sizeof(E.mat)); 64 for (int i=0;i<2;i++) 65 { 66 E.mat[i][i]=1; 67 } 68 } 69 Mat operator^(Mat A,int x) 70 { 71 Mat c; 72 builtE(); 73 c=E; 74 for ( ; x; x>>=1) 75 { 76 if ( x & 1) 77 c = c * A; 78 A = A * A; 79 } 80 return c; 81 } 82 int main(){ 83 while(cin>>d>>k){ 84 if (isquare()) { 85 printf("No answers can meet such conditions "); 86 continue; 87 } 88 search(); 89 Mat A; 90 A.mat[0][0]=x,A.mat[0][1]=d*y,A.mat[1][0]=y,A.mat[1][1]=x; 91 A=A^(k-1); 92 Mat B; 93 B.mat[0][0]=x,B.mat[1][0]=y; 94 B=A*B; 95 cout<<B.mat[0][0]<<endl; 96 } 97 return 0; 98 }
7、本原毕达哥拉斯三元组定理
1 /*fzu1669 2 毕达哥拉斯三元组的应用 3 4 本原毕达哥拉斯定理： 5 x=m^2-n^2; 6 y=2mn; 7 z=m^2+n^2; 8 m>n,n奇数m偶数或者n偶数m奇数 9 注意是本原，即(x,y,z)==1,可通过计算本原加上筛选求得一定范围内的所有(x,y,z) 10 11 题意： 12 a+b+c<=L, 13 a,b是直角边，c是斜边 14 给定L ,求满足上式的三角形的个数 15 思路：通过上面的定理，枚举出m,n(需满足gcd(n,m)==1)，找到本原，筛选可得解 16 具体实现：假设m>n,则m<sqrt(L),再枚举比m小的n即可，因为L≤2000000，m<1000,单次最多不会超过500000左右 17 */ 18 #include <iostream> 19 #include <string.h> 20 #include <stdio.h> 21 #include <cmath> 22 #define LL long long 23 24 using namespace std; 25 int L; 26 int gcd(int a,int b){ 27 if (b==0) return a; 28 else return gcd(b,a%b); 29 } 30 int main(){ 31 while(cin>>L){ 32 int up=(int)(sqrt(L+0.0)); 33 int ans=0; 34 for(int m=1;m<=up;m++){ 35 for(int n=m-1;n>=1;n--){ 36 if ((m+n)%2==0 || gcd(m,n)!=1) continue; 37 int x=m*m-n*n; 38 int y=2*n*m; 39 int z=m*m+n*n; 40 ans+=L/(x+y+z); 41 } 42 } 43 cout<<ans<<endl; 44 } 45 return 0; 46 }
8、欧拉函数的应用 POJ 3696
{参考:http://blog.csdn.net/xieshimao/article/details/6689780}
9、lucas定理 HDU3944