📄 bigintcpp.cpp
字号:
{
unsigned __int64 m;
short iwhich;
for(int i=0;i<=iIdx;i++)
iVal[i]=0;
m=0x1000000000000000;
if(which>0)
{
iwhich=(which-1)<<2;
iVal[iIdx]=m>>iwhich;
}
else
{
iIdx++;
iVal[iIdx]=1;
}
}
bool BIisZero(const BigInt & x)
{
//if the value is equal to zero return TURE,else FALSE
if(x.iIdx!=0)
return false;
if(x.iVal[0]!=0)
return false;
return true;
}
BigInt & BigInt::operator=(const BigInt & p1)
{
//operator=
iIdx=p1.iIdx;
for(int i=0;i<TOTALBOXS;i++)
iVal[i]=p1.iVal[i];
return *this;
}
//p444
BigInt & BigInt::operator=(const unsigned __int64 & p1)
{
iIdx=0;
iVal[0]=p1;
return *this;
}
inline bool operator>(const BigInt & p1,const BigInt & p2)
{
//operator>
short stop;
stop=p1.iIdx;
if(stop>p2.iIdx)
return true;
else if(stop<p2.iIdx)
return false;
else
{
while(stop>=0)
{
if(p1.iVal[stop]>p2.iVal[stop])
return true;
else if(p1.iVal[stop]<p2.iVal[stop])
return false;
stop--;
}
return false;
}
}
inline bool operator==(const BigInt & p1,const BigInt & p2)
{
//operator==
if(p1.iIdx!=p2.iIdx)
return false;
for(int i=0;i<=p1.iIdx;i++)
{
if(p1.iVal[i]!=p2.iVal[i])
return false;
}
return true;
}
//p445
//p446
BigInt operator-(const BigInt & p1,const BigInt & p2)
{
//operator-
int i;
BigInt result;
unsigned __int64 carry;
if(BIisZero(p2))
return p1;
else{
carry=0;
for(i=0;i<=p1.iIdx;i++)
{
if(p1.iVal[i]>=p2.iVal[i])
{
result.iVal[i]=p1.iVal[i]-p2.iVal[i];
if((result.iVal[i]==0x0)&&(carry==0x1))
result.iVal[i]=FULL;
else
{
result.iVal[i]-=carry;
carry=0x0;
}
}
else
{
result.iVal[i]=FULL-p2.iVal[i]+p1.iVal[i]+RMASK[1]-carry;
carry=RMASK[1];
}
}
result.iIdx=p1.iIdx;
while((result.iIdx>0)&&(result.iVal[result.iIdx]==0))
result.iIdx--;
}
return result;
}
//p447
//p449
//两大整数相除（参考BIdivision）
BigInt operator /(const BigInt & p1,const BigInt & p2)
{
//operator /
//result=a/b....d
BigInt d;
return BIdivision(p1,p2,d);
}
//返回余数的两大整数相除
BigInt Remainder(const BigInt & p1,const BigInt & p2,BigInt & p3)
{
//operator /
//result=a/b....d
// BigInt d;
return BIdivision(p1,p2,p3);
}
BigInt operator %(const BigInt & p1,const BigInt & p2)
{
//operator%
//C=p1/p2+result
BigInt a,mb;
int test;
a=p1;
if(BIcompare(a,p2)>0)
{
//Get a number(mb) that is largest multiple of (p2) less then(a).
mb=p2; //make mb close to a quickly
MakeThemClose(a,mb);//?????????????????????
while(BIcompare(a,mb)>0)
mb<<=1;
while(BIcompare(a,p2)>0)
{
while(BIcompare(a,mb)<0)
mb>>=1;
a=a-mb;
}
}
if (BIcompare(a,p2)==0)
a.setZero();
return a;
}
//p450
void BIswap(BigInt & p1,BigInt & p2)
{
BigInt tmp;
tmp=p1;
p1=p2;
p2=tmp;
}
BigInt BIdivision(const BigInt & p1,const BigInt & p2,BigInt & pr)
{
BigInt q,r,pb;
int test;
test=BIcompare(p1,p2);
if(test>0)
{
r=p1;//Get a number pb that is the largest multiple of(p2)less that (p1)
pb=p2;
MakeThemClose(r,pb);
while(BIcompare(r,pb)>0) pb<<=1;
while(BIcompare(r,pb)<0) pb>>=1;
while((BIcompare(p2,r)<=0)||(BIcompare(p2,pb)<=0))
{
if(BIcompare(r,pb)>=0)
{
r=r-pb;
q<<=1;
q++;
}
else
{
q<<=1;
}
pb>>=1;
}
pr=r;
}
else if(test==0)
{
q++;
pr.setZero();
}
else
{
q.setZero();
pr=p1;
}
return q;
}
//451
//建立预处理表（底数）
void MakeModBaseTable(const BigInt & pp)
{
TableModBase[0]=pp;
for(int i=1;i<36;i++)
{
TableModBase[i]=TableModBase[i-1];
TableModBase[i]<<=2;
}
}
//建立预处理标（乘法积）
void MakeModMultTable(const BigInt & pp)
{
TableModMult[0].setZero();
TableModMult[1]=pp;
for(int i=2;i<16;i++)
{
TableModMult[i]=TableModMult[i-1]+pp;
if(BIcompare(TableModMult[i],TableModBase[0])>=0)
TableModMult[i]=TableModMult[i]-TableModBase[0];
}
}
//两个大整数模乘
BigInt BImodmul(const BigInt & p1,const BigInt & p2,const BigInt & pp)
{
//operator *%
//result=(p1*p2)%pp.
BigInt a,b,c,result;
short i,j,kbit,kbit2;
if(BIisZero(p1))
return result;
if(BIisZero(p2))
return result;
if (TableModBase[0]!=pp)
MakeModBaseTable(pp);
//Suppose p1,p2<pp
a=p1%pp;
b=p2%pp;
kbit2=FindZeroBit(pp.iVal[pp.iIdx]);
for(i=b.iIdx;i>=0;i--)
{
result<<=64;
if(b.iVal[i]!=0)
{
c=a*b.iVal[i];
result=result+c;
}
if (BIcompare(result,pp)>=0)
{
kbit=kbit2-FindZeroBit(result.iVal[result.iIdx]);
if(result.iIdx>pp.iIdx) kbit+=64;
kbit>>=1;
for(j=kbit;j>=0;j--)
{
while(BIcompare(result,TableModBase[j])>=0)
result=result-TableModBase[j];
}
}
}
return result;
}
//p452
//两个大整数模乘（预处理方法）
BigInt BImodmul2(const BigInt & p1,const BigInt & p2,const BigInt & pb)
{
//operator *%
//result=(p1*p2)%pb
BigInt pre_a[16],a,b,pb2,pb4,pb8,result;
unsigned __int64 ma1,ma2;
int i,j,k,test;
if(BIisZero(p1))
return result;
if(BIisZero(p2))
return result;
a=p1;
b=p2;
test=BIcompare(a,pb);
if(test>0)
a=a%pb;
else if(test==0)
return result;
test=BIcompare(b,pb);
if (test>=0)
b=b%pb;
else if(test==0)
return result;
//pre-calculation,a[0000],a[0001],a[0010],....a[1111]
pre_a[0].setZero();
pre_a[1]=a;
for(i=2;i<16;i++)
{
j=i-1;
pre_a[i]=pre_a[j]+a;
if(BIcompare(pre_a[i],pb)>0)
pre_a[i]=pre_a[i]-pb;
}
pb2=pb;
pb2<<=1;
pb4=pb2;
pb4<<=1;
pb8=pb4;
pb8<<=1;
for(i=b.iIdx;i>=0;i--)
{
ma1=LMASK[4];
for(j=LenZ;j>=0;j--)
{
result<<=4;
ma2=b.iVal[i]&ma1;
if(ma2!=0)
{
k=j<<2;
ma2>>=k;
result=result+pre_a[ma2&0x0f];
}
while(BIcompare(result,pb8)>0)
result=result-pb8;
while(BIcompare(result,pb4)>0)
result=result-pb4;
while(BIcompare(result,pb2)>0)
result=result-pb2;
while(BIcompare(result,pb)>0)
result=result-pb;
ma1>>=4;
}
}
return result;
}
//p454
//大整数模幂
BigInt BImodexp(const BigInt & p1,const BigInt &p2,const BigInt &pb)
{
//operator^%
//result=(p1^p2)%pb;
BigInt a,b,result;
if(BIisZero(p1))
return result;
if(BIisZero(p2))
{
result++;
return result;
}
a=p1%pb;
b=p2;
//b=p2%pb;
result++;
do{
if((b.iVal[0]&RMASK[1])==RMASK[1])
result=BImodmul(a,result,pb);
b>>=1;
if(!BIisZero(b))
a=BImodmul(a,a,pb);
}while(!BIisZero(b));
return result;
}
BigInt BIgcd(const BigInt &p1,const BigInt &p2)
{
//gcd
//result=gcd(p1,p2)
BigInt a,b;
a=p1;
b=p2;
while(!BIisZero(b))
{
a=a%b;
BIswap(a,b);
}
return a;
}
//球逆运算
BigInt BIinverse(const BigInt & p1,const BigInt & p2)
{
//i*p1=1 mod p2
//i=1/p1
//i:range 1,...,n-1
BigInt y,g0,g1,g2,v0,v1,v2;
g0=p2;
g1=p1;
v1++;
while(!BIisZero(g1))
{
y=BIdivision(g0,g1,g2);
g0=g1;
g1=g2;
v2=y*v1;
while(BIcompare(v0,v2)<0)
v0=v0+p2;
v0=v0-v2;
BIswap(v0,v1);
}
return v0;
}
//p455
BigInt BImodpow2(const BigInt &p1,const BigInt &pp,const short &d)
{
//result=p1^(2^n) mod pp;
BigInt result;
short i;
result=p1;
i=d;
while(i>0)
{
result=BImodmul(result,result,pp);
i--;
}
return result;
}
//两个大整数模乘（预处理）
BigInt BImodmulPre(const BigInt &p2,const BigInt &pb)
{
//operator*%
//result=(p1*p2)%pp.
BigInt b,pb2,pb4,pb8,result;
unsigned __int64 ma1,ma2;
int i,j,k;
b=p2;
pb2=pb;
pb2<<=1;
pb4=pb2;
pb4<<=1;
pb8=pb4;
pb8<<=1;
for(i=p2.iIdx;i>=0;i--)
{
ma1=LMASK[4];
for(j=LenZ;j>=0;j--)
{
result<<=4;
ma2=b.iVal[i]&ma1;
if(ma2!=0)
{
k=j<<2;
ma2>>=k;
result=result+pre_base[ma2&0x0f];
}
while(BIcompare(result,pb8)>0)
result=result-pb8;
while(BIcompare(result,pb4)>0)
result=result-pb4;
while(BIcompare(result,pb2)>0)
result=result-pb2;
while(BIcompare(result,pb)>0)
result=result-pb;
ma1>>=4;
}
}
return result;
}
//p456
void MakePreTable(const BigInt& p1,const BigInt& pp)
{
int i;
PreTable[0].setZero();
PreTable[0]++;
PreTable[1]=p1;
//pre-calculation.a[0000],a[0001],a[0010],...,a[1111]
for(i=2;i<16;i++)
PreTable[i]=BImodmul(PreTable[i-1],p1,pp);
}
//p457
//模幂（预处理）
BigInt BImodexpPre2(const BigInt &p1,const BigInt &p2,const BigInt &pb)
{
BigInt a,b,result;
unsigned __int64 tmp;
short iz,Lt,count,kZero;
if(BIisZero(p1))
return result;
if(BIisZero(p2))
{
result++;
return result;
}
a=p1%pb;
b=p2%pb;
result++;
MakePreTable(a,pb);
Lt=b.iIdx;
count=TOTALBITS*(Lt+1);
//find total Zero bit to shift.
//Make sure the first bit must be 1.
kZero=FindZeroBit(b.iVal[Lt]);
b<<=kZero;
count-=kZero;
do{
tmp=b.iVal[Lt]&LMASK[4];
if(count>3)
{
tmp>>=60;
b<<=4;
count-=4;
}
else
{
tmp>>=(TOTALBITS-count);
b<<=count;
count-=count;
}
result=BImodmul(result,PreTable[tmp&0x0f],pb);
iz=0;
if(count!=0)
{
while(((b.iVal[Lt]&LMASK[1])!=LMASK[1])&&(count>0))
{
b<<=1;
count--;
iz++;
}
if(count>3)
iz+=4;
else
iz+=count;
result-BImodpow2(result,pb,iz);
}
}while(count!=0);
return result;
}
//p458
//单元乘法
void BImul(const unsigned __int64 &p1,const unsigned __int64 &p2,unsigned __int64 &high,unsigned __int64 &low)
{
//(aH+aL)*(bH+bL)=aH*bH+aH*bL+aL*bH+aL*bL
// -------HIGH--------????????
// -----------LOW-----????????
unsigned __int64 aH,aL,bH,bL;
unsigned __int64 m,m1,m2;
aH=p1>>32;
aL=p1 & 0xFFFFFFFF;
bH=p2>>32;
bL=p2 & 0xFFFFFFFF;
high=aH*bH;
low=aL*bL;
m1=aH*bL;
m2=aL*bH;
high+=(m1>>32);
high+=(m2>>32);
m=(m1 & 0xFFFFFFFF)+(m2 & 0xFFFFFFFF);
high+=(m>>32);
m=m<<32;
low+=m;
if(low<m) high++;
}
//////////////////////////////////////////////////////////////////////////////
//459
void MakeThemClose(const BigInt & p1,BigInt & p2)
{
short shift,Sa,Sb;
shift=TOTALBITS*(p1.iIdx-p2.iIdx);
if(shift!=0){
Sa=FindZeroBit(p1.iVal[p1.iIdx]);
Sb=FindZeroBit(p2.iVal[p2.iIdx]);
if(Sb>Sa)
shift=shift+Sb-Sa;
else
shift=shift+Sa-Sb;
shift=shift-1;
p2<<=shift;
}
}
short FindZeroBit(const unsigned __int64 & p1)
{
unsigned __int64 cc,t;
short Fbit=3;
t=p1;
cc=FULL>>4;
while(t<cc)
{
cc>>=4;
Fbit+=4;
}
cc=FULL>>Fbit;
while((cc&t)!=t)
{
Fbit--;
cc=FULL>>Fbit;
}
return Fbit;
}
//460
//Montgomery乘法
BigInt MontMult(BigInt & p1,BigInt & p2,BigInt & pp)
{
BigInt m,base,result;
BigInt pre_b[16],pre_p[16];
unsigned short mp,u;
int bblock;
base=0x10;
base.iIdx=0;
//Pre-
pre_b[0].setZero();
pre_b[1]=p2;
pre_p[0].setZero();
pre_p[0]=pp;
for(int i=2;i<16;i++)
{
pre_b[i]=pre_b[i-1]+p2;
pre_p[i]=pre_p[i-1]+pp;
}
m=BIinverse(pp,base);
mp=0x10-m.iVal[0];
bblock=16*pp.iIdx+15;
bblock-=pp.Get4BitBlockZero();
result.setZero();
short zero=0;
for(int j=0;j<=bblock;j++)
{
u=(result.Get4BitBlock(zero)+p1.Get4BitBlock(j)*p2.Get4BitBlock(zero))*mp;
u=u & 0xF;
result=result+pre_b[p1.Get4BitBlock(j)];
result=result+pre_p[u];
result>>=4;
}
if(BIcompare(result,pp)>=0)
result=result-pp;
return result;
}
//461
//Montgomery模幂
BigInt MontExpo(BigInt & p1,BigInt & p2,BigInt & pp)
{
BigInt zb,zR,zxp,result;
int bblock,kbit;
unsigned __int64 m;
short k,l;
bblock=pp.Get4BitBlockZero();
zR=pp;
zR.setUpperBlock(bblock);
zb=BImodmul(zR,zR,pp);
zxp=MontMult(p1,zb,pp);
result=zR%pp;
kbit=63-FindZeroBit(p2.iVal[p2.iIdx]);
kbit+=64*p2.iIdx;
for(int j=kbit;j>=0;j--)
{
k=j%64;
l=kbit/64;
m=p2.iVal[l]>>k;
result=MontMult(result,result,pp);
if((m&0x01)==0x01)
result=MontMult(result,zxp,pp);
}
zxp.setZero();
zxp.iVal[0]=1;
result=MontMult(result,zxp,pp);
return result;
}
// 大整数模幂（预处理）
BigInt BImodexpPre(const BigInt & p1,const BigInt & p2,const BigInt & pb)
{
//求p1的p2次方模pp；p1^p2 mod pb
//4 bits block
BigInt result;
//462
// unsigned __int64 tmp;
int bblock,kbit,k,l,b,j,jj;
unsigned m,m2;
if(BIisZero(p1))
return result;
if(BIisZero(p2))
{
result++;
return result;
}
result++;
MakePreTable(p1,pb);
kbit=63-FindZeroBit(p2.iVal[p2.iIdx]);
kbit+=64*p2.iIdx;
for(j=kbit;j>=0;j--)
{
k=j&0x3F;
l=j>>6;
m=p2.iVal[l]>>k;
if((m&0x1)==0x1)
{
m2=m&0x1;
b=1;
jj=j-1;
while((b<4)&(jj>=0))
{
k=jj&0x3F;
l=jj>>6;
m2<<=1;
m=p2.iVal[l]>>k;
m2=m2|(m&0x1);
b++;jj--;
}
j=j-b+1;
while(b>0)
{
result=BImodmul(result,result,pb);
b--;
}
result=BImodmul(result,PreTable[m2&0x0f],pb);
}
else
result=BImodmul(result,result,pb);
}
return result;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -