数论学习笔记

素数:

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
int p[10000000],a[100000001];
int main()
{
    int i,t=0,j,n;
    scanf("%d",&n);
    a[1]=1;
    for (i=2;i<=n;++i)
      {
        if (!a[i])
          p[++t]=i;
        for (j=1;j<=t&&i*p[j]<=n&&(a[i*p[j]]==0||a[i*p[j]]>p[j]);++j)
          {
            a[i*p[j]]=p[j];
            if (!i%p[j]) break;
          }
      }
}

快速幂:

int pow(int x,int a,int p)
{
    int t;
    if (a==0) return 1%p;
    if (a==1) return x%p;
    t=pow(x,a/2,p);
    if (a%2) return ((t*t)%p*x)%p;
    else return (t*t)%p;
}

欧拉函数:

欧拉函数?(n)是指1~n-1内与n互质的数的个数，通常用公式

?(n)=n?(1?1p1)?(1?1p2)?...?(1?1pk)

void Prime1()
{
     memset(a, 0, n * sizeof(a[0]));
     int num = 0, i, j;
     for(i = 2; i < n; ++i)
       if(!a[i])
         {
           p[num++] = i;
             for(j = i*i; j < n; j +=i)
               a[j] = 1;
         }
}
void Prime2()
{
    memset(a, 0, n*sizeof(a[0]));
    int num = 0, i, j;
    for(i = 2; i < n; ++i)
      {
        if(!(a[i]))
          p[num++] = i;
        for(j = 0; (j<num && i*p[j]<n && (p[j]<=a[i]||a[i]==0)); ++j)
            a[i*p[j]] = p[j];//a[j]可以记录每个数的最小质因子
      }
}

欧拉定理:

设gcd(x,p)=1 (若x<p，则称x为p的剩余系)则有

x?(p)=1(modp)

特别地 当p为质数时 ?(p)=p?1 则有

xp?1=1mod(p)

最小公倍数类：

gcd(a,b)?lcm(a,b)=a?b

欧几里德算法(运算次数最多为fib[n]与fib[n-1],辗转相除):

int gcd(int a,int b)
{
   if (a%b==0) return b;
   else return gcd(b,a%b);
}

高精度(Stein 算法):

gcd(x,y)=gcd(x,y/2)(x?>odd,y?>even)

gcd(x,y)=2?gcd(x/2,y/2)(x,y?>even)

gcd(x,y)=gcd(x,y?x)(x,y?>odd)

拓展欧几里德算法，用于求ax+by=gcd(a,b))的最小解(|x|+|y|最小):

void exgcd(int a,int b,int &d,int &x,int &y)
{
    if (b==0)
      {
        d=a;
        x=1;
        y=0;
      }
    else
      {
        exgcd(b,a%b,d,y,x);
        y-=(a/b)*x;
      }
}

同余类:

线性同余方程:

ax=b(modp)=>ax?py=b=>若gcd(a,p)|b 则有扩展欧几里德算法x=x?b/(gcd(a,p))

特别地当b=1且gcd(a,p)=1时，a与x互为模p意义下的逆元

求逆元:

基于扩展欧几里德：

ax?py=1 求得x，注意此时x并非最小正整数解x=(x+p)%p

int inv(int a,int p)
{
    int d,x,y;
    exgcd(a,p,d,x,y);
    if (d==1) return (x+p)%p;
    else return -1;
}

基于欧拉定理

a?(p)=1(modp)=>a?(p)?1?a=1(modp)=>x=a?(p)?1(modp)

线性同余方程组：

中国剩余定理(模数互质):

x=b[1](modm[1])

x=b[2](modm[2])

…

x=b[3](modm[3])

设M=Πni=1m[i]则有

x=∑i=1n(exgcd(M/m[i]%m[i],m[i])+m[i])%m[i]?b[i]?M/m[i]%M

#include<cstdio>
using namespace std;
int b[100],m[100];
int n,x,y,N;
void exgcd(int a,int b)
{
    int t;
    if (a%b==0)
      {
        x=0;
        y=1;
      }
    else
      {
        exgcd(b,a%b);
        t=x;
        x=y;
        y=t-(a/b)*y;
      }
}
int main()
{
    int i,now,ans;
    scanf("%d",&n); N=1; ans=0;
    for (i=1;i<=n;++i)
      {
        scanf("%d%d",&b[i],&m[i]);
        N=N*m[i];
      }
    for (i=1;i<=n;++i)
      {
        now=N/m[i];
        now=now%m[i];
        exgcd(now,m[i]);
        x=(x%m[i]+m[i])%m[i];
        x=(x*(N/m[i])*b[i])%N;
        ans=(ans+x)%N;
      }
    printf("%d\n",ans);
}

模数不互质:考虑合并方程x=b1(modm1)x=b2(modm2)

用扩展欧几里德求出x?m1+b1=y?m2+b2的解

使b=x?m1+b1,m=lcm(m1,m2)得到方程x=b(modm)

a^x=b (mod p):

BSGS法(p为质数或p,a互质,此时下式中p?1为?(p)):

设m=(√p),记am关于p的逆元为v=ap?1?m

暴力枚举 1≤x≤m 分别模p的值，记作e，此时考虑m+1≤x′≤2m?1

若有ax′=b(modp)则有ei?am=b(modp)=>ei=b?a?m(modp)=>ei=b?v(modp)

对于im≤x≤im+m?1循环i，迭代b即可。

int solve(int a,int b,int p)
{
    int m,v,e=1,i;
    m=(int)(sqrt(p));
    v=inv(pow(a,m,p),p); /*基于扩展欧几里德*/
    v=pow(a,n-1-m,p);  /*基于费马小定理*/
    map <int,int> x; x[1]=0;
    for (i=1;i<=m;++i)
      {
        e=(e*a)%p;
        if (!x.count(e)) x[e]=i;
      }
    for (i=0;i<m;++i)
      {
        if (x.count(b)) return i*m+x[b];
        b=(b*v)%n;
      }
    return -1;
}

扩展BSGS法(a,p不互质):

设g=gcd(gcd(a,b),p),原方程写为(a′g)x=b′g(modp′g)=>(a′)x?gx=b′g(modp′g)=>约去g得(a′)x?gx?1=b′(modp′),则用BSGS法即可

组合数及组合数取模类:

Cmnn!/m!(n?m)!

Cmn=Cm?1n?1+Cmn?1

Cmn=(n/m)?Cm?1n?1=Cm?1n?((n?m+1)/m)

组合数取模

当n,m较小p为质数:直接上逆元

当n,m较大p为质数:Lucas定理

Lucas定理:

n=nk?pk+nk?1?pk?1+...+n1?p+n0

m=mk?pk+mk?1?pk?1+...+m1?p+m0 则:

Cmn=Πki=0Cmini(modp)

再上逆元 (mi!(ni?mi)!)?1=(mi!(ni?mi)!)p?2

m次方求和公式：

设

Sm(n)=∑k=1nkm=1m+2m+3m+...+nm

普通求法:

nm+1?(n?1)m+1=anm+bnm?1+cnm?2+...+n

(n?1)m+1?(n?2)m+1=a(n?1)m+b(n?1)m?1+c(n?1)m?2+...+n?1

(n?2)m+1?(n?3)m+1=a(n?2)m+b(n?2)m?1+c(n?3)m?2+...+n?2

…

1m+1?0m+1=a1m+b1m?1+c1m?2+...+1

方程左右同时累加,再移项:

Sm(n)=nm+1?b?Sm?1(n)?c?Sm?2(n)?...?S1(n)a

伯努利数法:

Sm(n)=1m+1∑k=0mCkm+1Bknm+1?k

伯努利数：易得出Sm(0)=0易得出递推公式：

B0=1

Bn=?1n+1∑k=0n?1Ckn+1Bk

原根类:

定义:g1,g2,...,g?(p)能遍历p的剩余系,则称g为p的原根

原根数目为?(?(n))

原根性质 若ga=gb(modp)=>a=b(mod?(p))

如何求原根，枚举判断m，是否有d|phi(p)且md=1(modp)则m不是原根

证明，若m是原根，遍历序列中一定以1结尾(欧拉定理)，若其中存在另外一个1，那么原循环节一定是新循环节的整数倍，则d|phi(p),即只验证phi(p)约数即可。

利用原根求xa=b(modp)的解:

设p的原根为g,x=gy,b=gz,利用BSGS求出z

则有gay=gz(modp)=>ay=z(mod?(p)),利用扩展欧几里德算法，求出y

在用快速幂求得x即可

高斯消元：

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
using namespace std;
double a[1000][1000],b[1000],x[1000];
int n;
int main()
{
    int i,j,k;
    double xi,bz;
    scanf("%d",&n);
    for (i=1;i<=n;++i)
      {
        for (j=1;j<=n;++j)
          {
            scanf("%lf",&xi);
            a[i][j]=xi;
          }
        scanf("%lf",&b[i]);
      }
    for (k=1;k<=n;++k)
      {
        if (a[k][k]==0)
          {
            for (j=k+1;j<=n;++j)
              if (a[j][k]!=0)
                {
                  for (i=1;i<=n;++i)
                    swap(a[j][i],a[k][i]);
                  swap(b[j],b[k]);
                  break;
                }
          }
        for (j=k+1;j<=n;++j)
          {
            bz=a[j][k]/a[k][k];
            for (i=1;i<=n;++i)
              a[j][i]-=a[k][i]*bz;
            b[j]-=b[k]*bz;
          }
      }
    for (k=n;k>=1;--k)
      {
        x[k]=b[k]/a[k][k];
        for (j=k-1;j>=1;--j)
          {
            b[j]-=x[k]*a[j][k];
            a[j][k]=0;
          }
      }
    double a=5-(5/0.0*0.0);
    cout<<a<<endl;
    for (i=1;i<=n-1;++i)
      printf("%.3f ",x[i]);
    printf("%.3f\n",x[n]);
}