数论总结

1.欧几里得 求最大公约数,最小公倍数
(1)递归的写法:

int gcd(int a,int b) {return b?gcd(b,a%b):a;}

(2)最小公倍数: 

int lcm(int a,int b) {return a/gcd(a,b)*b;}

2.扩展欧几里得 求ax=b (mod m)=> ax+my=b 

如果d=gcd(a,m)且b%d==0,

则同余方程有解，其最小解为x=x0*(b/d);
ax+by=c 如d=gcd(a,b),则存在x,y,使xa+yb=d;

当x+=b,y-=a后仍然成立        
因为xa+yb+ab-ab=d;==>(x+b)a+(y-a)b=d

long long ext_gcd(long long a, long long b, long long & x, long long &y) {    if(!b)     {         x=1;         y=0;         return a; //d=a，x=1,y=0,此时等式d=ax+by成立      }     longlong d=ext_gcd(b,a%b,x,y);     long long xt=x;     x=y;     y=xt-a/b*y;     return d;      //对于拓展的欧几里德算法有：x=y0;y=x0-a/b*y0; } 

3.素数判定
(1)试除法:

bool isprime(int n){   int temp=(int)(sqrt(n*1.0)+1);   for(int i=2;i<=temp;i++)       if(n%i==0)              return  false;   return  true;}

(2)miller-rabin 算法
bool witness(i64 a,i64 n)
{

i64 x,d=1,i=ceil(log(n-1.0)/log(2.0))-1;
       for(;i>=0;i--)
       {
              x=d; d=(d*d)%n;
              if(d==1&&x!=1&&x!=n-1) return 1;
              if(((n-1)&(1<<i))>0) d=(d*a)%n;
       }
       return  d==1?0:1;
}

bool miller_rabin(i64 n)
{
  if(n==2)   return 1;
  if(n==1||(n&1)==0) return 0;
  i64 j,a;
  for(j=0;j<50;j++)
  {
   a=rand()*(n-2)/RAND_MAX+1;
   if(witness(a,n)) return 0;
  }
  return 1;
}

 

另一种写法，更好理解

bool witness(i64 a,i64 n)

{

 int i,j=0;

 i64 m=n-1,x,y;

 while(m%2==0)

 {

  m>>=1;

  j++;

 }

 x=pow(a,m,n);///快速幂取模

 for(i=1;i<=j;i++)

 {

  y=pow(x,2,n);

  if(y==1&&x!=1&&x!=n-1) return true;

  x=y;

 }

 return y==1?false:true;

}

bool miller_rabin(i64 n)

{

 if(n==2) return true;

 if(n==1||n%2==0) return false;

 for(int i=1;i<=10;i++)

 {

  i64 a=rand()%(n-1)+1;

  if(witness(a,n)) return false;

 }

 return true;

}

 

4.素数筛法  

(1)  O(n)

void get_prime(){    memset(tag,0,sizeof(tag));    for(int i=2;i<N;i++)    {        if(!tag[i]) p[cnt++]=i;        for(intj=0;j<cnt&&p[j]*i<N;j++)        {            tag[i*p[j]]=1;//tag[x]=1表示x不是素数            if(i%p[j]==0) break;        }    }}

(2) O(nlog(n))

void get_prime()   //打表法求出所有的素数 {     int m=sqrt(SIZE)+1;      memset(yes,0,sizeof(yes));     for(int i=2;i<=m;i++)         if(!yes[i])         {             for(int j=i*i;j<SIZE;j+=i)                 yes[j]=1;         } } 

(3) 素数的二次筛选

/************************************************intyes[SIZE]; *long longprim[SIZE]; *int k=0;     //prim[] 中存放着k个素数*intadp[1001000]  //二次需要筛选的（U-L）< 1001000*************************************************///获得基准素数voidgetPrim() {     memset(yes,0,sizeof(yes));     for(long long i=2;i<SIZE;i++)         if(!yes[i])         {             prim[k++]=i;             for(long longj=i*i;j<SIZE;j+=i)                 yes[j]=1;         } }  //二次筛选范围为：L~Uvoidtwice_choose(){    memset(adp,0,sizeof(adp));     for(int i=0;i<k&&prim[i]*prim[i]<=U;i++)     {         long longbeg=L/prim[i]+(L%prim[i]>0);         if(beg==1) beg++;         for(long longj=beg*prim[i];j<=U;j+=prim[i])  //调整j，使得L<=j<=U             adp[j-L]=1;     }}

 

5.整数分解

 (1)

/********************************************* prime[]数组保存筛选出的素数,长度为N;* p[]数组保存n的素数因子;* t[]保存的素数因子的个数,共tot个*********************************************/int tot=0;void split(intn,int *p,int *t){    for(int i=1;i<=N;i++)    {        int s=0;        while(n%prime[i]==0)        {            s++;            n/=prime[i];        }        if(s) {p[tot]=prime[i]; t[tot]=s;tot++;}        if(n==1) break;        if(n<prime[i]*prime[i])        {            p[tot]=n;            t[tot]=1;            n=1;            break;        }    }} 

(2)Pollard-pho大数分解

i64 Pollard(i64 n,int c)

{

 i64 i=1,k=2,x=rand()%n,y=x,d;

 srand(time(NULL));

 while(true)

 {

  i++;

  x=(mod_mult(x,x,n)+c)%n;

  d=gcd(y-x,n);

  if(d>1&&d<n) return d;

  if(y==x) return n;

  if(i==k)  { y=x; k<<=1; }

 }

}

6. 任意正整数都有且只有一种方式写出其素因子的乘积表达式:

A=(p1^k1)*(p2^k2)*(p3^k3)*....*(pn^kn) ，其中pi均为素数.

   有A的所有因子之和为

   S = (1+p1+p1^2+p1^3+...p1^k1) *(1+p2+p2^2+p2^3+….p2^k2) * (1+p3+ p3^3+…+p3^k3) * .... *(1+pn+pn^2+pn^3+...pn^kn)

(1).用递归二分求等比数列1+pi+pi^2+pi^3+...+pi^n：

1）若n为奇数,一共有偶数项，则：
      1 + p + p^2 + p^3 +...+ p^n

      = (1+p^(n/2+1))+ p * (1+p^(n/2+1)) +...+ p^(n/2) * (1+p^(n/2+1))
      = (1+ p + p^2 +...+ p^(n/2)) * (1 +p^(n/2+1))

上式中颜色加粗的前半部分恰好就是原式的一半，那么只需要不断递归二分求和就可以了，后半部分为幂次式，将在下面第4点讲述计算方法。

2）若n为偶数,一共有奇数项,则:
      1 + p + p^2 + p^3 +...+ p^n

      = (1+p^(n/2+1))+ p * (1+p^(n/2+1)) +...+ p^(n/2-1) * (1+p^(n/2+1)) + p^(n/2)
      = (1+ p + p^2 +...+ p^(n/2-1)) *(1+p^(n/2+1)) + p^(n/2);

   上式颜色加粗的前半部分恰好就是原式的一半，依然递归求解

long longMOD=9901;long longp_sum(int x, long long y) {     long long k1, k2, t;      if(y==0) return 1;     if(y==1) return (1+x);     if(y%2==1)    //求x的最高次数为y的等比数列。当y为奇数，可直接用公式。。。     {        k1=(p_mi(x,y/2+1)+1)%MOD;         k2=p_sum(x,y/2)%MOD;         return (k1*k2)%MOD;     }     else          //y为偶数时，则y-1为奇数，求出x^y+p_sum(x,y-1)     {         k1=p_mi(x,y)%MOD;         k2=p_sum(x,y-1)%MOD;         return (k1+k2)%MOD;     } }

 

//递归求x^y次幂 long longp_mi(int x, long long y) {     long long a;      if(y==0) return 1;     if(y==1) return x;     if(y%2==1)    //y为奇数时有，x^y=((x^((y-1)/2))^2)*x     {         a=p_mi(x,(y-1)/2)%MOD;         return ((a*a%9901)*x)%MOD;     }     else     {         a=p_mi(x,y/2)%MOD;         return (a*a)%MOD;     } }  

Code:

int euler_phi(int n){    int m=(int) sqrt(n+0.5);    int ans=n;       for(int i=2;i<=m;i++)        if(n%i==0)        {            ans=ans/i*(i-1);            while(n%i==0)                n/=i;        }    if(n>1) ans=ans/n*(n-1);} 

2)    求出1～n中所有数的欧拉函数值.

欧拉函数的递推关系：

假设素数p能整除n，那么

如果(n%p==0 && (n/p)%p==0) (p还能整除n / p) , PHI(n) = PHI(n / p) * p;

如果(n%p==0 && (n/p)%p==1) (p不能整除n / p), PHI(n) = PHI(n / p) * (p - 1);

可以用与筛选求素数非常类似的方法，在O(nloglogn)时间内计算完毕。

void phi_table(int n, int * phi){    for(int i=2;i<=n;i++)phi[i]=0;    phi[1]=1;       for(int i=2;i<=n;i++)        if(!phi[i])        {            for(intj=i;j<=n;j+=i)            {                if(!phi[j])phi[j]=j;               phi[j]=phi[j]/i*(i-1);            }        }}

8.求逆元ax=1 (mod m) , x是a的逆元 =>ax+my=1=> gcd(a,m)=1
/*

int mod(int x,int n)

{

    x%=n;   //如x是负数，进行处理

    if(x<0) x+=n;

    return x; //return (x%n+n)%n;会越界

}

*/

用扩展欧几里得求:

int Inv(int a,int m){    int r,x,y;    r=exgcd(a,m,x,y);    if(r==1) return (x%m+m)%m;    return -1;}

9.

(1)幂取模（ a^n mod m）

法一：

int pow_mod(inta, int n, int m){    int x=pow_mod(a,n/2,m);    long long ans=(long long)x*x%m;    if(n%2==1) ans=ans*a%m;    return (int)ans;}

法二：用快速幂取模（a^n mod m ）

int ModPow(inta,int n,int m){    int rec=1;    while(n)    {        if (n & 1)            rec = (rec * a) % m;        a = (a * a) % m;        n >>= 1;    }    return rec % m;}

（2）快速幂（ a^k）

int pow(int a,int k){    int rec = 1;    while( k )    {        if (k & 1)            rec *= a;        a *= a;        k >>= 1;    }    return rec;}

10.快速模乘 a*b%p  (位运算的思路)

int mul(int a, int b, int p){    intr=0;    while(b)    {        if(b&1) r=(r+a)%p;        a=(a<<1)%p;        b>>=1;    }    return r;}

11.求解模线性方程组(中国剩余定理)
  x=a1 mod m1
  x=a2 mod m2
  ......
  x=an mod mn  其中,a[],m[]已知,m[i]>0且m[i]与m[j]互质,求x.

同余模公式：

(a+b)%m=(a%m+b%m)%m

(a*b)%m=(a%m*b%m)%m

设m1,m2,...,mn是两两互素的正数,则对任意的整数a1,a2,...,an,同余方程组
其解为:X=((M_1*M1*a1)+(M_2*M2*a2)+...+(M_n*Mn*an)) mod m;
其中M=m1*m2*...*mn;  Mi=M/mi; M_i是Mi的逆元(mod mi).

int china(int*a,int *m,int n){    int M=1,ans=0,Mi,i,x,y;    for(i=0;i<n;i++) M*=m[i];    for(i=0;i<n;i++)    {        Mi=M/m[i];        exgcd(m[i],Mi,x,y);        ans=(ans+Mi*y*a[i])%M;    }    return (ans%M+M)%M;}

 

12.枚举出一个集合的子集。设原集合为mask，则下面的代码就可以列出它的所有子集：

for (i = mask ; i ; i = (i - 1) & mask) ;