数学-剩余系

最新推荐文章于 2021-04-24 19:54:46 发布

林尧彬

最新推荐文章于 2021-04-24 19:54:46 发布

阅读量242

点赞数

原文链接：http://www.cnblogs.com/DragoonKiller/p/4355855.html

版权

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

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\{\begin{matrix}
{x \equiv a_1 \; ( mod \; m_1 ) }\\
{x \equiv a_2 \; ( mod \; m_2 ) }\\
{ ...... }\\
{x \equiv a_n \; ( mod \; m_n ) }
\end{matrix}\right. $$

其中

$$ \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 裸的中国剩余定理.

 1 int n;
 2 ll a[20],m[20];
 3 ll t;
 4 
 5 struct value{ ll x,y; value(ll x,ll y):x(x),y(y){} };
 6 value ExEuclid(ll a,ll b)
 7 {
 8     if(!a) return value(1,1);
 9     value r=ExEuclid(b%a,a);
10     return value(r.y-b/a*r.x,r.x);
11 }
12 ll rev(ll k,ll m)
13 { return (ExEuclid(k%m,m).x%m+m)%m; }
14 
15 int main()
16 {
17     n=getint();
18     for(int i=0;i<n;i++)
19     scanf("%I64d%I64d",m+i,a+i);
20     
21     t=1;
22     for(int i=0;i<n;i++)
23     t*=m[i];
24     
25     ll res=0;
26     for(int i=0;i<n;i++)
27     res+=t/m[i]*rev(t/m[i],m[i])*a[i];
28     
29     printf("%I64d\n",res%t);
30     return 0;
31 }

View Code

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

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

考虑方程组

$$ x \equiv 0 \; (mod \; 2) \\ x\equiv 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 \times lcm(m_1,m_2,...,m_{i-1}) $$

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

亦即 $$ x_{i-1} + k \times 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 \times 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 \times lcm(m_1,m_2,...,m_i)) \; mod \; m_p \\ = x_i \; mod \; m_p + k \times 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 解一般的线性模方程组.

 1 inline ll modc(ll k,const ll MOD)
 2 { return (k%MOD+MOD)%MOD; }
 3 
 4 int n;
 5 ll t;
 6 
 7 struct value{ ll x,y; value(ll x,ll y):x(x),y(y){} };
 8 value ExEuclid(ll a,ll b)
 9 {
10     if(!a) return value(1,1);
11     value r=ExEuclid(b%a,a);
12     return value(r.y-b/a*r.x,r.x);
13 }
14 
15 ll rev(ll k,ll m)
16 { return modc(ExEuclid(k%m,m).x,m); }
17 
18 ll gcd(ll a,ll b)
19 { while(a) b=b%a,swap(a,b); return b; }
20 
21 int main()
22 {
23     while(scanf("%d",&n)>0)
24     {
25         bool ok=1;
26         ll a,m,x;
27         m=getll();
28         a=getll();
29         
30         x=a;
31         ll lcm=m;
32         
33         for(int i=1;i<n;i++)
34         {
35             m=getll();
36             a=getll();
37             
38             if(!ok) continue;
39             
40             ll d=gcd(lcm,m);
41             
42             if((a-x)%d!=0) { ok=false; continue; }
43             
44             ll k=ExEuclid(lcm,m).x;
45             k=(a-x)/d*k%(m/d);
46 
47             x+=k*lcm;
48             lcm=m/d*lcm;
49             x=x%lcm;
50         }
51         
52         if(!ok) printf("-1\n");
53         else printf("%lld\n",modc(x,lcm));
54     }
55     
56     return 0;
57 }

View Code

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

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

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

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

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

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

View Code

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}}\rceil$ 块,每一块包含元素

$$ 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)

  1 #include <cstdio>
  2 #include <fstream>
  3 #include <iostream>
  4  
  5 #include <cstdlib>
  6 #include <cstring>
  7 #include <algorithm>
  8 #include <cmath>
  9  
 10 #include <queue>
 11 #include <vector>
 12 #include <map>
 13 #include <set>
 14 #include <stack>
 15 #include <list>
 16  
 17 typedef unsigned int uint;
 18 typedef long long int ll;
 19 typedef unsigned long long int ull;
 20 typedef double db;
 21  
 22 using namespace std;
 23  
 24 inline int getint()
 25 {
 26     int res=0;
 27     char c=getchar();
 28     bool mi=false;
 29     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 30     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 31     return mi ? -res : res;
 32 }
 33 inline ll getll()
 34 {
 35     ll res=0;
 36     char c=getchar();
 37     bool mi=false;
 38     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 39     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 40     return mi ? -res : res;
 41 }
 42 
 43 //==============================================================================
 44 //==============================================================================
 45 //==============================================================================
 46 //==============================================================================
 47 
 48 ll MOD;
 49 
 50 struct pl{ ll x,y; pl(ll a,ll b):x(a),y(b){} };
 51 pl ExEuclid(ll a,ll b)
 52 {
 53     if(!a) return pl(1,1);
 54     pl r=ExEuclid(b%a,a);
 55     return pl(r.y-b/a*r.x,r.x);
 56 }
 57 
 58 ll FastMulti(ll a,ll x)
 59 {
 60     ll res=1;
 61     while(x)
 62     {
 63         if(x&1) (res*=a)%=MOD;
 64         (a*=a)%=MOD;
 65         x>>=1;
 66     }
 67     return res;
 68 }
 69 
 70 ll gcd(ll a,ll b)
 71 { if(a>b)swap(a,b); while(a)b=b%a,swap(a,b); return b; }
 72 
 73 ll rev(ll a)
 74 { return (ExEuclid(a,MOD).x%MOD+MOD)%MOD; }
 75 
 76 map<ll,int> g;
 77 
 78 int main()
 79 {
 80     int T=getint();
 81     int type=getint();
 82     if(type==1) while(T--)
 83     {
 84         ll a=getint();
 85         ll b=getint();
 86         MOD=getint();
 87         printf("%lld\n",FastMulti(a,b));
 88     }
 89     else if(type==2) while(T--)
 90     {
 91         ll a=getint();
 92         ll c=getint();
 93         MOD=getint();
 94         if(c%gcd(a,MOD)!=0) printf("Orz, I cannot find x!\n");
 95         else 
 96         printf("%lld\n",((c/gcd(a,MOD)*ExEuclid(a,MOD).x)%MOD+MOD)%MOD);
 97     }
 98     else if(type==3) while(T--)
 99     {
100         g.clear();
101         ll a=getint();
102         ll b=getint();
103         MOD=getint();
104         ll lim=(ll)(ceil(sqrt((long db)MOD)));
105         ll c=1;
106         for(int i=0;i<lim;i++)
107         g[(rev(c)*b%MOD+MOD)%MOD]=i,(c*=a)%=MOD;
108         ll d=1;
109         bool found=false;
110         ll res=0;
111         for(int i=0;i<=lim;i++)
112         {
113             map<ll,int>::iterator v=g.find(d);
114             if(v!=g.end()) 
115             { res=(lim*i+v->second)%MOD; found=true; break; } 
116             else (d*=c)%=MOD;
117         }
118         if(!found) printf("Orz, I cannot find x!\n");
119         else printf("%lld\n",(res%MOD+MOD)%MOD);
120     }
121     
122     return 0;
123 }

View Code

原根

考虑一个模 $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条目

转载于:https://www.cnblogs.com/DragoonKiller/p/4355855.html

林尧彬

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
数学-剩余系

欧几里得算法与扩展欧几里得算法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 Ex...
复制链接

扫一扫