数论总结

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~U
voidtwice_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) ;