《程序员密码学》之大数算数&Eratosthenes素数筛选

来源:互联网 发布:测试词汇量网站 知乎 编辑:程序博客网 时间:2024/06/05 20:53

10进制大数算法,支持+、-、*、/、%以及一般的比较运算符,支持字符串方式读入以及输出大数

#pragma once#include <iostream>#include <string>using namespace std;#define SWAP(x,y,t) ((t)=(x),(x)=(y),(y)=(x))#define MIN(x,y) (x)>(y)?(y):(x)#define MAX(x,y) (x)>(y)?(x):(y)#define BASE 10000       //大数以每位十进制数的方式存储//乘法实现,计算A*B,结果返回到 C,函数结果为乘法结果的长度#define MAX_SIZE 128//一般将算术和关系操作符定义非成员函数,而将赋值操作符定义为成员://大数以小端格式存储struct BigNum{int len;    int v[MAX_SIZE];BigNum(){    len=0;memset(v,0, sizeof(v));}BigNum(const char* num){*this=num;}BigNum(int num){len=0;for(; num>0; num/=BASE) v[len++] = num%BASE;}BigNum&operator=(const char* num){intl=strlen(num);intn=(l%4)?(l%4):4;len=0;charx[5]={0};memcpy(x,num, n );v[(l+4-1)/4-1]= atoi(x);len++;for(int i=(l+4-1)/4-2,j=n; i>=0; i--,j+=4){memcpy(x, num+j, 4);v[i]= atoi(x);len++;}return *this;}};booloperator<(const BigNum& left,const BigNum& right){int i;if(left.len!=right.len)return left.len<right.len;for(i=left.len-1; i>=0 && left.v[i]==right.v[i]; i--);return i>0? left.v[i]<right.v[i]: 0;}booloperator>(const BigNum& left,const BigNum& right){return right<left;}booloperator==(const BigNum& left,const BigNum& right){return (!(left<right)) && (!(right<left));}booloperator<=(const BigNum& left,const BigNum& right){return (left<right) || (left==right);}booloperator>=(const BigNum& left,const BigNum& right){return right<=left;}BigNumoperator+(const BigNum& left,const BigNum& right){int n = MAX(left.len,right.len);int x=0,i;BigNumresult;for(i=0; i<n; i++){if(i<left.len)  x+=left.v[i];if(i<right.len) x+=right.v[i];result.v[i] = x%BASE;x/=BASE;}if(x!=0)result.v[i++] = x;result.len= i;return result;}BigNumoperator-(const BigNum& left,const BigNum& right){BigNumresult;int i,x=0;for(i=0; i<left.len; i++){result.v[i] = left.v[i] - x;if(i<right.len) result.v[i] -= right.v[i];if(result.v[i]<0)result.v[i] += BASE, x = 1;elsex=0;}result.len= i;while(result.len>0&&result.v[result.len-1]==0)result.len--;return result;}BigNum  operator*(const BigNum& left,const BigNum& right){        BigNum  result;         int             c0,c1,c2;//存数临时结果         int             tx,ty;         int             len=left.len+right.len;         int             *templ,*tempr;        c0 = c1 = c2 = 0;         for(int i=0; i<len; i++)//i为本次计算在结果中所处的级别,这里一个整形为一个级别        {                ty              = MIN(i, right.len-1);                tx              = i-ty;                templ   = (int*)(left.v + tx);                tempr   = (int*)(right.v+ ty);                 int n = MIN( left.len-tx, ty+1);//本级别可计算的次数                c0 = c1; c1 = c2; c2 = 0; //上次计算结果向下移一个级别,进入本次计算                 for(int j=0; j<n; j++)                {                         int x=*templ++, y=*tempr--; c0 += x*y;            /*_asm   //mul operation            {                                   pushad                mov eax,x                mov ebx,y                mul ebx          //无符号乘法指令                add c0,eax                //adc c1,edx                //adc c2,0                popad                        }*/                }                result.v[i]=c0%BASE;                c1+=c0/BASE;                c2+=c1/BASE;        }         //计算长度         for(int i=len-1; i >= 0; i--)                 if (result.v[i]==0)                        len--;                 else                         break;        result.len      = len;         return result;}BigNumoperator/(const BigNum& a,const BigNum& b){BigNum tmp, mod, res;int i, lf, rg, mid,base=BASE;mod.v[0] = mod.len = 0;for (i = a.len - 1; i >= 0; i--) {//从最高4位开始,每次增加一个4位试商mod = mod * base + a.v[i];//二分法寻找mod/b的商for (lf = 0, rg = base -1; lf < rg; ) {mid = (lf + rg + 1) / 2;if (b * mid <= mod) lf = mid;else rg = mid - 1;}//存储结果并将剩余的数继续试商res.v[i] = lf;mod = mod - b * lf;}res.len = a.len;while (res.len > 0 && res.v[res.len - 1] == 0) res.len--;return res; // return mod 就是%运算}BigNumoperator%(const BigNum& a,const BigNum& b){BigNum tmp, mod, res;int i, lf, rg, mid,base=BASE;mod.v[0] = mod.len = 0;for (i = a.len - 1; i >= 0; i--) {//从最高4位开始,每次增加一个4位试商mod = mod * base + a.v[i];//二分法寻找mod/b的商for (lf = 0, rg = base -1; lf < rg; ) {mid = (lf + rg + 1) / 2;if (b * mid <= mod) lf = mid;else rg = mid - 1;}//存储结果并将剩余的数继续试商res.v[i] = lf;mod = mod - b * lf;}res.len = a.len;while (res.len > 0 && res.v[res.len - 1] == 0) res.len--;return mod;// %运算}ostream& operator<<(ostream& os,const BigNum& num){        printf( "%d",num.v[num.len-1]);         for(int i=num.len-2; i>=0; i--)                printf( "%04d",num.v[i]);         return os;}istream& operator>>(istream& in,BigNum& num){string str;in>>str;num=str.c_str();return in;}



X86版本的大数乘法简陋实现,没有书上的可移植性,结果以16进制存储

只是为了试试效果~

#include <stdio.h>#define MIN(x,y) (x)>(y)?(y):(x)/*#define MUL(x,y) _asm {\                          "pushad\r\n"\                          "mov eax,x\r\n" \                          "mov ebx,y\r\n" \                          "mul ebx\r\n"\                          "add c0,eax\r\n" \                          "adc c1,edx\r\n" \                          "adc c2,0\r\n"\                          "popad\r\n"\                          }*///乘法实现,计算A*B,结果返回到C,函数结果为乘法结果的长度,大数整体采用小端模式存储int mul(int *A,int lena, int *B,int lenb, int *C){   int    ix, iy, iz, tx, ty, pa;   int    c0, c1, c2, *tmpx, *tmpy;   c0=c1=c2=0;   /* get size of output and trim */   pa = lena+lenb;     for (ix = 0; ix < pa; ix++) {      /* get offsets into the two bignums */      ty = MIN(ix, lenb-1);      tx = ix - ty;      /* setup temp aliases */      tmpx = A + tx;      tmpy = B + ty;      /* this is the number of times the loop will iterrate, essentially its         while (tx++ < a->used && ty-- >= 0) { ... }       */      iy = MIN(lena-tx, ty+1);      /* execute loop */      c0=c1;c1=c2;c2=0;      for (iz = 0; iz < iy; ++iz) {          int i=*tmpx++,j= *tmpy--;                  _asm   //mul operation                  {                          pushad                          mov eax,i                          mov ebx,j                          mul ebx                          add c0,eax                          adc c1,edx                          adc c2,0                          popad                  }      }      /* store term */     C[ix]=c0;  }         return pa;}//平方计算,函数结果为乘法结果的长度int square(int *A,int lena, int *C){   int    ix, iy, iz, tx, ty, pa;   int    c0, c1, c2, *tmpx, *tmpy;   c0=c1=c2=0;   /* get size of output and trim */   pa = lena*2;    for (ix = 0; ix < pa; ix++) {      /* get offsets into the two bignums */      ty = MIN(ix, lena-1);      tx = ix - ty;      /* setup temp aliases */      tmpx = A + tx;      tmpy = A + ty;      /* this is the number of times the loop will iterrate, essentially its         while (tx++ < a->used && ty-- >= 0) { ... }       */      iy = MIN(lena-tx, ty+1);          iy = MIN(iy, (ty-tx+1)>>1);      /* execute loop */      c0=c1;c1=c2;c2=0;      for (iz = 0; iz < iy; ++iz) {          int i=*tmpx++,j= *tmpy--;                  _asm   // 由于是平方,加两次,省去一半的乘法指令                  {                          pushad                          mov eax,i                          mov ebx,j                          mul ebx                          add c0,eax                          adc c1,edx                          adc c2,0                                                  add c0,eax                                                  adc c1,edx                                                  adc c2,0                          popad                  }      }          if ((ix&1)==0)  // 偶数位的话需要加上中间未重复的结果          {          int i=A[ix>>1],j= A[ix>>1];          _asm   //mul operation          {                                     pushad                      mov eax,i                               mov ebx,j                         mul ebx              add c0,eax              adc c1,edx              adc c2,0              popad          }          }      /* store term */     C[ix]=c0;  }         return pa;}//求幂,lena:A的长度(4字节为单位),exponent:指数,C:存储结果//这个版本的求幂只能求2^n形式的幂int power(int *A,int lena, int exponent, int *C){         int len=lena,*temp,*dst=C;         int squ=0;// 可直接平方的次数         for(int t=exponent; t!=1; t>>=1)                squ++;         for(int i=0; i < squ; i++)        {                memset(C, 0 , 1024);                len = square(A,len, C);                SWAP(A,C,temp);        }         if (dst!=A) memcpy(dst,A,len);         return len;}void print_bignum(int *bignum,int u){         int first=1;        printf( "0x");         for(int i=u-1; i >= 0; i--)                 if ((bignum[i]==0)&&first)                         continue;                 else                {                        first = 0;                        printf( "%08X",bignum[i]);                }}int main(){         int a[] = {0x87654321,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,\                0x87654321,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,\                0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678};         int b[] = {0x87654321,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,\                0x87654321,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,\                0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678,0x12345678};         int c[256]={0};         int len=mul(a,sizeof (a)/sizeof( int),b,sizeof (b)/sizeof( int), c);        print_bignum(c,len);        getchar();         return 0;}

 


 
 
素数筛选
#include <iostream>#include <cmath>#include <vector>using namespace std;//查找指定范围内的所有素数//100以内的所有素数unsigned int prime_tab[] = {         2,  3,  5,  7,  11,  13,  17,        19, 23, 29, 31,  37,  41,  43,        47, 53, 59, 61,  67,  71,  73,        79, 83, 89, 97};//Eratosthenes 素数筛选,取得0-num内所有的素数,c++容器版int Eratosthenes(vector<int >& tab, unsigned int num){#define TEST(p,x) (((p)%(x))!=0)       //是否满足素数条件                 if(num <= 100)        {                 for(int i=0; i < sizeof(prime_tab)/ sizeof(int ); i++)                {                         if(prime_tab[i]<num)                                tab.push_back(prime_tab[i]);                         else                                 break;                }                 return tab.size();        }         unsigned int i=(int)sqrt(( double)num);        Eratosthenes(tab, i); /*递归初始化i以内的素数表 */         for(; i < num; i++)// 筛选        {                 bool flag=false ;                 int  n=tab.size();                 for(vector<int >::iterator iter=tab.begin(),end=iter+n; iter != end ; iter++)                {                         if( !TEST(i,*iter) )                                 break;                         if( (iter+1)==end )                                flag= true;                }                 if(flag) tab.push_back(i);        }#undef TEST(p,x)         return tab.size();}//Eratosthenes 素数筛选,取得 num内所有的素数int Eratosthenes2(unsigned int num){         unsigned int count=0;                 char *tab = new char[num+1];        memset(tab+2,1, num-2);         for(unsigned int i=2; i <= ( unsigned int )sqrt((double)num); i++)                 if(tab[i]==1)                         for(unsigned int j=i; j*i <= num; j++)                                tab[i*j] = 0;         for(unsigned int i=2; i <= num ; i++)                 if(tab[i]==1)                {                        printf( "%d\t",i);                        count++;                }        printf( "\n%d  ",count);         delete [] tab;         return count;}int main(){        vector< int> tab;        //素数表容器        tab.reserve(2000);       //预留空间        cout<<Eratosthenes(tab,150000); //函数返回统计个数         for(int i=0; i < tab.size(); i++)  //输出素数                printf( "%d\t", tab[i]);        getchar();         return 0;}



 

原创粉丝点击