基础数论 I:同余

update:

1. 数学概念

1.1 模意义下的逆元

ax1(modp) 的解。

该同余方程有解当且仅当 ap

对于 p 为质数,有 ap2a1(modp)

2.1 中将讲如何求解。

1.2 费马小定理与欧拉定理

1.2.1 费马小定理

无论学什么应该都是最先知道的定理 : )。

费马小定理: 任意 p质数,且 ap,则有 ap11(modp)

引理:对于一个模 p 的完全剩余系 r,若有 aZ,且 pa,则 ri×a(i[0,n1]) 也为模 p 的一个完全剩余系。

  • 证明:反证,假设存在 ij 使得 ri×arj×a(modp),因为 ap 即存在 a 的逆元 a1,既得 rirj(modp),矛盾。

首先有一组模 p 的完系 r,去掉 0r+,共有 p1 个元素,若有 ap,则 (r1×a)(r2×a)...(ap1×a)r1×r2×...×rp1(modp),即得 ap11(modp)

1.2.2 欧拉定理

欧拉定理: 任意 pZ,且 ap,则有 aφ(p)1(modp)

证明类似。

引理:对于一个模 p 的简化剩余系 r,若有 aZ,且 pa,则 ri×a(i[0,n1]) 也为模 p 的一个缩系。

  • 首先这 φ(p) 个数彼此不同余,可由 1.2.1 中证出,下证 ri×ap 互质,因为 rip,且 apri×amodp 也互质。

一个缩系 rap,可以得出 r1r2...rφ(p)(ar1)...(arφ(p))aφ(p)r1...rφ(p)(modp),即得 aφ(p)1(modp)

p 为质数,则 φ(p)=p1 ,即得费马小定理

1.3 阶与原根

非常厉害的东西。

1.3.1 阶

阶的定义:使得 ax1(modm) 的最小正整数 x 被称为 am 的阶,记作 δm(a)

amδm(a) 存在的充要条件。

性质 1: ax1(modm) 当且仅当 δm(a)|x

性质 2: 任意 i[1,δm(a)]aim 两两不同余。

  • 证明:考虑反证,我们假设存在两个数 i,j[1,δm(a)],使得 aiaj(modm),则有 a|ij|1(modm),显然 |ij|<δm(a),这与的最小性矛盾。

性质 3:δm(ak)=δm(a)gcd(δm(a),k)=lcm(δm(a),k)k 具体不会证,可以看初等数论上的充要条件。

  • 注意到 (ak)x=akx,因此 δm(a)|kx,得到 xmin=δm(a)gcd(δm(a),k)

利用 BSGS 可求,复杂度 O(p)

1.3.2 原根

定义:对于 mN+,aZam,若 δm(a)=φ(m),则称 a 为模 m 的原根。

原根判定定理:对于 m3,ama 是模 m 的原根的充要条件是对于任意 φ(m) 的质因子 p,均有 aφ(m)p1(modm)

  • 证明:因为 φ(m)p 的所有质因子取遍了 φ(m)真因数,所以若存在 ad1(modm),使得 d<φ(m),则必然存在 φ(m) 的质因子 p 使得 d|φ(m)p,这说明 aφ(m)p1(modm),矛盾。

原根存在定理:一个数 m 存在原根的充要条件是 m=2,4,pa,2pap 是奇质数)。
证明详见初等数论其实没必要

原根个数: 一个数 m 有原根,则其原根的个数为 φ(φ(m))

  • 证明:由原根的性质知:

δm(gk)=δm(a)gcd(δm(a),k)=φ(m)gcd(φ(m),k)

则若 φ(m)k,则有 δm(gk)=φ(m),即 gk 也是模 m 的原根,而满足 φ(m)kk 的个数为 φ(φ(m)),即原根有 φ(φ(m)) 个。

模版

中国数学家王元证明了最小原根的级别为 O(n0.25+ϵ)
所以我们可以直接枚举,复杂度 O(n0.25log(n))
进而求出所有原根。


2. 扩展欧几里得算法

也叫 exgcd

即求方程 ax+by=c 的解,当且仅当 gcd(a,b)|c 有无穷多项解。

我们考虑欧几里得法的最后一步,即 b=0 时,一定有解 x=1,y=0 ,我们可以归纳求解。

假设已知 gcd(b,amodb)=x0b+(amodb)y0,则我们直接拆开 gcd(a,b)=x0b+(amodb)y0=bx0+ay0baby0=ay0+b(x0aby0)

可得一解 x=y0,y=x0aby0,复杂度同欧几里得,O(logn),最后将 x,y 乘上 cgcd(a,b) 即可。

我们令一组特解为 x0,y0
通解即为 x=x0+bgcd(a,b)k,y=y0agcd(a,b)k(kZ)

模版

int t;
ll a,b,c,gc,x,y; 
void exgcd(ll a,ll b,ll &x,ll &y){
    if(b == 0){return x = 1,y = 0,void();}
    exgcd(b,a%b,y,x);
    y -= a / b * x;
} 
int main(){
    t = read();
    while(t--){
        a = read(),b = read(),c = read();
        gc = __gcd(a,b);
        ll dx = b/gc,dy = a/gc;
        if(c % gc != 0){//无解 
            printf("-1\n");
            continue;
        }
        exgcd(a,b,x,y);
        x = x * c / gc,y = y * c / gc;//特解 
        ll l = ceil(double(-x+1)/dx),r = floor((double)(y-1)/dy);// > 0 的区间 
        if(l > r)printf("%lld %lld\n",x + dx*l,y - dy*r);
        else printf("%lld %lld %lld %lld %lld\n",r-l+1,x + dx*l,y - dy*r,x + dx*r,y - dy*l);
    }

    return cerr<<endl<<"Time:"<<clock()<<"ms"<<endl,0;

}

线性同余方程axc(modp)

可以转化为求解不定方程 ax+kp=c,同余方程有解当且仅当 gcd(a,p)|c,解集是 x=x0+pgcd(a,p)k(kZ)

乘法逆元p 不为质数时可用。

3. 扩展中国剩余定理(ExCRT)

不讲中国剩余定理。

求解 n 个线性同余方程组 xai(modbi)(ai<bi) 的最小解 x

首先第一个同余方程的一个特解显然是 a1,我们考虑已经知道前 i 个同余方程组的解如何求出前 i+1 个同余方程组的解。

我们假设前 i 个同余方程组的解为 x0M=lcm(b1,b2,...bi),则显然前 i 个同余方程组的解集x0+kM(kZ),则我们只需考虑同余方程 x0+kMai+1(modbi+1),移项得 kMai+1x(modbi+1),可以用扩展欧几里得算法求出特解 k0,得出 k 的解集 k=k0+bi+1gcd(bi+1,M)ω(ωZ),求得最小 kmin,最后得出前 i+1 方程组最小解x0=x0+kmin×M

当有一次扩欧无解时整个同余方程无解。

模版

会爆 longlong,需要龟速乘

int n;
ll a[N],b[N];
void exgcd(ll a,ll b,ll &x,ll &y){
    if(b == 0){return x = 1,y = 0,void();}
    exgcd(b,a%b,y,x);
    y -= a / b * x;
} 
ll mul(ll x,ll y,ll p){return (x*y-(ll)(x/(long double)p*y+1e-3)*p+p) % p;}
ll lcm(ll x,ll y){return x / __gcd(x,y) * y;}
ll EXCRT(){
    ll x = a[1],M = b[1],k = 0,y = 0;
    F(i,2,n){
        ll gc = __gcd(M,b[i]);
        if((a[i]-x) % gc != 0)return -1;//无解
        exgcd(M,b[i],k,y);
        ll mod = b[i] / gc;
        k = (mul(k,(a[i] - x) / gc,mod) + mod) % mod;
        x = x + k * M;
        M = lcm(M,b[i]);
    }//EXCRT
    return x;
}
int main(){
    n = read();
    F(i,1,n)b[i] = read(),a[i] = read();
    printf("%lld\n",EXCRT());

    return cerr<<endl<<"Time:"<<clock()<<"ms"<<endl,0;

}

4. 离散对数

求解同余方程 axb(modp) 的解集。
该问题还未有 O(polylog n) 求法。

4.1 BSGS

ap 时,可以 O(p) 求解。

我们令 x=i×tj(t=p,j[0,t1])

则方程变为 aitjb(modp),即 (at)ib×aj(modp)

我们枚举所有的 j[0,t1],将 b×ajmodp插入到一个 Hash 表内。

然后我们再枚举所有的 i[0,t],计算 (at)imodp,在 Hash 表内查找是否存在对应的 j

复杂度 O(n)

模版

ll a,b,p;
ll ksm(ll x,ll y,ll mod){
    ll ans = 1;
    while(y){
        if(y & 1)ans = (ans * x) % mod;
        x = (x * x) % mod,y >>= 1;
    }return ans;
}
ll BSGS(ll a,ll b,ll p){
    map<int,int>mp;mp.clear();
    b %= p;ll t = sqrt(p)+1;
    for(int j = 0;j < t;j++,b = b * a % p)mp[b] = j;
    ll aa = ksm(a,t,p);a = 1;
    for(int i = 0;i <= t;i++,a = a * aa % p)
        if(mp.find(a) != mp.end() && i * t - mp[a] >= 0)return i * t - mp[a]; 
    return -1;
}
int main(){
    p = read(),a = read(),b = read();
    ll ans = BSGS(a,b,p);
    if(ans != -1)printf("%lld\n",ans);
    else printf("no solution\n");

    return cerr<<endl<<"Time:"<<clock()<<"ms"<<endl,0;

}


4.2 ExBSGS

考虑一般情况,我们尽量把问题转化为 ap 的情况,我们设 d=gcd(a,p),我们将等式同除 d,得到 ax1adbd(modpd),若 bd 则无解。

可以发现 adpd,即存在逆元,需要Exgcd求出,得到 ax1bd(ad)1(modpd),我们重复上述操作直到 apdk

最后式子为 axkbdk(adk)1(modpdk),我们可以BSGS求解。

复杂度 O(log2p+p)(需要求逆元)。

模版
模版2

ll a,b,p;
ll ksm(ll x,ll y,ll mod){
    ll ans = 1;
    while(y){
        if(y & 1)ans = (ans * x) % mod;
        x = (x * x) % mod,y >>= 1;
    }return ans;
}
void exgcd(ll a,ll b,ll &x,ll &y){
    if(b == 0)return x = 1,y = 0,void();
    exgcd(b,a%b,y,x);
    y -= a / b * x;
}
ll inv(ll x,ll p){return exgcd(x,p,x,*new ll),(x % p + p) % p;}
ll BSGS(ll a,ll b,ll p){
    map<int,int>mp;mp.clear();
    b %= p;ll t = sqrt(p)+1;//x = i * t - j
    for(int j = 0;j < t;j++,b = b * a % p)mp[b] = j;
    ll aa = ksm(a,t,p);a = 1;
    for(int i = 0;i <= t;i++,a = a * aa % p)
        if(mp.find(a) != mp.end() && i * t - mp[a] >= 0)return i * t - mp[a]; 
    return -1;
}
ll EXBSGS(ll a,ll b,ll p){
    ll d = __gcd(a,p),k = 0;
    b %= p,a %= p;
    if(b == 1 || p == 1)return k;
    while(d > 1){
        if(b == 1)return k;
        if(b % d)return -1;
        p /= d;b = b / d * inv(a / d,p) % p;
        d = __gcd(a % p,p),k++;
    }
    ll ans = BSGS(a,b,p);
    return ans == -1 ? -1 : ans + k;
}
int main(){
    while(1){
        a = read(),p = read(),b = read();
        if(!p)break;
        ll ans = EXBSGS(a,b,p);
        if(ans != -1)printf("%lld\n",ans);
        else printf("No Solution\n");
    }
    

    return cerr<<endl<<"Time:"<<clock()<<"ms"<<endl,0;

}

4.3 习题

I P2485 [SDOI2011] 计算器

第一问快速幂,第二问扩欧,第三问 BSGS

代码

II P3306 [SDOI2013] 随机数生成器

首先我们求一下通项。

xn=axn1+b,我们设 λA 使得 xn+λ=A(xn1+λ)

解得 A=a,λ=ba1.

xn=a(xn1+ba1)ba1=an1(x1+ba1)ba1

我们要求的即是:

an1t+ba1x1+ba1(modp)

这就很板了,BSGS 即可(注意 a=1,a=0 的特判)。

代码

5. 线性方程组

5.1 高斯消元法

求解 n 元一次方程组。

主要就是矩阵变换。
只需要知道将系数矩阵中的任意一行乘上一个常数加到另一行整个方程式不变的。

高斯消元就是把系数矩阵通过矩阵变换转化为对角矩阵进行求解。

模版
模版2

高斯法:

高斯-约旦法:

int n;
db a[N][N],b[N];
void Gauss(){
    for(int i = 1;i <= n;i++){
        int k = i;
        for(int j = i+1;j <= n;j++)
            if(fabs(a[j][i]) > fabs(a[k][i]))k = j;
        for(int j = 1;j <= n;j++)swap(a[i][j],a[k][j]);swap(b[i],b[k]);
        if(fabs(a[i][i]) < 1e-8)return printf("No Solution\n"),void();
        db c = a[i][i];
        for(int j = 1;j <= n;j++)a[i][j] /= c;b[i] /= c;
        for(int j = 1;j <= n;j++){
            if(i == j)continue;
            db c = a[j][i];
            for(int k = i;k <= n;k++)a[j][k] -= c * a[i][k];
            b[j] -= c * b[i];
        } 
    }
    for(int i = 1;i <= n;i++)printf("%.2lf\n",b[i]); 
}
int main(){
    n = read();
    for(int i = 1;i <= n;i++){
        for(int j = 1;j <= n;j++)a[i][j] = read();
        b[i] = read();
    }
    Gauss() ;
    
	return cerr << "Time : " << clock() << " ms" << endl, 0;
        
}

高斯消元只是一个解方程的方法。

5.2 习题

I P4035 [JSOI2008] 球形空间产生器

我们设每个每个点到球心的距离为 C,球心坐标为 (x1,x2,...,xn) 则我们需要解方程:

j=0n(ai,jxj)2=C2

这是一个 n+1 元二次方程,无法解决。
我们把相邻两项相减即可消掉二次项与常数项

j=1n2(ai,jai+1,j)xj=j=1n(ai,j2ai+1,j2)

n 元一次方程组,高斯消元即可。

代码

6. 卢卡斯定理

7. 二次剩余

参考文章:

初等数论学习笔记 I:同余相关 - Alex_wei

posted @   oXUo  阅读(23)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
网站统计
点击右上角即可分享
微信分享提示