数学-剩余系

欧几里得算法与扩展欧几里得算法

gcd. 用于求最大公约数.

int gcd(int a,int b)
{ return !a ? b : gcd(b%a,a); }

扩欧. 用于求方程 $ax+by=gcd(a,b)$ 的解$(x,y)$ .

struct value{ int x,y; value(int x,int y):x(x),y(y){} };
value ExtEuclid(int a,int b)
{
if(a==0) return value(0,1);
    value r=ExtEuclid(b%a,a);
    return value(r.y-b/a*r.x,r.x);
}

有解的充要条件:

$$ gcd(a,b)|c $$

它的解为

$$ (x,y)frac{c}{gcd(a,b)} $$

$(x,y)$ 由 $ExEuclid(a,b)$ 的返回值给出.

 

 

 

非递归写法

gcd

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

扩欧的不会.....

 

 

 

 

求剩余系中某数的乘法逆元

考虑剩余系下的除法 $ a/b=ab^{-1} space space (mod space MOD) $ ,有 $ a^{-1} a=1 space space (mod space MOD) $

这个算法求出$a^{-1}$ ........

注意只有模数为质数时,整个剩余系的除零以外其它元素才都会有逆元.

int getrev(int p)
{ return p==0 ? -1 : (ExtEuclid(p,MOD).x%MOD+MOD)%MOD; }

具体的,要求 $ ax equiv 1 space space (mod space MOD) $ 中的x,有 $ ax+MODy=1 $

如果a与MOD的gcd不为1,则逆元不存在.

直接用扩欧来算(x,y)就行了......

 

 

 

 

线性同余方程组

考虑方程组

$$ left{egin{matrix}
{x equiv a_1 ; ( mod ; m_1 ) }\
{x equiv a_2 ; ( mod ; m_2 ) }\
{ ...... }\
{x equiv a_n ; ( mod ; m_n ) }
end{matrix} ight. $$

 其中

$$ forall_{i,j} quad a_i equiv a_j ; (mod ; gcd(m_i,m_j))  quad (有解条件) $$

并且

$$ forall_{i,j} quad gcd(m_i,m_j)=1 quad (中国剩余定理应用条件) $$

 

• 使用中国剩余定理构造解x的方法如下

设 $ M=prod m_i $

设 $ M_i = M/m_i $

设 $ t_i = M_i ^{-1} (mod ; m_i) $ 即 $t_i$ 为 $M_i$ 在模 $m_i$ 剩余系下的逆元.

则 $ x=sum t_i a_i M_i $

AC VIJOS 1164 裸的中国剩余定理.

int n;
ll a[20],m[20];
ll t;

struct value{ ll x,y; value(ll x,ll y):x(x),y(y){} };
value ExEuclid(ll a,ll b)
{
    if(!a) return value(1,1);
    value r=ExEuclid(b%a,a);
    return value(r.y-b/a*r.x,r.x);
}
ll rev(ll k,ll m)
{ return (ExEuclid(k%m,m).x%m+m)%m; }

int main()
{
    n=getint();
    for(int i=0;i<n;i++)
    scanf("%I64d%I64d",m+i,a+i);
    
    t=1;
    for(int i=0;i<n;i++)
    t*=m[i];
    
    ll res=0;
    for(int i=0;i<n;i++)
    res+=t/m[i]*rev(t/m[i],m[i])*a[i];
    
    printf("%I64d
",res%t);
    return 0;
}

注意 $M$ 的值会非常大,考虑是否需要使用高精度.

 

 

应用的时候不要混淆解的判定条件以及定理应用条件!

考虑方程组

$$ x equiv 0 ; (mod ; 2) \ xequiv 2 ;(mod ; 4) $$

显然有解 $x=2,6,....$

而按照上述解法,有 $M_1=4$ ,因而 $M_1 ; mod ; m_1 = 0 $ ,这样 $M_1$ 就不存在逆元 $t_1$ .

 

• 如何得到一般方程组的解呢?

考虑把方程依次加入.

先看第一个方程 $$ x ; mod ; m_1 = a_1 $$

它的解是 $$ x=a_1+km_1 ; ; , ; k in N quad $$

现有解 $$ x_1 leftarrow a_1+km_1  $$

当我们加入第$i$个方程时,

$$ x_i leftarrow x_{i-1}+k imes lcm(m_1,m_2,...,m_{i-1}) $$

其中的 $k$ 满足 $$ x_{i-1} + k imes lcm(m_1,m_2,...,m_{i-1}) ; mod ; m_i = a_i  $$

亦即 $$ x_{i-1} + k imes lcm(m_1,m_2,...,m_{i-1}) + tm_i = a_i $$

使用扩展欧几里得算出k即可.

 

为什么算法是正确的? lcm的意义是什么?

@iwtwiioi 指出lcm是必要的,这里没有给出lcm的必要性证明.

呃,这个东西丢到周末去做吧.........在我彻底理解那玩意之后......

update(3.29):于是还是没改......

这里有归纳证明: $$ 若 x_i 满足了前i个方程,则由上述算法得出的 x_{i+1} 满足了前i+1个方程. $$

我们有 $$ k imes lcm(m_1,m_2,...,m_i) ; mod ; m_p = 0 $$ 这是因为 $1 leq p leq i$ ,从而 $m_p mid lcm(m_1,m_2,...,m_i)$

对于前$i$个方程中的某个方程 $p$ ,必定有 $$ x_i ; mod ; m_p = a_p $$ 这是因为 $x_i$ 满足前i个方程.

根据上边的式子以及递推过程我们得到

$$ x_{i+1} ; mod ; m_p \ = ( x_i + k imes lcm(m_1,m_2,...,m_i)) ; mod ; m_p \ = x_i ; mod ; m_p + k imes lcm(m_1,m_2,...,m_i) ; mod ; m_p \ = a_p + 0 \ = a_p $$

这样 $x_{i+1}$ 就满足前$i$个方程了.我们还剩第$i+1$个方程. 注意我们是按照第 $i+1$ 个方程的约束求出 $k$ 的,很明显 $k$ 的值会使 $x_{i+1}$ 满足方程 $i+1$.

这样就结束了$=omega=$

 

一题. AC POJ 2891 解一般的线性模方程组.

inline ll modc(ll k,const ll MOD)
{ return (k%MOD+MOD)%MOD; }

int n;
ll t;

struct value{ ll x,y; value(ll x,ll y):x(x),y(y){} };
value ExEuclid(ll a,ll b)
{
    if(!a) return value(1,1);
    value r=ExEuclid(b%a,a);
    return value(r.y-b/a*r.x,r.x);
}

ll rev(ll k,ll m)
{ return modc(ExEuclid(k%m,m).x,m); }

ll gcd(ll a,ll b)
{ while(a) b=b%a,swap(a,b); return b; }

int main()
{
    while(scanf("%d",&n)>0)
    {
        bool ok=1;
        ll a,m,x;
        m=getll();
        a=getll();
        
        x=a;
        ll lcm=m;
        
        for(int i=1;i<n;i++)
        {
            m=getll();
            a=getll();
            
            if(!ok) continue;
            
            ll d=gcd(lcm,m);
            
            if((a-x)%d!=0) { ok=false; continue; }
            
            ll k=ExEuclid(lcm,m).x;
            k=(a-x)/d*k%(m/d);

            x+=k*lcm;
            lcm=m/d*lcm;
            x=x%lcm;
        }
        
        if(!ok) printf("-1
");
        else printf("%lld
",modc(x,lcm));
    }
    
    return 0;
}

注意精度问题. 能取模就取模......

由于是同余恒等式,切记不要搞错符号......

ExEuclid的变量不是可以随便换的......

 

再来一题. AC HDU 1573 也是比较裸的同余方程,要求给定范围内解的个数.

没仔细看题导致WA了一个小时 TAT

X是正整数啊..... 我把0算进去了啊...... TAT

Baby Step Giant Step (BSGS)

考虑有一个完全剩余系:$$a^0,a^1,a^2 ; ... ;,a^{p-1} (mod ; p)$$

我们想求一个 $x$ 使得 $$a^x; equiv ; b ; ; (mod ; p)  $$

那么,考虑把剩余系按原顺序分成 $m=lceil{sqrt{p}} ceil$ 块,每一块包含元素

$$ a^{km},a^{km+1},...,a^{km+m-1} $$

,即

$$ a^{km},a^{km}a^1,...,a^{km}a^{m-1} $$

我们把 $$ ba^{-1},ba^{-2},...,ba^{-m+1} $$ 放入一个双射的表(hash,或者map).

然后枚举 $k$ 计算 $ a^{km} $ ,在表里面找到与之相等的元素 $ba^{-i}$ ,那么 $km+i$ 就是答案了.......

嗯,这么做只是因为有如下等式 $$ a^{km} ; equiv ; ba^{-i} ;;(mod; p) $$这是显然的吧.....

AC BZOJ 2242

涵盖了常用的数值算法....

task1: 快速幂

task2: 扩欧

task3: 分块(Baby Step Giant Step)

#include <cstdio>
#include <fstream>
#include <iostream>
 
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>
 
#include <queue>
#include <vector>
#include <map>
#include <set>
#include <stack>
#include <list>
 
typedef unsigned int uint;
typedef long long int ll;
typedef unsigned long long int ull;
typedef double db;
 
using namespace std;
 
inline int getint()
{
    int res=0;
    char c=getchar();
    bool mi=false;
    while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
    while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
    return mi ? -res : res;
}
inline ll getll()
{
    ll res=0;
    char c=getchar();
    bool mi=false;
    while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
    while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
    return mi ? -res : res;
}

//==============================================================================
//==============================================================================
//==============================================================================
//==============================================================================

ll MOD;

struct pl{ ll x,y; pl(ll a,ll b):x(a),y(b){} };
pl ExEuclid(ll a,ll b)
{
    if(!a) return pl(1,1);
    pl r=ExEuclid(b%a,a);
    return pl(r.y-b/a*r.x,r.x);
}

ll FastMulti(ll a,ll x)
{
    ll res=1;
    while(x)
    {
        if(x&1) (res*=a)%=MOD;
        (a*=a)%=MOD;
        x>>=1;
    }
    return res;
}

ll gcd(ll a,ll b)
{ if(a>b)swap(a,b); while(a)b=b%a,swap(a,b); return b; }

ll rev(ll a)
{ return (ExEuclid(a,MOD).x%MOD+MOD)%MOD; }

map<ll,int> g;

int main()
{
    int T=getint();
    int type=getint();
    if(type==1) while(T--)
    {
        ll a=getint();
        ll b=getint();
        MOD=getint();
        printf("%lld
",FastMulti(a,b));
    }
    else if(type==2) while(T--)
    {
        ll a=getint();
        ll c=getint();
        MOD=getint();
        if(c%gcd(a,MOD)!=0) printf("Orz, I cannot find x!
");
        else 
        printf("%lld
",((c/gcd(a,MOD)*ExEuclid(a,MOD).x)%MOD+MOD)%MOD);
    }
    else if(type==3) while(T--)
    {
        g.clear();
        ll a=getint();
        ll b=getint();
        MOD=getint();
        ll lim=(ll)(ceil(sqrt((long db)MOD)));
        ll c=1;
        for(int i=0;i<lim;i++)
        g[(rev(c)*b%MOD+MOD)%MOD]=i,(c*=a)%=MOD;
        ll d=1;
        bool found=false;
        ll res=0;
        for(int i=0;i<=lim;i++)
        {
            map<ll,int>::iterator v=g.find(d);
            if(v!=g.end()) 
            { res=(lim*i+v->second)%MOD; found=true; break; } 
            else (d*=c)%=MOD;
        }
        if(!found) printf("Orz, I cannot find x!
");
        else printf("%lld
",(res%MOD+MOD)%MOD);
    }
    
    return 0;
}

 

原根

考虑一个模 $p$ 剩余系. 我们知道对于 $leq a < p-1$ , $a^1,a^2,...,a^{phi(p)}$ 是 $a$ 在剩余系下的一个循环节.但不一定是最短的.

我们把这个剩余系下循环节长度刚好是 $phi(p)$ 的数 $r$ 称为这个数的原根.

模 $p$ 剩余系的充要条件: $p=2$ 或 $p=4$ 或 $p=s^k$ 或 $p=2s$ ,其中 $s$ 是奇素数, $k$ 是任意正整数.

如何判断 $x$ 是否是原根? 注意 $x^{phi(p)}=p$ 是必然的,不论 $x$ 是什么值.

于是乎,直接枚举 $d|phi(p)$ ,判断 $x^d=p$ 是否成立.如果有除了 $phi(d)$ 以外的数成立,那么这个数就不是原根.否则就是.

如何求原根? 从 $2$ 开始逐一枚举然后判断....... 原根一般都很小.....

 

 

 

未完待续

模意义下的高斯消元

 

 

---

吐槽

 写LaTeX要抓狂了啊啊啊啊啊 $QAQ$

谁告诉我公式左对齐怎么搞Wiki大法好.....

---

 

参考及拓展

1. EN-WIKI上的CRT条目