ACM数学模板.docx
《ACM数学模板.docx》由会员分享，可在线阅读，更多相关《ACM数学模板.docx（20页珍藏版）》请在冰豆网上搜索。
ACM数学模板
2.扩展欧几里得（求ax+by=gcd（a,b）的特解）
voide_gcd（LLa,LLb,LL&d,LL&x,LL&y）{
if（b==0）{
x=1;y=0;d=a;
}
else{
e_gcd（b,a%b,d,y,x）;
y-=x*（a/b）;
}
}
3.中国剩余定理
同余方程组
x ≡a1（modm1）
x ≡a2（modm2）
......
x ≡ak（modmk）
方程组所有的解的集合就是：
x1=N1*a1+N2*a2+...+Nk*ak
其中Nimodmi=1，Ni=ai*ti,可用欧几里得扩展定理求ti.其中M=m1*m2*m3····*mn;
//互质版
#include
usingnamespacestd;
//参数可为负数的扩展欧几里德定理
voidexOJLD（inta,intb,int&x,int&y）{
//根据欧几里德定理
if（b==0）{//任意数与0的最大公约数为其本身。
x=1;
y=0;
}else{
intx1,y1;
exOJLD（b,a%b,x1,y1）;
if（a*b<0）{//异号取反
x=-y1;
y=a/b*y1-x1;
}else{//同号
x=y1;
y=x1-a/b*y1;
}
}
}
//剩余定理
intcalSYDL（inta[],intm[],intk）{
intN[k];//这个可以删除
intmm=1;//最小公倍数
intresult=0;
for（inti=0;imm*=m[i];
}
for（intj=0;jintL,J;
exOJLD（mm/m[j],-m[j],L,J）;
N[j]=m[j]*J+1;//1
N[j]=mm/m[j]*L;//2【注】1和2这两个值应该是相等的。
result+=N[j]*a[j];
}
return（result%mm+mm）%mm;//落在（0,mm）之间，这么写是为了防止result初始为负数，本例中不可能为负可以直接写成：
returnresult%mm;即可。
}
intmain（）{
inta[3]={2,3,2};
intm[3]={3,5,7};
cout<<"结果:
"<}
//不互质版
/**
中国剩余定理（不互质）
*/
#include
#include
#include
usingnamespacestd;
typedeflonglongLL;
LLMod;
LLgcd（LLa,LLb）
{
if（b==0）
returna;
returngcd（b,a%b）;
}
LLExtend_Euclid（LLa,LLb,LL&x,LL&y）
{
if（b==0）
{
x=1,y=0;
returna;
}
LLd=Extend_Euclid（b,a%b,x,y）;
LLt=x;
x=y;
y=t-a/b*y;
returnd;
}
//a在模n乘法下的逆元，没有则返回-1
LLinv（LLa,LLn）
{
LLx,y;
LLt=Extend_Euclid（a,n,x,y）;
if（t!
=1）
return-1;
return（x%n+n）%n;
}
//将两个方程合并为一个
boolmerge（LLa1,LLn1,LLa2,LLn2,LL&a3,LL&n3）
{
LLd=gcd（n1,n2）;
LLc=a2-a1;
if（c%d）
returnfalse;
c=（c%n2+n2）%n2;
c/=d;
n1/=d;
n2/=d;
c*=inv（n1,n2）;
c%=n2;
c*=n1*d;
c+=a1;
n3=n1*n2*d;
a3=（c%n3+n3）%n3;
returntrue;
}
//求模线性方程组x=ai（modni）,ni可以不互质
LLChina_Reminder2（intlen,LL*a,LL*n）
{
LLa1=a[0],n1=n[0];
LLa2,n2;
for（inti=1;i {
LLaa,nn;
a2=a[i],n2=n[i];
if（!
merge（a1,n1,a2,n2,aa,nn））
return-1;
a1=aa;
n1=nn;
}
Mod=n1;
return（a1%n1+n1）%n1;
}
LLa[1000],b[1000];
intmain（）
{
inti;
intk;
while（scanf（"%d",&k）!
=EOF）
{
for（i=0;i scanf（"%I64d%I64d",&a[i],&b[i]）;
printf（"%I64d\n",China_Reminder2（k,b,a））;
}
return0;
}
4.欧拉函数（求一个数前面的所有与这个数互质的数的个数）
Euler函数表达通式：
euler（x）=x（1-1/p1）（1-1/p2）（1-1/p3）（1-1/p4）…（1-1/pn）,其中p1,p2……pn为x的所有素因数，x是不为0的整数。
euler
（1）=1（唯一和1互质的数就是1本身）。
Euler函数有几个性质：
1.如果q,p互质，则Euler（p*q）=Euler（p）*Euler（q）；
2.如果a=p^k，则Euler（a）=p^k-p^k-1;
//直接求解欧拉函数
inteuler（intn）{//返回euler（n）
intres=n,a=n;
for（inti=2;i*i<=a;i++）{
if（a%i==0）{
res=res/i*（i-1）;//先进行除法是为了防止中间数据的溢出
while（a%i==0）a/=i;
}
}
if（a>1）res=res/a*（a-1）;
returnres;
}
//线性筛选欧拉函数O（n）用到了一下性质：
//
（1）若（N%a==0&&（N/a）%a==0）则有:
E（N）=E（N/a）*a;
//
（2）若（N%a==0&&（N/a）%a!
=0）则有:
E（N）=E（N/a）*（a-1）;
//注意：
如果范围过大可能不适宜开数组来做
inteuler[maxN],vis[maxN],prime[maxN/5],e[maxN],cnt=0;
voidmake_euler（）{
memset（vis,0,sizeof（vis））;
euler[1]=1;
for（inti=2;i if（vis[i]==0）{
prime[cnt++]=i;
euler[i]=i-1;
}
for（intj=0;j vis[i*prime[j]]=1;
if（i%prime[j]==0）{
euler[i*prime[j]]=euler[i]*prime[j];
break;
}
elseeuler[i*prime[j]]=euler[i]*（prime[j]-1）;
}
}
}
5.求N以前N的约数个数
约数个数的性质，对于一个数N，N=p1^a1+p2^a2+...+pn^an。
其中p1,p2,p3...pn是N的质因数，a1,a2,a2,...an为相应的指数，则
div_num[N]=（p1+1）*（p2+1）*（p3+1）*...*（pn+1）;
结合这个算法的特点，在程序中如下运用：
对于div_num：
（1）如果i|prime[j]那么div_num[i*prime[j]]=div_sum[i]/（e[i]+1）*（e[i]+2） //最小素因子次数加1
（2）否则div_num[i*prime[j]]=div_num[i]*div_num[prime[j]] //满足积性函数条件
对于e：
（1）如果i|pr[j] e[i*pr[j]]=e[i]+1;//最小素因子次数加1
（2）否则e[i*pr[j]]=1; //pr[j]为1次
#include
#include
#defineM100000
usingnamespacestd;
intprime[M/3],e[M],div_num[M];//e[i]表示第i个素数因子的个数
boolflag[M];
voidget_prime（）
{
inti,j,k;
memset（flag,false,sizeof（flag））;
k=0;
for（i=2;iif（!
flag[i]）{
prime[k++]=i;
e[i]=1;
div_num[i]=2;//素数的约数个数为2
}
for（j=0;jflag[i*prime[j]]=true;
if（i%prime[j]==0）{
div_num[i*prime[j]]=div_num[i]/（e[i]+1）*（e[i]+2）;
e[i*prime[j]]=e[i]+1;
break;
}
else{
div_num[i*prime[j]]=div_num[i]*div_num[prime[j]];
e[i*prime[j]]=1;
}
}
}
}
6.莫比乌斯函数
一个讲得比较清楚的PPT：
线性筛打表莫比乌斯函数：
intmob[maxN],vis[maxN],prime[maxN],cnt=0;
voidmake_mobius（）{
mob[1]=1;
memset（vis,0,sizeof（vis））;
for（inti=2;iif（!
vis[i]）{
mob[i]=-1;
prime[cnt++]=i;
}
for（intj=0;jvis[i*prime[j]]=1;
if（i%prime[j]==0）{
mob[i*prime[j]]=0;
break;
}
elsemob[i*prime[j]]=-mob[i];
}
}
}
7.卢卡斯定理
//卢卡斯定理C（m,n）%p
LLLucas（LLm,LLn,LLp）{
LLres=1;
while（n&&m）{
LLn1=n%p,m1=m%p;
//费马小定理求逆元
res=res*fac[n1]*qsm（fac[m1]*fac[n1-m1]%p,p-2,p）%p;
//扩展欧几里得求逆元
//res=res*fac[n1]*reverse（fac[m1],p）*reverse（fac[n1-m1],p）%p;
n/=p;
m/=p;
}
return（res%p+p）%p;
}
8.卡特兰数
递推公式
9.伯努利数
伯努利数满足条件
，且有
那么继续得到
voidinit_ber（）{
ber[0]=1;
for（inti=1;iLLans=0;
for（intj=0;jans=（ans+comb[i+1][j]*ber[j]）%mod;
ans=-ans*inv（i+1,mod）%mod;
ber[i]=（ans%mod+mod）%mod;
}
}
10.乘法逆元
因为
可能会很大，超过int范围，所以在快速幂时要二分乘法。
逆元打表（O（n））：
inv[1] = 1;
for（int i=2;i {
if（i >= MOD） break;
inv[i] = （MOD - MOD / i） * inv[MOD % i]% MOD;
}
11.millerrabin算法pollardrho算法（概率高效判断素数，求因子）
#include
#include
#include
#include
#include
constintTimes=10;
constintN=5500;
usingnamespacestd;
typedeflonglongLL;
LLct,cnt,c,x,y;
LLfac[N],num[N],arr[N];
LLgcd（LLa,LLb）
{
returnb?
gcd（b,a%b）:
a;
}
LLmulti（LLa,LLb,LLm）
{
LLans=0;
while（b）
{
if（b&1）
{
ans=（ans+a）%m;
b--;
}
b>>=1;
a=（a+a）%m;
}
returnans;
}
LLquick_mod（LLa,LLb,LLm）
{
LLans=1;
a%=m;
while（b）
{
if（b&1）
{
ans=multi（ans,a,m）;
b--;
}
b>>=1;
a=multi（a,a,m）;
}
returnans;
}
boolMiller_Rabin（LLn）
{
if（n==2）returntrue;
if（n<2||!
（n&1））returnfalse;
LLa,m=n-1,x,y;
intk=0;
while（（m&1）==0）
{
k++;
m>>=1;
}
for（inti=0;i{
a=rand（）%（n-1）+1;
x=quick_mod（a,m,n）;
for（intj=0;j{
y=multi（x,x,n）;
if（y==1&&x!
=1&&x!
=n-1）returnfalse;
x=y;
}
if（y!
=1）returnfalse;
}
returntrue;
}
LLPollard_rho（LLn,LLc）
{
LLx,y,d,i=1,k=2;
y=x=rand（）%（n-1）+1;
while（true）
{
i++;
x=（multi（x,x,n）+c）%n;
d=gcd（（y-x+n）%n,n）;
if（1if（y==x）returnn;
if（i==k）
{
y=x;
k<<=1;
}
}
}
voidfind（LLn,intc）
{
if（n==1）return;
if（Miller_Rabin（n））
{
fac[ct++]=n;
return;
}
LLp=n;
LLk=c;
while（p>=n）p=Pollard_rho（p,c--）;
find（p,k）;
find（n/p,k）;
}
voiddfs（LLdept,LLproduct=1）
{
if（dept==cnt）
{
arr[c++]=product;
return;
}
for（inti=0;i<=num[dept];i++）
{
dfs（dept+1,product）;
product*=fac[dept];
}
}
voidSolve（LLn）
{
ct=0;
find（n,120）;
sort（fac,fac+ct）;
num[0]=1;
intk=1;
for（inti=1;i{
if（fac[i]==fac[i-1]）
++num[k-1];
else
{
num[k]=1;
fac[k++]=fac[i];
}
}
cnt=k;
}
constintM=1000005;
boolprime[M];
intp[M];
intk1;
voidisprime（）
{
k1=0;
inti,j;
memset（prime,true,sizeof（prime））;
for（i=2;i{
if（prime[i]）
{
p[k1++]=i;
for（j=i+i;j{
prime[j]=false;
}
}
}
}
intmain（）
{
LLn,t,record,tmp;
isprime（）;
while（cin>>n）
{
LLans=-1;
record=0;
while（true）
{
tmp=1;
if（Miller_Rabin（n））
{
Solve（n-1）;
c=0;
dfs（0,1）;
sort（arr,arr+c）;
boolflag=0;
for（inti=0;i{
if（quick_mod（10,arr[i],n）==1）
{
tmp=arr[i];
break;
}
if（i==c-2）flag=1;
}
if（flag）
{
if（anscout<break;
}
}
else
{
boolflag=false;
LLtmp1=（LL）sqrt（n*1.0）;
if（tmp1*tmp1==n&&Miller_Rabin（tmp1））
{
x=tmp1;
y=2;
flag=true;
}
else
{
LL
升级会员 