大数运算（能开方、求幂、排列组合、取对数、求三角值）

趁着寒假想用C写一个模拟星体运动的程序，发现星体的数据过于庞大，于是便开始写大数运算的程序。

目前大一，难免有些错误，还请各位猿友能指出。

#include<math.h>
#include<string>
#include<iostream>


#define pi asin(1)*2
#define lon 21//数据精度
#define d 1000//数据长度
/*大数表示演示：
* 12 345 678.123 456 7->long a[lon]={2,12,345,678,123,456,700}
* 负数的转换不同，如-123 456->{1,-124,544}即-124000+544=123456
* 我们把数组a里的一个数字成为一个节，如：第0节->2，第1节->12
* 数组a表示的数A可以这样来表示：
*{
* int A=0,i=0;
* for(i=0;i<lon;i++)A+=a[i+1]*pow(d,*a-i);
* }
* 定义函数时用b和n来表示大数（big num）和普通的数（num）
* 定义函数参数时用n和num表示大数和普通数
*/
long bn0[lon] = { 0 };//0大数
long bn1[lon] = { 1,11,110 };
long bn2[lon] = { 2 ,-1,-2,-3,-4 };
long bpi[lon] = { 0,3,141,592,653,589,793,115,997,963,468,544,185,161,590,576,171,875 };
long be[lon] = { 0,02,718,281,828,459,45,073,32,29,904,425,144,195,556,640,625 };
// 4,5,965000
//一段表示3位数

//复制大数，X用于储存返回值
void copyb(long* n, long X[]) {
    int i ;
    for (i=0;i<lon;X[i]=n[i++]);
}

//复制大数的某一段，st->start,fi->finish
void copybb(long* n, long X[],int st=1,int fi=lon) {
    st = st <= 0 ? 1 : st;
    st = st > lon ? lon - 1 : st;
    long n2[lon] = { 0 };int i, j;
    for (i = st;i < fi;i++)n2[i - st + 1] = n[i];
    *n2 = *n - st + 1;
    copyb(n2, X);
}
//大数归零
void zerob(long* n) {
    copyb(bn0, n);
}

//平的读音，整理大数使大数符合格式，还有复制大数的功能
void ping(long* n,long X[]=0) {
    int k, i, f, num = n[1];
    if (X == 0)X = n;
    for (k = 1;k < lon && !n[k];k++);k--;//查询头0在何处
    if (k == lon - 1) { copyb(n,X);return; };//如果头0在lon-1则是〇大数
    for (i = 1;i < lon - k;i++)n[i] = n[i + k];//除去头0
    *n = *n - k;
    num = n[1];
    for (f = 0;num /= d;f++);f;num = n[1];//检查头是否超出d
    for (i = lon - 1;i > f;i--)n[i] = n[i - f];*n = *n + f;
    for (i = f + 1;i >= 1;i--) { n[i] = num % d; num /= d; }

    for (i = lon - 1;i >= 2;i--) {
        if (n[i] >= d)
        {
            n[i - 1]++; n[i] -= d;
        }
        if (n[i] < 0) { n[i - 1]--; n[i] += d; }
    }
    if (!n[1] || abs(n[1]) >= d) ping(n,X);
    for (i = 2;i < lon;i++)
        if (abs(n[i]) >= d || n[i] < 0)ping(n, n); 
    {copyb(n, X);return;}
}

//与下面的prnum(print num)配套，补〇
void seentzer(long n) {
    int i = 2;
    (n /= 10) ? i-- : 0;(n /= 10) ? i-- : 0;
    (n /= 10) ? i-- : 0;
    for (;i>0;i--)printf("0");
}
//打印大数
void prnum(long* n,std::string str1="\n", std::string str2 = "\n") {
    std::cout << str1;
    int i = n[0];
    printf("[%d]", n[0]);
    if (-1 == *n)printf("0.");
    for (i = 0;i < lon-1;i++) {
        seentzer(n[i+1]);
        printf("%d,", n[i+1]);
        if (i == *n)printf(".");
    }
    std::cout << str2;
}
//大数相加
   void addbb(long* n1, long* n2,long X[]=0) {//1
    ping(n1,n1);ping(n2,n2);
    if (X == 0)X = n2;
    long in[lon * 2] = { 0 }, * a;
    int i, j;
    if (*n1 < *n2) { a = n1;n1 = n2;n2 = a; }
    *in = *n1;
    for (i = 1;i < lon ;i++) {//2
        if (i >*n1 - *n2 && i < lon) { in[i] = n1[i] + n2[i - *n1 + *n2]; }
        if (i <= *n1 - *n2)in[i] = n1[i];
        if (i >= lon)
            in[i] = n2[i + *n2 - *n1];
    }//2
    ping(in, X);
}//1

   
   //大数加普通数
   void addnb(long* n1, long double num,long X[]=0) {
       if (X == 0)X = n1;
       long a[lon] = { 0 }, * n2;
       int in, i, f = 1;long double fl;
       in = (int)num;
       fl = num;
       a[0] = 0;
       a[1] = in;
       //0.123,456,789,1
       for(i=2;i<lon;i++){fl = (fl-a[i-1])*d;a[i] = (int)fl; }
       ping(n1,n1);ping(a,a);
       addbb(n1, a, X);
   }

   //普通数加普通数
   void addnn(long double num1, long double num2, long X[]) {
       long n[lon] = { 0 };
       addnb(bn0, num1, n);
       addnb(n, num2, X);
   }
   //大数取负(minus sign)
   void musgb(long* n,long X[]=0) {
       if (X == 0)X = n;
       int i = 1;long* n1 = n;
       X[0] = n[0];
       for (i = 1;i < lon;i++) {
           X[i] = -n[i]; 
       }
   }


   //相减（minus），n=1表示n1-n2,n=0表示n2-n1
   void mnsbb(long* n1, long* n2,long X[]=0,int n=1) {
       if (X == 0)X = n2;
       long n3[lon * 2] = { 0 },*a;
       if (n) { musgb(n2, n3), a = n1; }
       else { musgb(n1, n3);a = n2; }
       ping(a,a);ping(n3,n3);
       addbb(a, n3, n3);ping(n3, X);
   }
   //相减，n=1表示n1-num;n=0表示num-n1
   void mnsnb(long* n1, long double num, long X[]=0, int n = 1) {
       if (X == 0)X = n1;
       ping(n1, n1);
       if (n) { addnb(n1, -num, X);return; }
       else { long n3[lon * 2] = { 0 };musgb(n1, n3);
           addnb(n3, num, X);
       }
   }

   //相乘（multiply）head表示d的head次幂
   void mtpnb(long* n1,long double num,long X[]=0,int head=0) {
       if (X == 0)X = n1;
       int in=(int)num;long double fl = num - in;
       int i,j;
       long sum[lon*2] = { 0 };
       
           long n2[lon] = { 0 };
           ping(n1, n2);
           *n2 += head;
           *sum = *n2;
           for (i = 1;i < lon;i++)sum[i] =n2[i]*in;
           ping(sum, sum);
           //addbb(sum, n2, sum);
       
       for (j = 1;j<lon;j++) {
           copyb(n1, n2);n2[0] = *n2-j+head;
           fl *= d;
           for (i = 1;i < lon;i++) { n2[i] *=(int)fl; }
           
           addbb(sum, n2, sum);
           fl = (fl-(int)fl) ;
           
       }
       ping(sum, X);
   }

   void mtpbb(long* n1, long* n2, long X[]=0) {
       float cmpnb(long* n1, long double num);
       if (X == 0)X = n2;
       long n3[lon]={0},sum[lon]={0};
       long n4[lon] = { 0 };
       int i, j;ping(n1, n4);
       for (i = 1;i < lon;i++) {
           ;n4[0] =*n1+ 1-i;
           mtpnb(n4, n2[i], n3);
           addbb(n3, sum, sum);
       }
       ping(sum, X);
   }
   void mtpnn(long double num1, long double num2, long X[]) {
       long n[lon] = { 0,1 };
       mtpnb(n, num1, n);
       mtpnb(n, num2, X);
   
   }
   //相除（divide）
   void dvdbb(long* n1, long* n2, long x[]=0, int n = 1) {
       float cmpnb(long* n1, long double num);
       long X[lon] = {0};
       
       long sum[lon] = { 0 }, n3[lon] = { 0 }, n4[lon] = { 0 };
       int i,j;ping(n1, n4);ping(n2, n3);
       int f=1;
       if (cmpnb(n3, 0) < 0) { musgb(n3, n3); f *= -1; }
       if (cmpnb(n4, 0) < 0) { musgb(n4, n4); f *= -1; }
       X[0] = *n1 - *n2;
       *n4 = 0;*n3 = 0;
       for (j = 1;j < lon;j++) {
           for (i = 1;i < d;i++) {
               mnsbb(n4, n3, n4);
               ping(n4, n4);
               if (n4[1] < 0||i>=d) { i-=1;break; }}
           addbb(n4, n3, n4);
           n3[0] -= 1;
           X[j] = i; }
       if(x)mtpnb(X,f,x);
       else mtpnb(X, f, n2);
       }
   
   void dvdnn(long double num1, long double num, long X[], int n = 1) {

       long sto;int i,j,k;
       long n1[lon] = { 0 }, n2[lon] = { 0 };
       if (!n) { num1 += num;num = num1 - num;num1 -=  num; } 
       addnb(bn0, num1,n1);addnb(bn0, num, n2);
       dvdbb(n1, n2, X);

   }

   void dvdnb(long *n1,long double num,long X[]=0,int n=1) {
       if (X == 0)X = n1;

       int i;long sum[lon] = { 0 },n2[lon]={0},n3[lon]={0};
        addnb(bn0, num, n2);
       if (n) { dvdbb(n1, n2, X); }
       else { dvdbb(n2, n1, X); }
   }

  
   long double mysqr(long double a, int m) {
       long double x = 1;int i,j = 300;

       for (;j--;printf("\n迭代%d次->%.36f\n", j, x))
           for (i = j, x = 1;i--;)
           //    x = (a * (m - 1) / x + pow(x, m -1)) / m;
            x = (  a *pow(x, 1 - m) +(m - 1) * x) / m;
       return x;
   }
  

   //整数型次幂（即num为int型）平方,n等于0的那一功能还没设计好
   void powin(long* n1, int num, long X[] = 0, int n = 1) {
       if (X == 0)X = n1;long n2[lon] = { 0 }, x[lon] = { 0 };
       if (num < 0) { num *= -1;dvdnb(n1, 1, n2, 0); }
       else ping(n1, n2);

       if (n) {
           if (num == 0) { zerob(X);addnb(X, 1, X);return; }
           for (ping(n2, x);--num;)mtpbb(n2, x, x);
       }
       else {
           long n2[lon] = { 0 };ping(n1, n2);
           for (;;)mtpnb(X, num, X);
       }ping(x, X);
   }
 
   //整数型次方（即num为int型）平方
 void sqrin(long* n1, int num, long X[] = 0) {
       int i, j = 100;long x[lon] = {0,1}, re[4][lon] = {0};
       if (X == 0)X = n1;
       for (;j--;)
       {// x = (  a *pow(x, 1 - m) +(m - 1) * x) / m;
           powin(x, 1 - num, re[1]);//2
           mtpbb(re[1], n1, re[1]);//3
           mtpnb(x, num - 1, re[2]);//4
           addbb(re[1], re[2], re[3]);//5
           dvdnb(re[3], num, x);//6
       }
       ping(x, X);
   }
 //考录到幂很容易把大数弄炸导致精度不高，而且也不怎么常用，故未写小数型幂和开方

 //取余数（remainder）
       void rmdbb(long* n1, long *n2, long X[] = 0,int n=1) {
            long *n3,*n4,re[lon]={0};
            if (n) { n3 = n1;n4 = n2; }
            else { n3 = n2;n4 = n1; }
            ping(n3, n3);ping(n4, n4);
            dvdbb(n3, n4, re);
            copybb(re,re,re[0]+2);
            mtpbb(re, n4, re);
           // mnsbb(n3, re, re);
            ping(re, X); 
   }
       void rmdnb(long* n1, long double num, long X[] = 0, int n = 1) {
           long n2[lon] = { 0 };
           addnb(bn0, num, n2);
           rmdbb(n1, n2, X, n);
       }
       //比较大小
       float cmpbb(long* n1, long* n2) {
           ping(n1, n1);ping(n2, n2);
           int f, i;
           for (i = 0;i < lon;i++)  if (n1[i] != n2[i]) {
               if (n1[i] > n2[i])return i + 0.1;
               else return -i - 0.1;
           }
           return 0;
       }
       float cmpnb(long* n1, long double num) {
           long n2[lon] = { 0 };
           addnb(bn0, num, n2);
           return cmpbb(n1, n2);
       }

       //大数转普通数（big to num）
       long double bg2nu(long *n) {
           long double n1=0;
           int i;
           for (i = 1;i < lon;i++) { n1 += n[i]*pow(d,*n-i+1); }
           return n1;
       }
       //sin
       long double sinb(long* n, long X[]=0) {
           long double n1 = 0;long n2[lon] = { 0 };
           mtpnb(bpi, 2, n2);
           rmdbb(n, n2, n2);
           n1 = sin(bg2nu(n2));
           if(X)addnb(bn0, sin(bg2nu(n2)), X);
           return n1;
       }

       long double cosb(long* n, long X[] = 0) {
           long double n1 = 0;long n2[lon] = { 0 };
           mtpnb(bpi, 2, n2);
           rmdbb(n, n2, n2);
           n1 = cos(bg2nu(n2));
           if (X)addnb(bn0, cos(bg2nu(n2)), X);
           return n1;
       }

       long double tanb(long* n, long X[] = 0) {
           long double n1 = 0;long n2[lon] = { 0 };
         //  mtpnb(bpi, 2, n2);
           rmdbb(n, bpi, n2);
           n1 = tan(bg2nu(n2));
           if (X)addnb(bn0, tan(bg2nu(n2)), X);
           return n1;
       }

       long double atanb(long* n, long X[] = 0) {
           long double num;
           if (cmpnb(n, pow(10, 18)) > 0) { if (X)ping(bpi, X);return bg2nu(bpi) / 2; }
           num = atan(bg2nu(n));
           if (X) addnb(bn0, num, X);
           return num;
       }

      
       //不完全排列A（n1,n2）n1>n2
       void Abb(long* n1, long* n2, long X[], int n = 1) {
           long* n3, * n4, re[lon] = { 0 ,1 };
           long n5[lon] = { 0 };
           int f = 1;
           if (n) { n3 = n1;n4 = n2; }
           else { n3 = n2;n4 = n1; }
           if (cmpbb(n3, n4) < 0) {copybb(bn0, X);return; }
           if (cmpnb(n3, 0) < 0) { musgb(n3, n3); f *= -1; }
           if (cmpnb(n4, 0) < 0) { musgb(n4, n4); f *= -1; }
           mnsbb(n3, n4, n5);addnb(n5, 1, n5);
           prnum(n5, "n5->");
           ping(n3, n3);ping(n4, n4);
           for (;cmpbb(n5,n3)<=0;addnb(n5, 1, n5)) {
               prnum(re, "re->");
               mtpbb(re, n5, re);
           }
           mtpnb(re,f, X);
       }

       //全排列（full permutation）
       void fpmtb(long*n1,long X[]) {
           Abb(n1, n1, X);
       }
       void Cbb(long* n1, long* n2, long X[], int n = 1) {
           long* n3, * n4, re[3][lon] = { 0 };
           long n5[lon] = { 0 };
           int f = 1;
           if (n) { n3 = n1;n4 = n2; }
           else { n3 = n2;n4 = n1; }
           Abb(n1, n2, re[0]);
           fpmtb(n2, re[1]);
           dvdbb(re[0], re[1], X);
       }

       void lnb(long* n1, long X[]) {
           long n3[lon]={0}, re[lon] = {1 ,1};
           long double num3,num4;
           ping(n1, n3);
           mtpnn(3*log(10), *n3, re);
           *n3 = 0 ;
           addnb(re, log(bg2nu(n3)),X);
       }

       //以10为底的对数
       long double lg(long double num) {
           return(log(num) / log(10));
       }
       void lgb(long* n1, long X[]) {
           long n3[lon] = { 0 }, re[lon] = { 1 ,1 };
           long double num3, num4;
           ping(n1, n3);
           mtpnn(3 * lg(10), *n3, re);
           *n3 = 0;
           addnb(re, lg(bg2nu(n3)), X);
       }

       long double lognn(long double num1, long double num2,  int n = 1,long X[]=0) {
           long double num3;
           num3=n? log(num1) / log(num2):log(num2) / log(num1);
           addnn(0, num3, X);
           return num3;
       }
       void lognb(long* n1, long double num, long X[], int n = 1) {
           long n2[lon] = { 0 };
           lnb(n1, n2);
           if(n)dvdnb(n2, log(num), X);
           else dvdnb(n2, log(num), X,0);
       }
       void logbb(long* n1, long* n2, long X[], int n = 1) {
           long n4[lon] = { 0 },n3[lon] = { 0 };
           lnb(n1, n3);lnb(n2, n4);
           
           if (n)dvdbb(n3, n4, X);
           else dvdbb(n4, n3, X);
       }