数论基础（浅谈数论的部分实现）

最近写到一些基础数论题，
发现一个可怕的事实
基础数论的理论我都懂，但是连最基础的板子都有可能敲错
所以特意停下手中的题，进行基础数论的实现

First.欧几里得(辗转相除)

int gcd(int a,int b)
{
    int r=a%b;
    while (r)
    {
        a=b;b=r;r=a%b;
    }
    return b;
}

Second.扩展欧几里得

ax+by=gcd(a,b)
ax+by=bx’+(a%b)y’
ax+by=bx’+(a-(a/b)*b)y’
ax+by=ay’+b(x’-(a/b)y’)

x=y’
y=x’-(a/b)y’

使用条件

等号右边一定是gcd(a,b)*k (k!=0)
如果是求拟元，则gcd(a,b)==1

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

Third.KSM+费马小定理

使用条件

a,p互质

ll KSM(ll a,int b,ll p)
{
    ll t=1;
    a%=p;
    while (b)
    {
        if (b&1)
           t=(t%p*a%p)%p;
        b>>=1;
        a=(a%p*a%p)%p;
    }
    return t%p;
}

ll fm(ll x,ll p)
{
    return KSM(x,p-2,p);
}

Forth.线性求拟元

a在mod p意义下的拟元
p%a=p-(p/a)*a
p%a=-(p/a)*a
a=(p%a)*(-p/a)^-1
a^-1=inv[p%a]*(p-p/a)

int inv[N];

void INV(int n,int p)
{
    inv[0]=0;
    inv[1]=1;
    for (int i=2;i<=n;i++)
        inv[i]=(p-(p/i))*inv[p%i]%p;
}

Fifth.bsgs

map<ll,int> mp;

int bsgs(ll x,ll z,ll p)
{
    x%=p; z%=p;
    mp.clear();
    if (x==0&&z==0) return 0;
    if (x==0) return -1;
    ll m=(ll)ceil(sqrt((double)p)),now=1;
    mp[1]=m+1;
    for (int i=1;i<m;i++)
    {
        now=(now%p*x%p)%p;
        if (!mp[now]) mp[now]=i;
    }
    ll inv=1,tmp=KSM(x,p-m-1,p);
    for (int k=0;k<m;k++)
    {
        int i=mp[(z%p*inv%p)%p];
        if (i)
        {
            if (i==m+1) i=0;
            return k*m+i;
        }
        inv=(inv%p*tmp%p)%p;
    }
    return -1;
}

Sixth.组合数的递推

int C[N][N];

void doit(int n,int m)
{
    int i,j;
    C[1][1]=1;
    for (i=2;i<=n;i++)
        for (j=1;j<=i;j++)
            C[i][j]=C[i-1][j]+C[i-1][j-1];
    for (i=1;i<=n;i++)
    {
        for (j=1;j<=i;j++)
            printf("%d ",C[i][j]);
        puts("");
    }
}

Seventh.Lucas

使用条件

p模数是素数，
相当于把n，m变成p进制数
C组合数可以预处理
因为p是质数，所以可以用费马小定理求拟元，
当然如果预处理了C就没有这个问题了

int inv(int x,int p)
{
    return KSM(x,p-2,p);
}

int C(int n,int m)
{
    if (m>n) return 0;
    int FZ=1,FM=1;
    for (int i=n-m+1;i<=n;i++) FZ=(FZ*i)%p;
    for (int i=2;i<=m;i++) FM=(FM*i)%p;
    return (FZ%p*inv(FM,p)%p)%p;
}

int Lucas(int n,int m,int p)
{
    if (n<m) return 0;
    int ans=1;
    while (m)
    {
        ans=(ans%p*C(n%p,m%p)%p)%p;
        n/=p;
        m/=p;
    }
    return ans;
}

Eighth.线性筛素数

int sshu[N],tot=0;
bool no[N];

void make(int n)
{
    memset(no,0,sizeof(no));
    for (int i=2;i<=n;i++)
    {
        if (!no[i])
           sshu[++tot]=i;
        for (int j=1;j<=tot&&sshu[j]*i<=n;j++)
        {
            no[sshu[j]*i]=1;
            if (i%sshu[j]==0) break;   //i%sshu[j] 
        }
    }
}

Ninth.欧拉函数（phi）

计算式：
phi(i)=i*∏((j-1)/j) {j是素数且i%j==0}

注意

先除后乘防止炸掉

int phi[N];

void makephi(int n)   //phi[i]小于等于i且与i互质的数的个数 
{
    int i,j;
    for (i=1;i<=n;i++) phi[i]=i;
    for (i=1;i<=tot&&sshu[i]<=n;i++)
        for (j=sshu[i];j<=n;j+=sshu[i])
        {
            phi[j]=phi[j]/sshu[i];
            phi[j]=phi[j]*(sshu[i]-1);
        } 
}

实际上我们是可以把phi的计算放到线性筛中的

void makephi()
{
    phi[1]=1;
    for (int i=2;i<N;i++)
    {
        if (!no[i])
        {
            sshu[++tot]=i;
            phi[i]=i-1;    //素数 
        }
        for (int j=1;j<=tot&&sshu[j]*i<N;j++)
        {
            no[sshu[j]*i]=1;
            if (i%sshu[j]==0)
            {
                phi[i*sshu[j]]=phi[i]*sshu[j];   //积性函数 
                break;
            }
            phi[i*sshu[j]]=phi[i]*phi[sshu[j]];
        }
    }
}

Tenth.莫比乌斯函数（mu）

μ(1)=1;
μ(素数)=-1
μ(分解质因数后，每个质因子<=1个)=-1^（质因子个数）；
μ(其他)=0

莫比乌斯函数完整定义的通俗表达:
1)莫比乌斯函数μ(n)的定义域是N
2)μ(1)=1
3)当n存在平方因子时，μ(n)=0
4)当n是素数或奇数个不同素数之积时，μ(n)=-1
5)当n是偶数个不同素数之积时，μ(n)=1

void makemu(int n)
{
    mu[1]=1;
    memset(no,0,sizeof(no));
    for (int i=1;i<=n;i++)
    {
        if (!no[i])
        {
            sshu[++tot]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=tot&&sshu[j]*i<=n;j++)
        {
            no[sshu[j]*i]=1;
            if (i%sshu[j]==0)
            {
                mu[i*sshu[j]]=0;
                break;
            }
            mu[i*sshu[j]]=-mu[i];
        }
    }
}

Eleventh.中国剩余定理（CRT）

有若干同余方程：
x≡a1 (%m1)
x≡a2 (%m2)
…
（mi互质）
设M=m1*m2*m3*…
对于每一个mi，我们计算M/mi在%mi意义下的逆元k
答案就是ΣM/mi*ki*ai （%M意义下的）

ll china(int n,int *a,int *m)
{
    ll M=1; 
    ll an=0;
    for (int i=1;i<=n;i++) M*=m[i];
    for (int i=1;i<=n;i++)
    {
        ll w=M/m[i];
        exgcd(w,m[i]);
        an+=w*x*a[i]; an%=M;       //x是w%m[i]意义下的逆元 
    }
    return an;
}

Twelfth.约数函数（d）