数论学习笔记

2016.2.14
就在情人节的这天，kzoacn给我们讲了数论，懵逼一天。。
把我懂的先做个总结。。
day6下载链接：
http://pan.baidu.com/s/1eQOYRxo
具体概念不介绍了，写点有用的东西。。
线性筛欧拉函数：
欧拉函数ϕ(n)是指1~n-1内与n互质的数的个数，通常用公式：

void phi()
{
    int i,j;
    f[1]=1;
    for (i=2;i<=n;i++)
      {
        if (!flag[i])
          {
            tot++;
            prime[tot]=i;
            f[i]=i-1;
          }
        for (j=1;j<=tot;j++)
          {
            if (i*prime[j]>n)
              break;
            flag[i*prime[j]]=true;
            if (i%prime[j]==0)
              {
                f[i*prime[j]]=f[i]*prime[j];
                break;
              }
            else
              f[i*prime[j]]=f[i]*(prime[j]-1);
          }
      }
}

欧拉定理：

素数：（普通的n²算法肯定不太好用了。。）

#include<cstdio>
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;
          }
      }
    for (i=1;i<=t;i++)
      printf("%d ",p[i]);
    return 0;
}

快速幂：

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

最大公约数：
欧几里得算法：

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意义下的逆元
求逆元：

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

基于拓展欧几里得：
ax-py=1求得x，注意此时x并非最小正整数解x=（x+p）%p

#include<cstdio>
using namespace std;
int b[100],m[100],n,x,y,sum;

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);
    sum=1;
    ans=0;
    for (i=1;i<=n;i++)
      {
        scanf("%d %d",&b[i],&m[i]);
        sum*=m[i];
      }
    for (i=1;i<=n;i++)
      {
        now=sum/m[i];
        now%=m[i];
        exgcd(now,m[i]);
        x=(x%m[i]+m[i])%m[i];
        x=(x*(sum/m[i])*b[i])%sum;
        ans=(ans+x)%sum;
      }
    printf("%d",ans);
    return 0;
}

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

原根类:
定义: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|ϕ(p),即只验证ϕ(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]);
}

快速乘：

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