周教授大數模板【完整版】

/*這是周尚超教授多年前寫的，希望看過此文的能記住有這麼一位元學家一直在默默地奉獻！！簡單的介紹一下我們教授吧，文革期間自學外語，精通英語，會 日，俄，德，法五種外語，早期的優秀知識分子，一個沒有讀過大學的青年於1983年自學考上了武漢大學的研究生，畢業後從事數學事業，在國外多家學術報發表過各種學術文章，在學術上頗有建樹，有著國際數學家的榮譽，帶出了一大批的優秀人才，不扯淡了....上次有幸到教授家吃飯，真的很感謝...今天聖誕了，希望教授的身體能夠好起來...這輩子能成為教授的學生，也是一種幸福...*//*高精度計算模板函數 2011-9-5  本檔案名稱 bigint.h整數 大於 2^64-1, 稱為大整數，要用高精度計算例如 a=123444455556666777788889999 是27位元用數組表示它 ，最低4位用 a[1]=9999表示，稱a[1]是10000進位數的個位元。a[2]=8888;a[3]=7777，a[4]=6666,a[5]=5555,a[6]=4444,a[7]=123;a[0]=7;表示它是萬進位 數的7位元。Base=10000進位 。大整數用字串讀入， 再用函數 tran 轉為數組本檔案有下列函數trans,comp,PN,copy,add,sub,mult,mult,Div,Div,gcd程式中用到高精度計算時，   (1)將本檔案複製到 int main()的前面, 然後將不需要的函數刪除，也可以不刪除.   (2) 也可以加上 #include<bigint.h>例如 hdu1002 A + B Problem II，只用到加法,可將sub,mult,mult,Div刪除 要用高精度計算的題hdu1002 A + B Problem IIhdu1502 Regular Wordshdu1042 n!    n<=10000hdu1133  Buy the Ticket   http://acm.hdu.edu.cn/showproblem.php?pid=1133zju2061  http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=2061hdu1134 Game of Connections           catalan數 h[1]=1,h[n]=h[n-1]*(4n-2) /(n+1) ,n<=100 hdu1250 Hat's Fibonacci http://acm.hdu.edu.cn/showproblem.php?pid=1250pku2084 Game of Connections = hdu1134pku2515 Birthday Cake           計算 1+2^M +...+n^M,     n<=10^41, M<=100pku1220 NUMBER BASE CONVERSIONpku2680 Computer Transformationpku3181 Dollar Dayzpku3199 Uncle Jackpku3324 Lucas-Lehmer Testpku3331 The Idiot of the Year Contest*/#include <vector>#include <list>#include <map>#include <set>#include <deque>#include <stack>#include <bitset>#include <algorithm>#include <functional>#include <numeric>#include <utility>#include <sstream>#include <iostream>#include <iomanip>#include <cstdio>#include <cmath>#include <cstdlib>#include <ctime>#include<stdio.h>#include<string.h>using namespace std;const int Base=10000;//萬進位 int 最大,如果只用 +,-,可加大Base=1000000000int ONE[]={1,1}; //大整數 1int ZERO[]={1,0}; //大整數 0int comp(int *a, int *b);//大整數比較大小int comp(int *a, int *b,int t);//大整數比較大小void trans(char *ch,int *A);//串->數 base=10000,A[1]是最低位,A[0]是位元(段數) void PN(int *a);//printf 大整數avoid copy(int *a, int *b);// a=bvoid add(int *a,int *b,int *c);//大整數+大整數    c=a+bvoid sub(int *a, int *b, int *c);//大整數a-大整數 b,a>bvoid mult(int *a, int b, int *d);//大整數*整數   c=a*bvoid mult(int *a,int *b,int *c);//大整數*大整數  c=a*bint Div(int *a,int b,int *c);// c=a/b   k=a%bvoid Div(int *a,int *b,int *c,int *r);// c=a/b   r=a%bvoid gcd(int *a,int *b,int *d);// d=gcd(a,b)/*串-->數 base=10000,A[1]是最低位,A[0]是位元(段數)  串12355558888-->A[1]=8888,A[2]=5555,A[3]=123  A[0]=3,每段 seg=4個數字      */void trans(char *ch,int *A){int L,i,j,k,s,i9,seg=4;  //seg=6;//Base==1000000  L=strlen(ch);  s=L-seg;k=1;  for(i=s;i>=0;i-=seg,k++)  {i9=i+seg;A[k]=ch[i]-'0';   for(j=i+1;j<i9;j++)A[k]=A[k]*10+ch[j]-48;  }  i+=seg;  A[k]=0;  for(j=0;j<i;j++)A[k]=A[k]*10+ch[j]-48;  if(A[k])A[0]=k;else A[0]=k-1;}//大整數比較大小int comp(int *a, int *b) {    int i;    if (a[0] > b[0])return 1;    if (a[0] < b[0])return -1;    for (i = a[0]; i >= 1; i--) {        if (a[i] > b[i])return 1;        if (a[i] < b[i])return -1;    }    return 0;}int comp(int *a, int *b,int t) {    int i,j=1;    for (i = b[0]; i >= 1; i--,j++) {        if (a[t] > b[i])return j;        if (a[t] < b[i])return -1;        t--;    }    return 0;}void PN(int *a){ int i;//輸出 "%04d" 10000進位 printf("%d", a[a[0]]); for(i=a[0]-1;i>=1;i--)printf("%04d", a[i]);printf("\n");}void copy(int *a, int *b){ int i;  //a=b  for (i = 0;i<= b[0]; ++i)a[i] = b[i];}//大整數+大整數   c=a+b,Base進位 ,Base=10000// c[0] 是段數,c[1]是最低位。// 例9912348888, c[0]=3,c[1]=8888,c[2]=1234,c[3]=99 void add(int *a, int *b, int *c) {  int s, i,t,p,d[1000];  if((b[0]==1)&&(b[1]==0)){copy(c,a);return;}  if((a[0]==1)&&(a[1]==0)){copy(c,b);return;}  if (a[0] >= b[0]) { copy(c,a);copy(d,b);}  else { copy(c,b);copy(d,a);}  s=c[0];t=d[0];//c[0]>=d[0]  c[s + 1] = 0;  for (i = 1; i <= t; i++)   {c[i]+=d[i];    if (c[i]>=Base){c[i]-=Base;c[i+1]++;}   }  for (; i <= s; i++){if(c[i]>=Base){c[i] -= Base;c[i+1]++;}  else break;}  if (c[s+1]>0)c[0] = s + 1;//處理最後一位}//大整數a-大整數 b, a>b，c=a-bvoid sub(int *a, int *b, int *c){ int i, j,p=0;  if(b[0]==1&&(b[1]==0)){copy(c,a);return;}  if(comp(a,b)<=0){copy(c,ZERO);return;}  if(comp(a,ZERO)==0){copy(c,ZERO);return;}  for (i = 1; i <= b[0]; i++)  { c[i] = a[i] - b[i] - p;    if (c[i] < 0)    { c[i] += Base;  p = 1;    } else p = 0;  }  for (; i <= a[0]; i++)  { c[i] = a[i] - p;    if (c[i] < 0)    { c[i] += Base; p = 1;    } else p = 0;  }  for(i=a[0];i>=1;i--){    if (c[i]) {c[0] = i;break;}  }}//大整數a乘常數 b b<Base, d=a*bvoid mult(int *a,int b,int *d){ int w, i,p;  if(b==1){copy(d,a);return;}  if((b==0)||(a[0]==1&&a[1]==0)){d[0]=1;d[1]=0;return;}  w = a[0];  p = 0; //  for (i = 1; i <= w; i++)  { d[i] = a[i] * b + p;    if (d[i] >= Base)    { p = d[i] / Base; d[i] %= Base;}else p=0;  }  if (p) {w++;d[w] = p;}  d[0] = w;}//大整數a*大整數 b, c=a*bvoid mult(int *a,int *b,int *c){int i,j,s,m,n,k,p;//Base<=10000 m=a[0];n=b[0];k=m+n-1;/*m位元 *n位元 = m+n-1位 或 m+n 位元  執行個體  12345678*123456789012 =1524157764056090136  a[0]=2  a[1]=5678;a[2]=1234;  b[0]=3  b[1]=9012;  b[2]=5678;b[3]=1234;    */ p=0; for(i=0;i<=k;i++)c[i]=0;c[i]=0;//  a[1],a[2]....,a[m]//  b[1]//  c[1],c[2],...,c[m] for(j=1;j<=n;j++) { s=j;   if(b[j])for(i=1;i<=m;i++,s++)   {c[s]+=a[i]*b[j];    if(c[s]>=Base){c[s+1]+=c[s]/Base;c[s]%=Base;}   } } for(i=1;i<=k;i++)if(c[i]>=Base) {  c[i+1]+=c[i]/Base;c[i]%=Base; } if(c[k+1])k++; c[0]=k;}// c=a/b   k=a%bint Div(int *a,int b,int *c){int i,j,s,t; int k=0,q;// long long k=0,q;  Base>=100000 if(comp(a,ZERO)==0){copy(c,ZERO);return 0;} for(i=a[0];i>=1;i--)   {q=k*Base+a[i];  c[i]=q/b;k=q%b;   } if(c[a[0]]==0) { c[0]=a[0]-1; }else c[0]=a[0]; return k;}//  大整數除大整數  c=a/b, r=a%bvoid Div(int *a,int *b,int *c,int *r){int m=a[0],n=b[0],i,j,k,s,t,p,q,x,y; int Base1000=Base*1000,Base100=Base*100,Base10=Base*10; //m位元除 n位元= k=m-n+1位元 if(comp(a,b)<0){copy(r,a);copy(c,ZERO);return;}    k=m-n+1; for(s=k;s>=0;s--)c[s]=0; for(s=k,t=m;s>=1,t>=n;) {  p=comp(a,b,t); // a[t]-b[n] /////////////////  if(p<0 && t==n)break;  if(p==0){c[s]++;   for(j=t;j>=t-n+1;j--)a[j]=0;     t-=n;if(t<n)break;s-=n;  }                if(p>1)  {   c[s]++;for(j=t,q=n;j>=t-n+1;j--,q--)a[j]-=b[q];   for(j=t-n+1;j<=t;j++){if(a[j]<0){a[j]+=Base;a[j+1]--;}   //////////////////   }   for(j=t;j>=t-n+1;j--)if(a[j])break;   if(j<n)break;   s-=t-j;   t=j;  }  if(p<0 && t>n){a[t-1]+=a[t]*Base;a[t]=0;t--;s--;}  if(p==1)  //   9999 1234 /  1999 9999  x=4  {    x=a[t]/(b[n]+1);    c[s]+=x;//x< 4999    for(j=t,q=n;q>=1;q--,j--)a[j]-=x*b[q];                                   for(j=t-n+1;j<=t;j++)    {     while(a[j]<=-Base1000){a[j]+=Base1000;a[j+1]-=1000;}     while(a[j]<=-Base100){a[j]+=Base100;a[j+1]-=100;}     while(a[j]<=-Base10){a[j]+=Base10;a[j+1]-=10;}     while(a[j]<0){a[j]+=Base;a[j+1]--;}    }  }              } if(c[k]==0)c[0]=k-1; else c[0]=k; for(i=m;i>=1;i--)if(a[i])break; if(i==0)copy(r,ZERO); else { r[0]=i; while(i>=1){r[i]=a[i];i--;} }}// d=gcd(a,b)void gcd(int *a,int *b,int *d){int p,x[10000],y[10000],c[10000]; p=comp(a,b); if(p<0){copy(x,b);copy(y,a);}    else {copy(x,a);copy(y,b);}    // x>=b if(comp(y,ZERO)==0) {copy(d,x);return;} copy(d,y); while(1) {  Div(x,d,c,y);// c=x/d; y=x%d  if(comp(y,ZERO)==0) {return;}  if(comp(y,ONE)==0) {copy(d,ONE);return;}   copy(x,d);  copy(d,y); }}