• 数论


    一.扩展欧几里德(exgcd)

    扩展欧几里德的基本形式是a*x+b*y=gcd(x,y)

    设 a>b。

    1,显然当 b=0,gcd(a,b)=a。此时 x=1,y=0;

    2,ab!=0 时

    设 ax1+by1=gcd(a,b);

    bx2+(a mod b)y2=gcd(b,a mod b);

    根据欧几里德有 gcd(a,b)=gcd(b,a mod b);

    则:ax1+by1=bx2+(a mod b)y2;

    即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2=ay2+b(x2-(a/b)*y2);

    根据恒等定理得:x1=y2; y1=x2-(a/b)*y2;

    这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.

    上面的思想是以递归定义的,因为 gcd 不断的递归求解一定会有个时候 b=0,所以递归可以结束。

    对于一般的不定式 a*x+b*y==c;有整数解得充分必要条件是(c % gcd(a,b)==0)。

    若已经求出 a*x+b*y==gcd(a,b)的解。则原不定式的解是x1=x*(c/gcd(a,b)),y1=y*(c/gcd(a,b))。

    那么原不定式的所有解xi=x1+b/gcd(a,b)*t; yi=y1-a/gcd(a,b)*t;其中t为任意常整数。

    当gcd(a,b)==1时,形象化理解就是当你的x要+b时,你的y要-a,所以,把x变成非负整数的方法是:x=(x%b+b)%b。

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

    二.乘法逆元(模数为p)

    求a^(-1)(mod p意义下)。

    以下仅讨论gcd(a,p)==1的情况,否则不存在逆元。

    在  mod p的意义下我们把x的乘法逆元写作 x^(-1)。

    乘法逆元有如下的性质:x*x^(-1)≡1(mod p)。

    应用:x/y≡x*y^(-1)(mod p)

    1.费马小定理。

    欧拉定理:当gcd(a,p)==1时,a^(φ(p))≡1(mod p)

    a^φ(p)1(mod p)

    a*a^(φ(p)-1)≡1(mod p)

    ∴a^(φ(p)-1)==a^(-1)。

    用快速幂即可出解。

    通常用在p为质数的情况:当p为质数时,φ(p)=p-1。

    a^(p-2)==a^(-1)

    2.exgcd。

    求ax≡1(mod p)的x。

    原方程可化为:ax+py=1。

    由于gcd(a,p)==1,所以用exgcd求出x即可,最后把x弄成非负整数。

    通常用此计算逆元。

    3.线性求逆元(inv[i]表示i的逆元)。

    要求p是质数。我们要求inv[1~p-1]。

    设当前已知inv[1~i-1]。记a=p/i,b=p%i。则p=ai+b。

    ai+b≡0(mod p)

    方程两边同乘i^(-1)*b^(-1)

    (ai+b)*i^(-1)*b^(-1)≡0(mod p)

    a*b^(-1)+i^(-1)≡0(mod p)

    i^(-1)≡-a*b^(-1)(mod p)

    ∴inv[i]=(-p/i*inv[p%i])%p=((p-p/i)*inv[p%i])%p。

    初始化:inv[1]=1,inv[0]=1。

    #include<cstdio>
    long long inv[10000000];
    int main()
    {
        int n,p;scanf("%d%d",&n,&p);inv[0]=inv[1]=1;
        for(int i=2;i<=n;i++)inv[i]=((long long)(p-p/i)*inv[p%i])%(long long)p;
        for(int i=1;i<=n;i++)printf("%d
    ",inv[i]);
        return 0;
    }

    通过这种方法,我们还可以用递归的方式用O(log p)的时间求某个数的逆元。

    即我们要查inv[i],就dfs调用inv[p%i]。

    三.BSGS

    以下仅讨论p为质数且gcd(a,p)==1的情况

    我们要求a^x≡b(mod p)(给出a,b)的最小非负整数的x

    我们先考虑暴力枚举x。那么我们枚举的x应该何时停止?

    由费马小定理得:a^(p1)1(mod p)

    a^x/a^(p1)^m≡a^x(mod p)(m为非负整数)

    a^(x-m*(p-1))≡a^x(mod p)

    ∵x-m*(p-1)==x%(p-1)

    ∴a^(x%(p-1))≡a^x(mod p)

    所以,我们的x只要从0枚举到p即可。

    设x=ki-j。

    则a^(ki-j)≡b(mod p)

    a^ki≡b*a^j(mod p)

    同样的,我们的ki-j只要枚举到p即可。

    那么我们令k=√p(上取整),把i从1枚举到k,j从0枚举到k,这样我们的ki-j就只会枚举到p了。

    可是这样还是O(k^2)=O(p)的啊。

    以下的b*a^j与a^ki均膜p。

    我们先枚举j,把所有枚举出来的b*a^j压入hash表,如果有两个j,b*a^j的值是一样的,那么我们取大的一个j,因为这样ki-j最小。

    然后再枚举i,对于每个a^ki,我们查表,如果查到一个j满足a^ki==b*a^j,那么我们停止搜索,此时的ki-j即为答案。

    时间复杂度(可以使用hash表来代替map,降低时间复杂度,map好写)为O(√plogp)

    typedef long long ll;
    map<ll,ll>mp;
    ll mod;
    ll ksm(ll a,ll p)
    {
        ll ans=1;
        while(p)
        {
            if(p&1)ans=(ans*a)%mod;
            a=(a*a)%mod;p>>=1;
        }
        return ans;
    }
    ll BSGS(ll a,ll b)
    {
        if(a%mod==0)return -1;
        mp.clear();
        ll m=(ll)sqrt((double)mod)+1;ll am=ksm(a,m);
        ll now=1;
        for(int j=0;j<=m;j++){mp[(b*now)%mod]=j;now=(now*a)%mod;}now=1;
        for(int i=1;i<=m;i++){now=(now*am)%mod;if(mp[now])return i*m-mp[now];}
        return -1;
    }

    四.中国剩余定理(CRT)

    我们要求一个数x,满足以下这些同余方程(要求p1,p2,...,pn互质)。

    x≡a1(mod p1)

    x≡a2(mod p2)

    ...

    x≡an(mod pn)

    我们设x1是第一个同余方程的一个解,x2是第二个同余方程的一个解...xn是第n个同余方程的一个解

    当x2~xn都是p1的倍数,x1,x3~xn都是x2的倍数......x1~xn-1都是xn的倍数时

    我们的答案x=x1+x2+...+xn,此时的x可以符合每一个同余方程。

    我们不妨假设此时的a1<p1,a2<p2...an<pn(否则可以mod p)

    x1 mod p1==a1且是p2*p3*...*pn的公倍数。

    x2 mod p2==a2且是p1*p3*...*pn的公倍数。

    ...

    xn mod pn==an且是p1*p2*...*pn-1的公倍数。

    我们知道:如果a mod b==1,那么ac mod b==c(c为正整数且c<b)

    x1 mod p1==1且是p2*p3*...*pn的公倍数。x1*=a1。

    x2 mod p2==1且是p1*p3*...*pn的公倍数。x2*=a2。

    ...

    xn mod pn==1且是p1*p2*...*pn-1的公倍数。xn*=an。

    记M=p1*p2*...*pn。

    显然,当我们的xi是M/pi的倍数时,我们满足了第2个条件。

    我们要找到一个ki,使得M/pi *ki≡1(mod pi),这样就能满足第1个条件。

    显然,ki是M/pi的逆元(mod pi意义下)!

    于是,xi=M/pi * (M/pi)^(-1) * ai。x=x1+x2+...+xn。

    记Mi=M/pi。则x=Σai*Mi*Mi^(-1)。

    这只是一个解,其余解可以表示为Σai*Mi*Mi^(-1)   +kM(k为整数)(因为我们不管是加上还是减去M,所得的x仍会符合以上的同余方程)。

    所以最小非负整数解x'=Σai*Mi*Mi^(-1)   %M

         扩展中国剩余定理(ExCRT)

    考虑合并

    x≡a1(mod p1)

    x≡a2(mod p2)

    x=a1+k1p1=a2+k2p2

    k1p1-k2p2=a2-a1

    考虑Xp1+Yp2=gcd(p1,p2)

    k1=X(a2-a1)/gcd(p1,p2),k2=-Y(a2-a1)/gcd(p1,p2)

    因为x=a1+k1p1 >=0 ,所以X为满足Xp1+Yp2=gcd(p1,p2)最小的>=0的数

    这样a1+k1p1=a2+k2p2

    所以,当x=a1+k1p1=a2+k2p2时可以满足上面两个方程,且x=a1+k1p1 + K*lcm(p1,p2)的x也都可以满足上面两个方程(K为任意自然数)

    所以x≡a1+k1p1(mod lcm(p1,p2))的x,就能满足上面两个方程,且不满足该条件的x,不能满足上面两个方程

  • 相关阅读:
    UltraEdit的配置
    字符编码笔记:ASCII,Unicode和UTF-8
    Hello World Hexo
    好久不见,味素
    记一次springboot+dubbo+zookeeper吐血的问题
    [深度学习]模型部署之优化
    [深度学习]pytorch部署之onnx
    count(*)、count(1)、count(column)
    like %和-的区别与使用
    SQL语言分类
  • 原文地址:https://www.cnblogs.com/lher/p/8135197.html
Copyright © 2020-2023  润新知