[理学]ACM必须的算法数论.doc

ACM必须的算法欧几里德#include<iostream>using namespace std;int hcf(int a,int b) int r=0; while(b!=0) r=a%b; a=b; b=r; return(a); lcd(int u,int v,int h) /u=a，v=b，h为最小公约数=hcf(a,b); return(u*v/h);int main()int a,b,x,y;cin>>a>>b;x=hcf(a,b);y=lcd(a,b,x);cout<<x<<" "<<y<<endl;return 0;扩展欧几里德#include <iostream>using namespace std;_int64 ext_euclid(_int64 a,_int64 b, _int64 &x, _int64 &y) int t;_int64 d; if (b=0) x=1;y=0;return a; d=ext_euclid(b,a %b,x,y); t=x; x=y; y=t-a/b*y; return d;void modular_equation(_int64 a,_int64 b,_int64 c)/ax = b(mod n) _int64 d; _int64 x,y; d = ext_euclid(a,b,x,y); if ( c % d != 0 )printf("No answern"); else x = (x * c/d) % b ;/ 第一次求出的x ; _int64 t = b/d;x = (x%t + t)%t;printf("%I64dn",x);/最小的正数的值for (int i=0;i<d;i+)printf("The %dth answer is : %ldn",i+1,(x+i*(b/d)%b); /所有的正数值 /*函数返回值为gcd(a,b),并顺带解出ax+by=gcd(a,b)的一个解x,y, 对于不定方程ax+by=c的通解为： x=x*c/d+b/d*t y=y*c/d+a/d*t 当c%gcd(a,b)!=0时,不定方程无解.*/ 中国剩余定理#include <iostream>using namespace std; int ext_euclid(int a,int b,int &x,int &y) /求gcd(a,b)=ax+by int t,d; if (b=0) x=1;y=0;return a; d=ext_euclid(b,a %b,x,y); t=x; x=y; y=t-a/b*y; return d;int China(int W,int B,int k) /W为按多少排列，B为剩余个数 W>B K为组数 int i; int d,x,y,a=0,m,n=1; for (i = 0;i<k;i+) n *= Wi; for (i=0;i<k;i+)m=n/Wi; d=ext_euclid(Wi,m,x,y); a=(a+y*m*Bi)%n; if (a>0)return a; elsereturn(a+n);int main()int B100,W100; 求int k ; a = 2( mod 5 )cin >> k ; a = 3( mod 13)for(int i = 0 ; i < k ;i+) 的解 2cin >> Wi; 5 2cin >> Bi; 13 3 输出 42cout<<China(W,B,k)<<endl;return 0;欧拉函数求小于n的所有欧拉数#include <iostream>using namespace std;int phi1000; /数组中储存每个数的欧拉数void genPhi(int n)/求出比n小的每一个数的欧拉数(n-1的)int i, j, pNum = 0 ;memset(phi, 0, sizeof(phi) ;phi1 = 1 ;for(i = 2; i < n; i+)if(!phii)for(j = i; j < n; j += i)if(!phij)phij = j;phij = phij / i * (i - 1);求n的欧拉数int eular(int n)int ret=1,i;for (i=2;i*i<=n;i+)if (n%i=0)n/=i,ret*=i-1;while (n%i=0)n/=i,ret*=i;if (n>1)ret*=n-1;return ret; /n的欧拉数行列式计算#include <iostream>using namespace std;int js(int s100100,int n) int z,j,k,r,total=0; int b100100; /*bNN用于存放，在矩阵sNN中元素s0的余子式*/ if(n>2) for(z=0;z<n;z+) for(j=0;j<n-1;j+) for(k=0;k<n-1;k+) if(k>=z) bjk=sj+1k+1; else bjk=sj+1k; if(z%2=0) r=s0z*js(b,n-1); /*递归调用*/ else r=(-1)*s0z*js(b,n-1); total=total+r; else if(n=2)total=s00*s11-s01*s10; return total; 排列long A(long n,long m) /n>m long a=1; while(m!=0) a*=n;n-;m-; return a; 组合long C(long n,long m) /n>m long i,c=1; i=m; while(i!=0) c*=n;n-;i-;while(m!=0) c/=m;m-; return c; 大数乘大数#include <iostream>#include <string>using namespace std;char a1000,b1000,s10000;void mult(char a,char b,char s) /a被乘数，b乘数，s为积 int i,j,k=0,alen,blen,sum=0,res6565=0,flag=0; char result65; alen=strlen(a);blen=strlen(b); for (i=0;i<alen;i+) for (j=0;j<blen;j+) resij=(ai-'0')*(bj-'0'); for (i=alen-1;i>=0;i-) for (j=blen-1;j>=0;j-) sum=sum+resi+blen-j-1j; resultk=sum%10; k=k+1; sum=sum/10; for (i=blen-2;i>=0;i-) for (j=0;j<=i;j+) sum=sum+resi-jj; resultk=sum%10; k=k+1; sum=sum/10; if (sum!=0) resultk=sum;k=k+1; for (i=0;i<k;i+) resulti+='0' for (i=k-1;i>=0;i-) si=resultk-1-i; sk='0' while(1) if (strlen(s)!=strlen(a)&&s0='0') strcpy(s,s+1); else break; int main()cin>>a>>b;mult(a,b,s);cout<<s<<endl;return 0;大数乘小数#include <iostream>using namespace std;char a100,t1000;void mult(char c,int m,char t) / c为大数,m<=10,t为积 int i,l,k,flag,add=0; char s100; l=strlen(c); for (i=0;i<l;i+) sl-i-1=ci-'0' for (i=0;i<l;i+)k=si*m+add;if (k>=10) si=k%10;add=k/10;flag=1; else si=k;flag=0;add=0; if (flag) l=i+1;si=add; else l=i; for (i=0;i<l;i+) tl-1-i=si+'0' tl='0'int main()int i;cin>>a>>i;mult(a,i,t);cout<<t<<endl;return 0;大数加法#include <iostream>#include <string>using namespace std;char a1000,b1000,s10000;void add(char a,char b,char s)/a被加数，b加数，s和 int i,j,k,up,x,y,z,l; char *c; if (strlen(a)>strlen(b) l=strlen(a)+2; else l=strlen(b)+2; c=(char *) malloc(l*sizeof(char); i=strlen(a)-1; j=strlen(b)-1; k=0;up=0; while(i>=0|j>=0)if(i<0) x='0' else x=ai;if(j<0) y='0' else y=bj;z=x-'0'+y-'0'if(up) z+=1;if(z>9) up=1;z%=10; else up=0;ck+=z+'0'i-;j-; if(up) ck+='1' i=0; ck='0' for(k-=1;k>=0;k-) si+=ck; si='0' int main()cin>>a>>b;add(a,b,s);cout<<s<<endl;return 0;大数减法#include <iostream>using namespace std;char a1000,b1000,s10000;void sub(char a,char b,char s) int i,l2,l1,k; l2=strlen(b);l1=strlen(a); sl1='0'l1-; for (i=l2-1;i>=0;i-,l1-) if (al1-bi>=0) sl1=al1-bi+'0' else sl1=10+al1-bi+'0' al1-1=bl1-1-1; k=l1; while(ak<0) ak+=10;ak-1-=1;k-; while(l1>=0) sl1=al1;l1-;loop: if (s0='0') l1=strlen(a); for (i=0;i<l1-1;i+) si=si+1; sl1-1='0' goto loop; if (strlen(s)=0) s0='0's1='0' int main()cin>>a>>b;sub(a,b,s);cout<<s<<endl;return 0;大数阶乘#include <iostream>#include <cmath>using namespace std;long a10000;int factorial(int n) /n为所求阶乘的n!的nint i,j,c,m=0,w; a0=1; for(i=1;i<=n;i+) c=0; for(j=0;j<=m;j+) aj=aj*i+c; c=aj/10000; aj=aj%10000; if(c>0) m+;am=c; w=m*4+log10(am)+1;printf("%ld",am); / 输出for(i=m-1;i>=0;i-) /printf("%4.4ld",ai);/printf("n");return w; /返回值为阶乘的位数储存方法很巧，每一个ai中存四位，不足四位在前加0补齐大数求余int mod(int B) /A为大数，B为小数int i = 0,r = 0;while( Ai != '0' )r=(r*10+Ai+-'0')%B; return r ; /r为余数高精度任意进制转换#include <iostream>#include <string>using namespace std;char s1000,s21000; / s:原进制数字，用字符串表示，s2:转换结果，用字符串表示long d1,d2; / d1:原进制数，d2:需要转换到的进制数void conversion(char s,char s2,long d1,long d2) long i,j,t,num; char c; num=0; for (i=0;si!='0'i+) if (si<='9'&&si>='0') t=si-'0' else t=si-'A'+10; num=num*d1+t; i=0; while(1) t=num%d2; if (t<=9) s2i=t+'0' else s2i=t+'A'-10; num/=d2; if (num=0) break; i+; for (j=0;j<=i/2;j+)c=s2j;s2j=s2i-j;s2i-j=c; s2i+1='0'int main()while (1)cin>>s>>d1>>d2;conversion(s,s2,d1,d2);cout<<s2<<endl;return 0;判断一个数是否素数#include <iostream>/基本方法，n为所求数，返回1位素数，0为合数#include <cmath>using namespace std;int comp(int n)int i,flag=1; for (i=2;i<=sqrt(n);i+) if (n%i=0) flag=0;break; if (flag=1) return 1; else return 0;素数表int prime(int a,int n) /小于n的素数 int i,j,k,x,num,*b; n+; n/=2; b=(int *)malloc(sizeof(int)*(n+1)*2); a0=2;a1=3;num=2; for(i=1;i<=2*n;i+) bi=0; for(i=3;i<=n;i+=3) for(j=0;j<2;j+) x=2*(i+j)-1; while(bx=0) anum+=x; for(k=x;k<=2*n;k+=x) bk=1; return num; /小于n的素数的个数bool flag1000000;void prime(int M) /01表int i , j;int sq = sqrt(double(M);for(i = 0 ;i < M ;i +)flagi = true;flag1 = false;flag0 = false;for(i = 2 ;i <= sq ;i+)if(flagi)for(j = i*i ;j < M ;j += i)flagj = false; Miller_Rabin随机素数测试算法 说明：这种算法可以快速地测试一个数是否 满足素数的必要条件，但不是充分条件。不 过也可以用它来测试素数，出错概率很小， 对于任意奇数n>2和正整数s,该算法出错概率 至多为2(-s)，因此，增大s可以减小出错概 率，一般取s=50就足够了。#include<iostream>#include <cmath>using namespace std;int Witness(int a, int n) int i, d = 1, x; for (i = ceil( log( (float) n - 1 ) / log(2.0) ) - 1; i >= 0; i-) x = d; d = (d * d) % n; if ( (d = 1) && (x != 1) && (x != n-1) ) return 1; if ( ( (n - 1) & ( 1<<i ) ) >0 )d = (d * a) % n; return (d = 1 ? 0 : 1); int Miller_Rabin(int n, int s) int j, a; for (j = 0; j < s; j+) a = rand() * (n - 2) / RAND_MAX + 1; if (Witness(a, n) return 0; return 1; int main()int x;cin>>x;cout<<Miller_Rabin(x , 50)<<endl;return 0;KMP#include<iostream> #include<string> using namespace std; char t1010,p100; /t为大字符串，p为小字符串 int next100;void sn()int j,k; next0=-1; j=1; do k=nextj-1; while(k!=-1 && pk!=pj) k=nextk; nextj=k+1; j+=1; while(pj!='0'); int match(int x) /x是大字符串下标,从头开始为0，j为小字符串下标 int k=x,j=0; if(t0='0')return 0; while(tk!='0') while(j!=-1 && pj!=tk) j=nextj; k+=1;j+=1; if(pj='0')return k-j; /返回p所在的下标 return -1; /搜索失败返回-1 int main()int x=0;sn();int r=match( x );cout<<r<<endl;（一）巴什博奕（Bash Game）：只有一堆n个物品，两个人轮流从这堆物品中取物，规定每次至少取一个，最多取m个。最后取光者得胜。     显然，如果n=m+1，那么由于一次最多只能取m个，所以，无论先取者拿走多少个，后取者都能够一次拿走剩余的物品，后者取胜。因此我们发现了如何取胜的法则：如果n=（m+1）r+s，（r为任意自然数，sm),那么先取者要拿走s个物品，如果后取者拿走k（m)个，那么先取者再拿走m+1-k个，结果剩下（m+1）（r-1）个，以后保持这样的取法，那么先取者肯定获胜。总之，要保持给对手留下（m+1）的倍数，就能最后获胜。     这个游戏还可以有一种变相的玩法：两个人轮流报数，每次至少报一个，最多报十个，谁能报到100者胜。（二）威佐夫博奕（Wythoff Game）：有两堆各若干个物品，两个人轮流从某一堆或同时从两堆中取同样多的物品，规定每次至少取一个，多者不限，最后取光者得胜。     这种情况下是颇为复杂的。我们用（ak，bk）（ak bk ,k=0，1，2，.,n)表示两堆物品的数量并称其为局势，如果甲面对（0，0），那么甲已经输了，这种局势我们称为奇异局势。前几个奇异局势是：（0，0）、（1，2）、（3，5）、（4，7）、（6，10）、（8，13）、（9，15）、（11，18）、（12，20）。     可以看出,a0=b0=0,ak是未在前面出现过的最小自然数,而 bk= ak + k，奇异局势有如下三条性质：     1。任何自然数都包含在一个且仅有一个奇异局势中。     由于ak是未在前面出现过的最小自然数，所以有ak > ak-1 ，而 bk= ak + k > ak-1 + k-1 = bk-1 > ak-1 。所以性质1。成立。     2。任意操作都可将奇异局势变为非奇异局势。     事实上，若只改变奇异局势（ak，bk）的某一个分量，那么另一个分量不可能在其他奇异局势中，所以必然是非奇异局势。如果使（ak，bk）的两个分量同时减少，则由于其差不变，且不可能是其他奇异局势的差，因此也是非奇异局势。     3。采用适当的方法，可以将非奇异局势变为奇异局势。     假设面对的局势是（a,b），若 b = a，则同时从两堆中取走 a 个物体，就变为了奇异局势（0，0）；如果a = ak ，b > bk，那么，取走b   - bk个物体，即变为奇异局势；如果 a = ak ，   b < bk ,则同时从两堆中拿走 ak - ab - ak个物体,变为奇异局势（ ab - ak , ab - ak+ b - ak）；如果a > ak ，b= ak + k,则从第一堆中拿走多余的数量a - ak 即可；如果a < ak ，b= ak + k,分两种情况，第一种，a=aj （j < k）,从第二堆里面拿走 b - bj 即可；第二种，a=bj （j < k）,从第二堆里面拿走 b - aj 即可。     从如上性质可知，两个人如果都采用正确操作，那么面对非奇异局势，先拿者必胜；反之，则后拿者取胜。     那么任给一个局势（a，b），怎样判断它是不是奇异局势呢？我们有如下公式：     ak =k（1+5）/2，bk= ak + k   （k=0，1，2，.,n 方括号表示取整函数)奇妙的是其中出现了黄金分割数（1+5）/2 = 1。618.,因此,由ak，bk组成的矩形近似为黄金矩形，由于2/（1+5）=（5-1）/2，可以先求出j=a（5-1）/2，若a=j（1+5）/2，那么a = aj，bj = aj + j，若不等于，那么a = aj+1，bj+1 = aj+1+ j + 1，若都不是，那么就不是奇异局势。然后再按照上述法则进行，一定会遇到奇异局势。（三）尼姆博奕（Nimm Game）：有三堆各若干个物品，两个人轮流从某一堆取任意多的物品，规定每次至少取一个，多者不限，最后取光者得胜。     这种情况最有意思，它与二进制有密切关系，我们用（a，b，c）表示某种局势，首先（0，0，0）显然是奇异局势，无论谁面对奇异局势，都必然失败。第二种奇异局势是（0，n，n），只要与对手拿走一样多的物品，最后都将导致（0，0，0）。仔细分析一下，（1，2，3）也是奇异局势，无论对手如何拿，接下来都可以变为（0，n，n）的情形。     计算机算法里面有一种叫做按位模2加，也叫做异或的运算，我们用符号（+）表示这种运算。这种运算和一般加法不同的一点是1+1=0。先看（1，2，3）的按位模2加的结果：1 =二进制012 =二进制103 =二进制11 （+）0 =二进制00 （注意不进位）     对于奇异局势（0，n，n）也一样，结果也是0。     任何奇异局势（a，b，c）都有a（+）b（+）c =0。如果我们面对的是一个非奇异局势（a，b，c），要如何变为奇异局势呢？假设 a < b< c,我们只要将 c 变为 a（+）b,即可,因为有如下的运算结果: a（+）b（+）(a（+）b)=(a（+）a)（+）(b（+）b)=0（+）0=0。要将c 变为a（+）b，只要从 c中减去 c-（a（+）b）即可。     例1。（14，21，39），14（+）21=27，39-27=12，所以从39中拿走12个物体即可达到奇异局势（14，21，27）。     例2。（55，81，121），55（+）81=102，121-102=19，所以从121中拿走19个物品就形成了奇异局势（55，81，102）。     例3。（29，45，58），29（+）45=48，58-48=10，从58中拿走10个，变为（29，45，48）。     例4。我们来实际进行一盘比赛看看：         甲:(7,8,9)->(1,8,9)奇异局势         乙:(1,8,9)->(1,8,4)         甲:(1,8,4)->(1,5,4)奇异局势         乙:(1,5,4)->(1,4,4)         甲:(1,4,4)->(0,4,4)奇异局势         乙:(0,4,4)->(0,4,2)         甲:(0.4,2)->(0,2,2)奇异局势         乙:(0,2,2)->(0,2,1)         甲:(0,2,1)->(0,1,1)奇异局势         乙:(0,1,1)->(0,1,0)         甲:(0,1,0)->(0,0,0)奇异局势         甲胜。叉乘法求任意多边形面积语法：result=polygonarea(Point *polygon,int N);参数：*polygon：多变形顶点数组N：多边形顶点数目返回值：多边形面积注意：  支持任意多边形，凹、凸皆可 多边形顶点输入时按顺时针顺序排列源程序： typedef struct     double x,y; Point; double polygonarea(Point *polygon,int N)    int i,j;    double area = 0;    for (i=0;i<N;i+)         j = (i + 1) % N;        area += polygoni.x * polygonj.y;        area -= polygoni.y * polygonj.x;            area /= 2;    return(area < 0 ? -area : area);求三角形面积语法：result=area3(float x1,float y1,float x2,float y2,float x3,float y3);参数：x13：三角形3个顶点x坐标y13：三角形3个顶点y坐标返回值：三角形面积注意：  需要 math.h源程序：  float area3(float x1,float y1,float x2,float y2,float x3,float y3)    float a,b,c,p,s;    a=sqrt(x1-x2)*(x1-x2)+(y1-y2)*(y1-y2);    b=sqrt(x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);    c=sqrt(x3-x2)*(x3-x2)+(y3-y2)*(y3-y2);    p=(a+b+c)/2;    s=sqrt(p