数论学习笔记

素数:

#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;
}

欧拉函数:

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

$$\phi(n)=n*(1-\frac{1}{p1})*(1-\frac{1}{p2})*...*(1-\frac{1}{pk})$$

m=floor(sqrt(n));
phi[n]=n;
for (int i=2;i<=m;++i)
  if (n%i==0)
    {
      phi[n]=phi[n]/i*(i-1);
      while (n%i==0) n/=i;
    }
if (n>1) phi[n]=phi[n]/n*(n-1);
int phi[maxn];


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

欧拉定理:
设$gcd(x,p)=1$ (若$x<p$，则称$x$为$p$的剩余系)则有

$$x^{\phi(p)}=1 (mod p)$$
特别地 当p为质数时 $\phi(p)=p-1$ 则有
$$x^{p-1}=1 mod (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(mod p)$=$>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^{\phi(p)}=1(mod p)$=>$a^{\phi(p)-1}*a=1(mod p)$=>$x=a^{\phi(p)-1}(mod p)$

线性同余方程组：
中国剩余定理(模数互质):
$x=b[1](mod m[1])$
$x=b[2](mod m[2])$
…
$x=b[3](mod m[3])$
设M=$\Pi_{i=1}^{n} m[i]$则有

$$x=\sum_{i=1}^{n} (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(mod m1) x=b2(mod m2)$
用扩展欧几里德求出$x*m1+b1=y*m2+b2$的解
使$b=x*m1+b1,m=lcm(m1,m2)$得到方程$x=b(mod m)$

a^x=b (mod p):

BSGS法(p为质数或p,a互质,此时下式中$p-1$为$\phi(p)$):
设$m=\sqrt(p)$,记$a^m$关于$p$的逆元为$v=a^{p-1-m}$
暴力枚举 $1\le x\le m$ 分别模p的值，记作e，此时考虑$m+1\le x'\le2m-1$
若有$a^{x'}=b(mod p)$则有$e_i*a^m=b(mod p)$=>$ei=b*a^{-m}(mod p)$=>$ei=b*v(mod p)$
对于$im\le x \le 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(mod p'g)$=>$(a')^x*g^x=b'g(mod p'g)$=>约去g得$(a')^x*g^{x-1}=b'(mod p')$,则用BSGS法即可

组合数及组合数取模类:

$C_{n}^{m}n!/m!(n-m)!$
$C_{n}^{m}=C_{n-1}^{m-1}+C_{n-1}^{m}$
$C_{n}^{m}=(n/m)*C_{n-1}^{m-1}=C_{n}^{m-1}*((n-m+1)/m)$

组合数取模
当$n$,$m$较小$p$为质数:直接上逆元
当$n$,$m$较大$p$为质数:Lucas定理
Lucas定理:
$n=n_k*p^{k}+n_{k-1}*p^{k-1}+...+n_1*p+n0$
$m=m_k*p^{k}+m_{k-1}*p^{k-1}+...+m_1*p+m0$ 则:

$$C_{n}^{m}=\Pi_{i=0}^{k} C_{ni}^{mi} (mod p)$$
再上逆元 $(mi!(ni-mi)!)^{-1}=(mi!(ni-mi)!)^{p-2}$

m次方求和公式：

设

$$S_m(n)=\sum_{k=1}^{n} k^m=1^m+2^m+3^m+...+n^m$$

普通求法:

$$n^{m+1}-(n-1)^{m+1}=an^m+bn^{m-1}+cn^{m-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$$
…
$$1^{m+1}-0^{m+1}=a1^m+b1^{m-1}+c1^{m-2}+...+1$$
方程左右同时累加,再移项:
$$S_m(n)=\frac{n^{m+1}-b*S_{m-1}(n)-c*S_{m-2}(n)-...-S_1(n)}{a}$$

伯努利数法:

$$S_m(n)=\frac{1}{m+1}\sum_{k=0}^{m} C_{m+1}^{k} B_k n^{m+1-k}$$
伯努利数：易得出$S_m(0)=0$易得出递推公式：
$$B_0=1$$
$$B_n=-\frac{1}{n+1} \sum_{k=0}^{n-1} C_{n+1}^{k}B_k$$

原根类:

定义:$g^1,g^2,...,g^{\phi(p)}$能遍历$p$的剩余系,则称$g$为$p$的原根
原根数目为$\phi(\phi(n))$
原根性质 若$g^a=g^b(mod p)$=>$a=b(mod \phi(p))$

如何求原根，枚举判断m，是否有$d|phi(p)且m^d=1(mod p)$则m不是原根
证明，若m是原根，遍历序列中一定以1结尾(欧拉定理)，若其中存在另外一个1，那么原循环节一定是新循环节的整数倍，则$d|\phi(p)$,即只验证$\phi(p)$约数即可。

利用原根求$x^a=b(mod p)$的解:
设$p$的原根为$g$,$x=g^y$,$b=g^z$,利用BSGS求出$z$
则有$g^{ay}=g^{z}(mod p)$=>$ay=z(mod \phi(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]);
}